Python scipy.misc.derivative() Examples

def calculate_engine_temperature_derivatives(
        times, engine_coolant_temperatures):
    """
    Calculates the derivative of the engine temperature [°C/s].

    :param times:
        Time vector [s].
    :type times: numpy.array

    :param engine_coolant_temperatures:
        Engine coolant temperature vector [°C].
    :type engine_coolant_temperatures: numpy.array

    :return:
        Derivative of the engine temperature [°C/s].
    :rtype: numpy.array
    """
    from statsmodels.nonparametric.smoothers_lowess import lowess
    par = dfl.functions.calculate_engine_temperature_derivatives
    temp = lowess(
        engine_coolant_temperatures, times, is_sorted=True,
        frac=par.tw * len(times) / (times[-1] - times[0]) ** 2, missing='none'
    )[:, 1].ravel()
    return _derivative(times, temp) 

def test_UNIFAC_misc():
    from scipy.misc import derivative
    from math import log
    T = 273.15 + 60
    
    def gE_T(T):
        xs = [0.5, 0.5]
        gammas = UNIFAC(chemgroups=[{1:2, 2:4}, {1:1, 2:1, 18:1}], T=T, xs=xs)
        return R*T*sum(xi*log(gamma) for xi, gamma in zip(xs, gammas))
    
    def hE_T(T):
        to_diff = lambda T: gE_T(T)/T
        return -derivative(to_diff, T,dx=1E-5, order=7)*T**2

    # A source gives 854.758 for hE, matching to within a gas constant
    assert_allclose(hE_T(T), 854.771631451345)
    assert_allclose(gE_T(T), 923.6408846044955) 

def test_PR78MIX():
    # Copied and pasted example from PR78.
    eos = PR78MIX(Tcs=[632], Pcs=[5350000], omegas=[0.734], zs=[1], T=299., P=1E6)
    three_props = [eos.V_l, eos.H_dep_l, eos.S_dep_l]
    expect_props = [8.351960066075052e-05, -63764.64948050847, -130.737108912626]
    assert_allclose(three_props, expect_props)

    # Fugacities
    eos = PR78MIX(T=115, P=1E6, Tcs=[126.1, 190.6], Pcs=[33.94E5, 46.04E5], omegas=[0.6, 0.7], zs=[0.5, 0.5], kijs=[[0,0],[0,0]])
    # Numerically test fugacities at one point, with artificially high omegas
    def numerical_fugacity_coefficient(n1, n2=0.5, switch=False, l=True):
        if switch:
            n1, n2 = n2, n1
        tot = n1+n2
        zs = [i/tot for i in [n1,n2]]
        a = PR78MIX(T=115, P=1E6, Tcs=[126.1, 190.6], Pcs=[33.94E5, 46.04E5], omegas=[0.6, 0.7], zs=zs, kijs=[[0,0],[0,0]])
        phi = a.phi_l if l else a.phi_g
        return tot*log(phi)

    phis = [[derivative(numerical_fugacity_coefficient, 0.5, dx=1E-6, order=25, args=(0.5, i, j)) for i in [False, True]] for j in [False, True]]
    assert_allclose(phis, [eos.lnphis_g, eos.lnphis_l]) 

def calculate_derivative(self, T, method, order=1):
        r'''Method to calculate a derivative of a property with respect to 
        temperature, of a given order  using a specified method. Uses SciPy's 
        derivative function, with a delta of 1E-6 K and a number of points 
        equal to 2*order + 1.

        This method can be overwritten by subclasses who may perfer to add
        analytical methods for some or all methods as this is much faster.

        If the calculation does not succeed, returns the actual error
        encountered.

        Parameters
        ----------
        T : float
            Temperature at which to calculate the derivative, [K]
        method : str
            Method for which to find the derivative
        order : int
            Order of the derivative, >= 1

        Returns
        -------
        derivative : float
            Calculated derivative property, [`units/K^order`]
        '''
        return derivative(self.calculate, T, dx=1e-6, args=[method], n=order, order=1+order*2) 

def calculate_derivative_P(self, P, T, method, order=1):
        r'''Method to calculate a derivative of a temperature and pressure
        dependent property with respect to pressure at constant temperature,
        of a given order using a specified method. Uses SciPy's derivative 
        function, with a delta of 0.01 Pa and a number of points equal to 
        2*order + 1.

        This method can be overwritten by subclasses who may perfer to add
        analytical methods for some or all methods as this is much faster.

        If the calculation does not succeed, returns the actual error
        encountered.

        Parameters
        ----------
        P : float
            Pressure at which to calculate the derivative, [Pa]
        T : float
            Temperature at which to calculate the derivative, [K]
        method : str
            Method for which to find the derivative
        order : int
            Order of the derivative, >= 1

        Returns
        -------
        d_prop_d_P_at_T : float
            Calculated derivative property at constant temperature, 
            [`units/Pa^order`]
        '''
        f = lambda P: self.calculate_P(T, P, method)
        return derivative(f, P, dx=1e-2, n=order, order=1+order*2) 

def T_dependent_property_derivative(self, T, order=1):
        r'''Method to obtain a derivative of a property with respect to 
        temperature, of a given order. Methods found valid by 
        `select_valid_methods` are attempted until a method succeeds. If no 
        methods are valid and succeed, None is returned.

        Calls `calculate_derivative` internally to perform the actual
        calculation.
        
        .. math::
            \text{derivative} = \frac{d (\text{property})}{d T}

        Parameters
        ----------
        T : float
            Temperature at which to calculate the derivative, [K]
        order : int
            Order of the derivative, >= 1

        Returns
        -------
        derivative : float
            Calculated derivative property, [`units/K^order`]
        '''
        if self.method:
            # retest within range
            if self.test_method_validity(T, self.method):
                try:
                    return self.calculate_derivative(T, self.method, order)
                except:  # pragma: no cover
                    pass
        sorted_valid_methods = self.select_valid_methods(T)
        for method in sorted_valid_methods:
            try:
                return self.calculate_derivative(T, method, order)
            except:
                pass
        return None 

def interval(self, alpha, *args, **kwds):
        """
        Confidence interval with equal areas around the median.

        Parameters
        ----------
        alpha : array_like of float
            Probability that an rv will be drawn from the returned range.
            Each value should be in the range [0, 1].
        arg1, arg2, ... : array_like
            The shape parameter(s) for the distribution (see docstring of the
            instance object for more information).
        loc : array_like, optional
            location parameter, Default is 0.
        scale : array_like, optional
            scale parameter, Default is 1.

        Returns
        -------
        a, b : ndarray of float
            end-points of range that contain ``100 * alpha %`` of the rv's
            possible values.

        """
        alpha = asarray(alpha)
        if np.any((alpha > 1) | (alpha < 0)):
            raise ValueError("alpha must be between 0 and 1 inclusive")
        q1 = (1.0-alpha)/2
        q2 = (1.0+alpha)/2
        a = self.ppf(q1, *args, **kwds)
        b = self.ppf(q2, *args, **kwds)
        return a, b


##  continuous random variables: implement maybe later
##
##  hf  --- Hazard Function (PDF / SF)
##  chf  --- Cumulative hazard function (-log(SF))
##  psf --- Probability sparsity function (reciprocal of the pdf) in
##                units of percent-point-function (as a function of q).
##                Also, the derivative of the percent-point function. 

def TP_dependent_property_derivative_T(self, T, P, order=1):
        r'''Method to calculate a derivative of a temperature and pressure
        dependent property with respect to temperature at constant pressure,
        of a given order. Methods found valid by `select_valid_methods_P` are 
        attempted until a method succeeds. If no methods are valid and succeed,
        None is returned.

        Calls `calculate_derivative_T` internally to perform the actual
        calculation.
        
        .. math::
            \text{derivative} = \frac{d (\text{property})}{d T}|_{P}

        Parameters
        ----------
        T : float
            Temperature at which to calculate the derivative, [K]
        P : float
            Pressure at which to calculate the derivative, [Pa]
        order : int
            Order of the derivative, >= 1

        Returns
        -------
        d_prop_d_T_at_P : float
            Calculated derivative property, [`units/K^order`]
        '''
        sorted_valid_methods_P = self.select_valid_methods_P(T, P)
        for method in sorted_valid_methods_P:
            try:
                return self.calculate_derivative_T(T, P, method, order)
            except:
                pass
        return None 

def TP_dependent_property_derivative_P(self, T, P, order=1):
        r'''Method to calculate a derivative of a temperature and pressure
        dependent property with respect to pressure at constant temperature,
        of a given order. Methods found valid by `select_valid_methods_P` are 
        attempted until a method succeeds. If no methods are valid and succeed,
        None is returned.

        Calls `calculate_derivative_P` internally to perform the actual
        calculation.
        
        .. math::
            \text{derivative} = \frac{d (\text{property})}{d P}|_{T}

        Parameters
        ----------
        T : float
            Temperature at which to calculate the derivative, [K]
        P : float
            Pressure at which to calculate the derivative, [Pa]
        order : int
            Order of the derivative, >= 1

        Returns
        -------
        d_prop_d_P_at_T : float
            Calculated derivative property, [`units/Pa^order`]
        '''
        sorted_valid_methods_P = self.select_valid_methods_P(T, P)
        for method in sorted_valid_methods_P:
            try:
                return self.calculate_derivative_P(P, T, method, order)
            except:
                pass
        return None 

def calculate_derivative_T(self, T, P, zs, ws, method, order=1):
        r'''Method to calculate a derivative of a mixture property with respect 
        to temperature at constant pressure and composition
        of a given order using a specified  method. Uses SciPy's derivative 
        function, with a delta of 1E-6 K and a number of points equal to 
        2*order + 1.

        This method can be overwritten by subclasses who may perfer to add
        analytical methods for some or all methods as this is much faster.

        If the calculation does not succeed, returns the actual error
        encountered.

        Parameters
        ----------
        T : float
            Temperature at which to calculate the derivative, [K]
        P : float
            Pressure at which to calculate the derivative, [Pa]
        zs : list[float]
            Mole fractions of all species in the mixture, [-]
        ws : list[float]
            Weight fractions of all species in the mixture, [-]
        method : str
            Method for which to find the derivative
        order : int
            Order of the derivative, >= 1

        Returns
        -------
        d_prop_d_T_at_P : float
            Calculated derivative property at constant pressure, 
            [`units/K^order`]
        '''
        return derivative(self.calculate, T, dx=1e-6, args=[P, zs, ws, method], n=order, order=1+order*2) 

def calculate_derivative_P(self, P, T, zs, ws, method, order=1):
        r'''Method to calculate a derivative of a mixture property with respect 
        to pressure at constant temperature and composition
        of a given order using a specified method. Uses SciPy's derivative 
        function, with a delta of 0.01 Pa and a number of points equal to 
        2*order + 1.

        This method can be overwritten by subclasses who may perfer to add
        analytical methods for some or all methods as this is much faster.

        If the calculation does not succeed, returns the actual error
        encountered.

        Parameters
        ----------
        P : float
            Pressure at which to calculate the derivative, [Pa]
        T : float
            Temperature at which to calculate the derivative, [K]
        zs : list[float]
            Mole fractions of all species in the mixture, [-]
        ws : list[float]
            Weight fractions of all species in the mixture, [-]
        method : str
            Method for which to find the derivative
        order : int
            Order of the derivative, >= 1

        Returns
        -------
        d_prop_d_P_at_T : float
            Calculated derivative property at constant temperature, 
            [`units/Pa^order`]
        '''
        f = lambda P: self.calculate(T, P, zs, ws, method)
        return derivative(f, P, dx=1e-2, n=order, order=1+order*2) 

def property_derivative_P(self, T, P, zs, ws, order=1):
        r'''Method to calculate a derivative of a mixture property with respect
        to pressure at constant temperature and composition,
        of a given order. Methods found valid by `select_valid_methods` are 
        attempted until a method succeeds. If no methods are valid and succeed,
        None is returned.

        Calls `calculate_derivative_P` internally to perform the actual
        calculation.
        
        .. math::
            \text{derivative} = \frac{d (\text{property})}{d P}|_{T, z}

        Parameters
        ----------
        T : float
            Temperature at which to calculate the derivative, [K]
        P : float
            Pressure at which to calculate the derivative, [Pa]
        zs : list[float]
            Mole fractions of all species in the mixture, [-]
        ws : list[float]
            Weight fractions of all species in the mixture, [-]
        order : int
            Order of the derivative, >= 1

        Returns
        -------
        d_prop_d_P_at_T : float
            Calculated derivative property, [`units/Pa^order`]
        '''
        sorted_valid_methods = self.select_valid_methods(T, P, zs, ws)
        for method in sorted_valid_methods:
            try:
                return self.calculate_derivative_P(P, T, zs, ws, method, order)
            except:
                pass
        return None 

def fugacities_partial_derivatives(self, xs=None, ys=None):
        if self.phase in ['l', 'l/g']:
            if xs is None:
                xs = self.zs
            self.dphis_dni_l = [derivative(self._dphi_dn, xs[i], args=[i, 'l'], dx=1E-7, n=1) for i in self.cmps]
            self.dfugacities_dni_l = [derivative(self._dfugacity_dn, xs[i], args=[i, 'l'], dx=1E-7, n=1) for i in self.cmps]
            self.dlnphis_dni_l = [dphi/phi for dphi, phi in zip(self.dphis_dni_l, self.phis_l)]
        if self.phase in ['g', 'l/g']:
            if ys is None:
                ys = self.zs
            self.dphis_dni_g = [derivative(self._dphi_dn, ys[i], args=[i, 'g'], dx=1E-7, n=1) for i in self.cmps]
            self.dfugacities_dni_g = [derivative(self._dfugacity_dn, ys[i], args=[i, 'g'], dx=1E-7, n=1) for i in self.cmps]
            self.dlnphis_dni_g = [dphi/phi for dphi, phi in zip(self.dphis_dni_g, self.phis_g)]
            # confirmed the relationship of the above 
            # There should be an easy way to get dfugacities_dn_g but I haven't figured it out 

def fugacities_partial_derivatives_2(self, xs=None, ys=None):
        if self.phase in ['l', 'l/g']:
            if xs is None:
                xs = self.zs
            self.d2phis_dni2_l = [derivative(self._dphi_dn, xs[i], args=[i, 'l'], dx=1E-5, n=2) for i in self.cmps]
            self.d2fugacities_dni2_l = [derivative(self._dfugacity_dn, xs[i], args=[i, 'l'], dx=1E-5, n=2) for i in self.cmps]
            self.d2lnphis_dni2_l = [d2phi/phi  - dphi*dphi/(phi*phi) for d2phi, dphi, phi in zip(self.d2phis_dni2_l, self.dphis_dni_l, self.phis_l)]
        if self.phase in ['g', 'l/g']:
            if ys is None:
                ys = self.zs
            self.d2phis_dni2_g = [derivative(self._dphi_dn, ys[i], args=[i, 'g'], dx=1E-5, n=2) for i in self.cmps]
            self.d2fugacities_dni2_g = [derivative(self._dfugacity_dn, ys[i], args=[i, 'g'], dx=1E-5, n=2) for i in self.cmps]
            self.d2lnphis_dni2_g = [d2phi/phi  - dphi*dphi/(phi*phi) for d2phi, dphi, phi in zip(self.d2phis_dni2_g, self.dphis_dni_g, self.phis_g)]
        # second derivative lns confirmed 

def test_EQ127_more():
    # T derivative
    coeffs = (3.3258E4, 3.6199E4, 1.2057E3, 1.5373E7, 3.2122E3, -1.5318E7, 3.2122E3)
    diff_1T = derivative(EQ127, 50,  dx=1E-3, order=21, args=coeffs)
    diff_1T_analytical = EQ127(50., *coeffs, order=1)
    assert_allclose(diff_1T, diff_1T_analytical, rtol=1E-3)
    assert_allclose(diff_1T, 0.000313581049006, rtol=1E-4)
    
    # Integral
    int_50 = EQ127(50., *coeffs, order=-1) 
    int_20 = EQ127(20., *coeffs, order=-1)
    numerical_1T = quad(EQ127, 20, 50, args=coeffs)[0]
    assert_allclose(int_50 - int_20, numerical_1T)
    assert_allclose(numerical_1T, 997740.00147014)
    
    # Integral over T
    T_int_50 = EQ127(50., *coeffs, order=-1j)
    T_int_20 = EQ127(20., *coeffs, order=-1j)
    
    to_int = lambda T :EQ127(T, *coeffs)/T
    numerical_1_over_T = quad(to_int, 20, 50)[0]
    assert_allclose(T_int_50 - T_int_20, numerical_1_over_T)
    assert_allclose(T_int_50 - T_int_20, 30473.9971912935)

    with pytest.raises(Exception):
        EQ127(20., *coeffs, order=1E100) 

def test_EQ116_more():
    # T derivative
    coeffs = (647.096, 17.863, 58.606, -95.396, 213.89, -141.26)
    diff_1T = derivative(EQ116, 50,  dx=1E-3, order=21, args=coeffs)
    diff_1T_analytical = EQ116(50., *coeffs, order=1)
    assert_allclose(diff_1T, diff_1T_analytical, rtol=1E-3)
    assert_allclose(diff_1T_analytical, 0.020379262711650914)
    
    # Integral
    int_50 = EQ116(50., *coeffs, order=-1) 
    int_20 = EQ116(20., *coeffs, order=-1)
    numerical_1T = quad(EQ116, 20, 50, args=coeffs)[0]
    assert_allclose(int_50 - int_20, numerical_1T)
    assert_allclose(int_50 - int_20, 1636.962423782701)
    
    # Integral over T
    T_int_50 = EQ116(50., *coeffs, order=-1j)
    T_int_20 = EQ116(20., *coeffs, order=-1j)
    
    to_int = lambda T :EQ116(T, *coeffs)/T
    numerical_1_over_T = quad(to_int, 20, 50)[0]
    assert_allclose(T_int_50 - T_int_20, numerical_1_over_T)
    assert_allclose(T_int_50 - T_int_20, 49.95109104018752)
    
    with pytest.raises(Exception):
        EQ116(20., *coeffs, order=1E100) 

def test_EQ114_more():
    # T derivative
    coeffs = (33.19, 66.653, 6765.9, -123.63, 478.27)
    diff_1T = derivative(EQ114, 20,  dx=1E-3, order=21, args=coeffs)
    diff_1T_analytical = EQ114(20., *coeffs, order=1)
    assert_allclose(diff_1T, diff_1T_analytical, rtol=1E-3)
    assert_allclose(diff_1T, 1135.38618941)
    
    # Integral
    int_50 = EQ114(30., *coeffs, order=-1) 
    int_20 = EQ114(20., *coeffs, order=-1)
    numerical_1T = quad(EQ114, 20, 30, args=coeffs)[0]
    assert_allclose(int_50 - int_20, numerical_1T)
    assert_allclose(int_50 - int_20, 295697.48978888744)
    
#     Integral over T
    T_int_50 = EQ114(30., *coeffs, order=-1j)
    T_int_20 = EQ114(20., *coeffs, order=-1j)
    
    to_int = lambda T :EQ114(T, *coeffs)/T
    numerical_1_over_T = quad(to_int, 20, 30)[0]
    assert_allclose(T_int_50 - T_int_20, numerical_1_over_T)
    assert_allclose(T_int_50 - T_int_20, 11612.331762721366)
    
    with pytest.raises(Exception):
        EQ114(20., *coeffs, order=1E100) 

def test_EQ102_more():
    # T derivative
    coeffs = (1.7096E-8, 1.1146, 1.1, 2.1)
    diff_1T = derivative(EQ102, 250,  dx=1E-3, order=21, args=coeffs)
    diff_1T_analytical = EQ102(250., *coeffs, order=1)
    assert_allclose(diff_1T, diff_1T_analytical, rtol=1E-3)
    assert_allclose(diff_1T, 3.5861274167602139e-08)
    
    # Integral
    int_250 = EQ102(250., *coeffs, order=-1) 
    int_220 = EQ102(220., *coeffs, order=-1)
    numerical_1T = quad(EQ102, 220, 250, args=coeffs)[0]
    assert_allclose(int_250 - int_220, numerical_1T)
    assert_allclose(int_250 - int_220, 0.00022428562125110119)
    
#     Integral over T
    T_int_250 = EQ102(250., *coeffs, order=-1j)
    T_int_220 = EQ102(220., *coeffs, order=-1j)
    
    to_int = lambda T :EQ102(T, *coeffs)/T
    numerical_1_over_T = quad(to_int, 220, 250)[0]
    assert_allclose(T_int_250 - T_int_220, numerical_1_over_T)
    assert_allclose(T_int_250 - T_int_220, 9.5425212178091671e-07)
#    
    with pytest.raises(Exception):
        EQ102(20., *coeffs, order=1E100) 

def test_EQ100_more():
    # T derivative
    coeffs = (276370., -2090.1, 8.125, -0.014116, 0.0000093701)
    diff_1T = derivative(EQ100, 250,  dx=1E-3, order=21, args=coeffs)
    diff_1T_analytical = EQ100(250., *coeffs, order=1)
    assert_allclose(diff_1T, diff_1T_analytical, rtol=1E-3)
    assert_allclose(diff_1T, -88.7187500531)
    
    # Integral
    int_250 = EQ100(250., *coeffs, order=-1) 
    int_220 = EQ100(220., *coeffs, order=-1)
    numerical_1T = quad(EQ100, 220, 250, args=coeffs)[0]
    assert_allclose(int_250 - int_220, numerical_1T)
    assert_allclose(int_250 - int_220, 2381304.7021859996)
    
#     Integral over T
    T_int_250 = EQ100(250., *coeffs, order=-1j)
    T_int_220 = EQ100(220., *coeffs, order=-1j)
    
    to_int = lambda T :EQ100(T, *coeffs)/T
    numerical_1_over_T = quad(to_int, 220, 250)[0]
    assert_allclose(T_int_250 - T_int_220, numerical_1_over_T)
    assert_allclose(T_int_250 - T_int_220, 10152.09780143667)

    with pytest.raises(Exception):
        EQ100(20., *coeffs, order=1E100) 

def test_EQ104_more():
    # T derivative
    coeffs = (0.02222, -26.38, -16750000, -3.894E19, 3.133E21)
    diff_1T = derivative(EQ104, 250,  dx=1E-3, order=21, args=coeffs)
    diff_1T_analytical = EQ104(250., *coeffs, order=1)
    assert_allclose(diff_1T, diff_1T_analytical, rtol=1E-3)
    assert_allclose(diff_1T, 0.0653824814073)
    
    # Integral
    int_250 = EQ104(250., *coeffs, order=-1) 
    int_220 = EQ104(220., *coeffs, order=-1)
    numerical_1T = quad(EQ104, 220, 250, args=coeffs)[0]
    assert_allclose(int_250 - int_220, numerical_1T)
    assert_allclose(int_250 - int_220, -127.91851427119406)
    
#     Integral over T
    T_int_250 = EQ104(250., *coeffs, order=-1j)
    T_int_220 = EQ104(220., *coeffs, order=-1j)
    
    to_int = lambda T :EQ104(T, *coeffs)/T
    numerical_1_over_T = quad(to_int, 220, 250)[0]
    assert_allclose(T_int_250 - T_int_220, numerical_1_over_T)
    assert_allclose(T_int_250 - T_int_220, -0.5494851210308727)

    with pytest.raises(Exception):
        EQ104(20., *coeffs, order=1E100) 

def test_fuzz_dV_dT_and_d2V_dT2_derivatives():
    from thermo import eos
    eos_list = list(eos.__all__); eos_list.remove('GCEOS')
    eos_list.remove('ALPHA_FUNCTIONS'); eos_list.remove('VDW')
    
    phase_extensions = {True: '_l', False: '_g'}
    derivative_bases_dV_dT = {0:'V', 1:'dV_dT', 2:'d2V_dT2'}
    
    def dV_dT(T, P, eos, order=0, phase=True, Tc=507.6, Pc=3025000., omega=0.2975):
        eos = globals()[eos_list[eos]](Tc=Tc, Pc=Pc, omega=omega, T=T, P=P)
        phase_base = phase_extensions[phase]
        attr = derivative_bases_dV_dT[order]+phase_base
        return getattr(eos, attr)
    
    x, y = [], []
    for eos in range(len(eos_list)):
        for T in np.linspace(.1, 1000, 50):
            for P in np.logspace(np.log10(3E4), np.log10(1E6), 50):
                T, P = float(T), float(P)
                for phase in [True, False]:
                    for order in [1, 2]:
                        try:
                            # If dV_dx_phase doesn't exist, will simply abort and continue the loop
                            numer = derivative(dV_dT, T, dx=1E-4, args=(P, eos, order-1, phase))
                            ana = dV_dT(T=T, P=P, eos=eos, order=order, phase=phase)
                        except:
                            continue
                        x.append(numer)
                        y.append(ana)
    assert allclose_variable(x, y, limits=[.009, .05, .65, .93],rtols=[1E-5, 1E-6, 1E-9, 1E-10]) 

def test_fuzz_dV_dP_and_d2V_dP2_derivatives():
    from thermo import eos
    eos_list = list(eos.__all__); eos_list.remove('GCEOS')
    eos_list.remove('ALPHA_FUNCTIONS'); eos_list.remove('VDW')
    
    phase_extensions = {True: '_l', False: '_g'}
    derivative_bases_dV_dP = {0:'V', 1:'dV_dP', 2:'d2V_dP2'}
    
    def dV_dP(P, T, eos, order=0, phase=True, Tc=507.6, Pc=3025000., omega=0.2975):
        eos = globals()[eos_list[eos]](Tc=Tc, Pc=Pc, omega=omega, T=T, P=P)
        phase_base = phase_extensions[phase]
        attr = derivative_bases_dV_dP[order]+phase_base
        return getattr(eos, attr)
    
    
    x, y = [], []
    for eos in range(len(eos_list)):
        for T in np.linspace(.1, 1000, 50):
            for P in np.logspace(np.log10(3E4), np.log10(1E6), 50):
                T, P = float(T), float(P)
                for phase in [True, False]:
                    for order in [1, 2]:
                        try:
                            # If dV_dx_phase doesn't exist, will simply abort and continue the loop
                            numer = derivative(dV_dP, P, dx=15., args=(T, eos, order-1, phase))
                            ana = dV_dP(T=T, P=P, eos=eos, order=order, phase=phase)
                        except:
                            continue
                        x.append(numer)
                        y.append(ana)
    assert allclose_variable(x, y, limits=[.02, .04, .04, .05, .15, .45, .95],
                            rtols=[1E-2, 1E-3, 1E-4, 1E-5, 1E-6, 1E-7, 1E-9]) 

def test_fuzz_dPsat_dT():
    from thermo import eos
    eos_list = list(eos.__all__); eos_list.remove('GCEOS')
    eos_list.remove('ALPHA_FUNCTIONS'); eos_list.remove('eos_list')
    eos_list.remove('GCEOS_DUMMY')
    
    Tc = 507.6
    Pc = 3025000
    omega = 0.2975
    
    e = PR(T=400, P=1E5, Tc=507.6, Pc=3025000, omega=0.2975)
    dPsats_dT_expect = [938.7777925283981, 10287.225576267781, 38814.74676693623]
    assert_allclose([e.dPsat_dT(300), e.dPsat_dT(400), e.dPsat_dT(500)], dPsats_dT_expect)
    
    # Hammer the derivatives for each EOS in a wide range; most are really 
    # accurate. There's an error around the transition between polynomials 
    # though - to be expected; the derivatives are discontinuous there.
    dPsats_derivative = []
    dPsats_analytical = []
    for eos in range(len(eos_list)):
        for T in np.linspace(0.2*Tc, Tc*.999, 50):
            e = globals()[eos_list[eos]](Tc=Tc, Pc=Pc, omega=omega, T=T, P=1E5)
            anal = e.dPsat_dT(T)
            numer = derivative(e.Psat, T, order=9)
            dPsats_analytical.append(anal)
            dPsats_derivative.append(numer)
    assert allclose_variable(dPsats_derivative, dPsats_analytical, limits=[.02, .06], rtols=[1E-5, 1E-7]) 

def test_exponential_integral(x, y):
    z = x + 1j*y

    # Compare with Scipy implementation
    assert np.isclose(E1(z), exp1(z), rtol=1e-3)

    # Test property (A3.5) of the function according to [Del, p.367]. 
    if y != 0.0:
        assert np.isclose(E1(np.conjugate(z)), np.conjugate(E1(z)), rtol=1e-3)

    # BROKEN...
    # def derivative_of_complex_function(f, z, **kwargs):
    #     direction = 1
    #     return derivative(lambda eps: f(z + eps*direction), 0.0, **kwargs)
    # if abs(x) > 1 and abs(y) > 1:
    #     assert np.isclose(
    #         derivative_of_complex_function(Delhommeau_f90.initialize_green_wave.gg, z),
    #         Delhommeau_f90.initialize_green_wave.gg(z) - 1.0/z,
    #         atol=1e-2)


# CHECK OF THE THEORY, NOT THE CODE ITSELF
# Not necessary to run every time...
#
# @given(r=floats(min_value=0, max_value=1e2),
#        z=floats(min_value=-16, max_value=-1),
#        k=floats(min_value=0.1, max_value=10)
#        )
# def test_zeta(r, z, k):
#     """Check expression k/π ∫ 1/ζ(θ) dθ = - 1/r """

#     def one_over_zeta(theta):
#         return np.real(1/(k*(z + 1j*r*np.cos(theta))))

#     int_one_over_zeta, *_ = quad(one_over_zeta, -pi/2, pi/2)

#     assert np.isclose(k/pi * int_one_over_zeta,
#                       -1/(np.sqrt(r**2 + z**2)),
#                       rtol=1e-4) 

def _derivative(times, temp):
    import scipy.misc as sci_misc
    import scipy.interpolate as sci_itp
    par = dfl.functions.calculate_engine_temperature_derivatives
    func = sci_itp.InterpolatedUnivariateSpline(times, temp, k=1)
    return sci_misc.derivative(func, times, dx=par.dx, order=par.order) 

def isothermal_compressibility(V, dV_dP):
    r'''Calculate the isothermal coefficient of a compressibility, given its 
    molar volume at a certain `T` and `P`, and its derivative of molar volume
    with respect to `P`.

    .. math::
        \kappa = -\frac{1}{V}\left(\frac{\partial V}{\partial P} \right)_T

    Parameters
    ----------
    V : float
        Molar volume at `T` and `P`, [m^3/mol]
    dV_dP : float
        Derivative of molar volume with respect to `P`, [m^3/mol/Pa]

    Returns
    -------
    kappa : float
        Isothermal coefficient of a compressibility, [1/Pa]
        
    Notes
    -----
    For an ideal gas, this expression simplified to:
    
    .. math::
        \kappa = \frac{1}{P}

    Examples
    --------
    Calculated for hexane from the PR EOS at 299 K and 1 MPa (liquid):
    
    >>> isothermal_compressibility(0.000130229900873546, -2.72902118209903e-13)
    2.095541165119158e-09

    References
    ----------
    .. [1] Poling, Bruce E. The Properties of Gases and Liquids. 5th edition.
       New York: McGraw-Hill Professional, 2000.
    '''
    return -dV_dP/V 

def _pdf(self, x, *args):
        return derivative(self._cdf, x, dx=1e-5, args=args, order=5)

    ## Could also define any of these 

def norm_lls_grad(y, params):
    '''Jacobian of normal loglikelihood wrt mean mu and variance sigma2

    Parameters
    ----------
    y : array, 1d
        normally distributed random variable
    params: array, (nobs, 2)
        array of mean, variance (mu, sigma2) with observations in rows

    Returns
    -------
    grad : array (nobs, 2)
        derivative of loglikelihood for each observation wrt mean in first
        column, and wrt variance in second column

    Notes
    -----
    this is actually the derivative wrt sigma not sigma**2, but evaluated
    with parameter sigma2 = sigma**2

    '''
    mu, sigma2 = params.T
    dllsdmu = (y-mu)/sigma2
    dllsdsigma2 = ((y-mu)**2/sigma2 - 1)/np.sqrt(sigma2)
    return np.column_stack((dllsdmu, dllsdsigma2)) 

def normgrad(y, x, params):
    '''Jacobian of normal loglikelihood wrt mean mu and variance sigma2

    Parameters
    ----------
    y : array, 1d
        normally distributed random variable with mean x*beta, and variance sigma2
    x : array, 2d
        explanatory variables, observation in rows, variables in columns
    params: array_like, (nvars + 1)
        array of coefficients and variance (beta, sigma2)

    Returns
    -------
    grad : array (nobs, 2)
        derivative of loglikelihood for each observation wrt mean in first
        column, and wrt scale (sigma) in second column
    assume params = (beta, sigma2)

    Notes
    -----
    TODO: for heteroscedasticity need sigma to be a 1d array

    '''
    beta = params[:-1]
    sigma2 = params[-1]*np.ones((len(y),1))
    dmudbeta = mean_grad(x, beta)
    mu = np.dot(x, beta)
    #print(beta, sigma2)
    params2 = np.column_stack((mu,sigma2))
    dllsdms = norm_lls_grad(y,params2)
    grad = np.column_stack((dllsdms[:,:1]*dmudbeta, dllsdms[:,:1]))
    return grad 

def norm_dlldy(y):
    '''derivative of log pdf of standard normal with respect to y
    '''
    return -y