Python scipy.integrate.quad() Examples

def test_roots_chebys():
    weightf = orth.chebys(5).weight_func
    verify_gauss_quad(sc.roots_chebys, orth.eval_chebys, weightf, -2., 2., 5)
    verify_gauss_quad(sc.roots_chebys, orth.eval_chebys, weightf, -2., 2., 25)
    verify_gauss_quad(sc.roots_chebys, orth.eval_chebys, weightf, -2., 2., 100)

    x, w = sc.roots_chebys(5, False)
    y, v, m = sc.roots_chebys(5, True)
    assert_allclose(x, y, 1e-14, 1e-14)
    assert_allclose(w, v, 1e-14, 1e-14)

    muI, muI_err = integrate.quad(weightf, -2, 2)
    assert_allclose(m, muI, rtol=muI_err)

    assert_raises(ValueError, sc.roots_chebys, 0)
    assert_raises(ValueError, sc.roots_chebys, 3.3) 

def get_heat_capacity(temp, td):
        # credits to Dr. Joseph Montoya, LBNL
        t_ratio = temp / td

        def integrand(x):
            return (x ** 4 * pd.np.exp(x)) / (pd.np.exp(x) - 1) ** 2

        if isinstance(t_ratio, int) or isinstance(t_ratio, float):
            cv_p = 9 * R * (t_ratio ** 3) * quad(integrand, 0, t_ratio ** -1)[0]
        else:
            cv_p = []
            for i in range(len(t_ratio)):
                cv_i = (
                    9 * R * (t_ratio[i] ** 3) * quad(integrand, 0, t_ratio[i] ** -1)[0]
                )
                cv_p = pd.np.append(cv_p, cv_i)
        return cv_p * 5 

def _entropy(self, *args):
        def integ(x):
            val = self._pdf(x, *args)
            return entr(val)

        # upper limit is often inf, so suppress warnings when integrating
        olderr = np.seterr(over='ignore')
        h = integrate.quad(integ, self.a, self.b)[0]
        np.seterr(**olderr)

        if not np.isnan(h):
            return h
        else:
            # try with different limits if integration problems
            low, upp = self.ppf([1e-10, 1. - 1e-10], *args)
            if np.isinf(self.b):
                upper = upp
            else:
                upper = self.b
            if np.isinf(self.a):
                lower = low
            else:
                lower = self.a
            return integrate.quad(integ, lower, upper)[0] 

def mvstdtprob(a, b, R, df, ieps=1e-5, quadkwds=None, mvstkwds=None):
    '''probability of rectangular area of standard t distribution

    assumes mean is zero and R is correlation matrix

    Notes
    -----
    This function does not calculate the estimate of the combined error
    between the underlying multivariate normal probability calculations
    and the integration.

    '''
    kwds = dict(args=(a,b,R,df), epsabs=1e-4, epsrel=1e-2, limit=150)
    if not quadkwds is None:
        kwds.update(quadkwds)
    #print kwds
    res, err = integrate.quad(funbgh2, *chi.ppf([ieps,1-ieps], df),
                          **kwds)
    prob = res * bghfactor(df)
    return prob

#written by Enzo Michelangeli, style changes by josef-pktd
# Student's T random variable 

def cdf(self):
        """
        Returns the cumulative distribution function evaluated at the support.

        Notes
        -----
        Will not work if fit has not been called.
        """
        _checkisfit(self)
        kern = self.kernel
        if kern.domain is None: # TODO: test for grid point at domain bound
            a,b = -np.inf,np.inf
        else:
            a,b = kern.domain
        func = lambda x,s: kern.density(s,x)

        support = self.support
        support = np.r_[a,support]
        gridsize = len(support)
        endog = self.endog
        probs = [integrate.quad(func, support[i - 1], support[i],
                    args=endog)[0] for i in range(1, gridsize)]
        return np.cumsum(probs) 

def entropy(self):
        """
        Returns the differential entropy evaluated at the support

        Notes
        -----
        Will not work if fit has not been called. 1e-12 is added to each
        probability to ensure that log(0) is not called.
        """
        _checkisfit(self)

        def entr(x,s):
            pdf = kern.density(s,x)
            return pdf*np.log(pdf+1e-12)

        kern = self.kernel

        if kern.domain is not None:
            a, b = self.domain
        else:
            a, b = -np.inf, np.inf
        endog = self.endog
        #TODO: below could run into integr problems, cf. stats.dist._entropy
        return -integrate.quad(entr, a, b, args=(endog,))[0] 

def length(self, t0=0, t1=1, error=LENGTH_ERROR, min_depth=LENGTH_MIN_DEPTH):
        """Calculate the length of the path up to a certain position"""
        if t0 == 0 and t1 == 1:
            if self._length_info['bpoints'] == self.bpoints() \
                    and self._length_info['error'] >= error \
                    and self._length_info['min_depth'] >= min_depth:
                return self._length_info['length']

        # using scipy.integrate.quad is quick
        if _quad_available:
            s = quad(lambda tau: abs(self.derivative(tau)), t0, t1,
                            epsabs=error, limit=1000)[0]
        else:
            s = segment_length(self, t0, t1, self.point(t0), self.point(t1),
                               error, min_depth, 0)

        if t0 == 0 and t1 == 1:
            self._length_info['length'] = s
            self._length_info['bpoints'] = self.bpoints()
            self._length_info['error'] = error
            self._length_info['min_depth'] = min_depth
            return self._length_info['length']
        else:
            return s 

def length(self, t0=0, t1=1, error=LENGTH_ERROR, min_depth=LENGTH_MIN_DEPTH):
        """The length of an elliptical large_arc segment requires numerical
        integration, and in that case it's simpler to just do a geometric
        approximation, as for cubic bezier curves."""
        assert 0 <= t0 <= 1 and 0 <= t1 <= 1

        if t0 == 0 and t1 == 1:
            h = hash(self)
            if self.segment_length_hash is None or self.segment_length_hash != h:
                self.segment_length_hash = h
                if _quad_available:
                    self.segment_length = quad(lambda tau: abs(self.derivative(tau)),
                                               t0, t1, epsabs=error, limit=1000)[0]
                else:
                    self.segment_length = segment_length(self, t0, t1, self.point(t0),
                                                         self.point(t1), error, min_depth, 0)
            return self.segment_length

        if _quad_available:
            return quad(lambda tau: abs(self.derivative(tau)), t0, t1,
                        epsabs=error, limit=1000)[0]
        else:
            return segment_length(self, t0, t1, self.point(t0), self.point(t1),
                                  error, min_depth, 0) 

def calculate_c_baseline(c_L, Au_C, Au_L, deltaz):
    """Equations in the New_CST.pdf. Calculates the upper chord in order for
       the cruise and landing airfoils ot have the same length."""

    def integrand(psi, Au, delta_xi):
        return np.sqrt(1 + dxi_u(psi, Au, delta_xi)**2)

    def f(c_C):
        """Function dependent of c_C and that outputs c_C."""
        y_C, err = quad(integrand, 0, 1, args=(Au_C, deltaz/c_C))
        y_L, err = quad(integrand, 0, 1, args=(Au_L, deltaz/c_L))
        return c_L*y_L/y_C
    c_C = optimize.fixed_point(f, [c_L])
    # In case the calculated chord is really close to the original, but the
    # algorithm was not able to make them equal
    if abs(c_L - c_C) < 1e-7:
        return c_L
    # The output is an array so it needs the extra [0]
    return c_C[0] 

def calculate_psi_goal(psi_baseline, Au_baseline, Au_goal, deltaz,
                       c_baseline, c_goal):
    """Find the value for psi that has the same location w on the upper
    surface of the goal as psi_baseline on the upper surface of the
    baseline"""

    def integrand(psi_baseline, Au, deltaz, c):
        return c*np.sqrt(1 + dxi_u(psi_baseline, Au, deltaz/c)**2)

    def equation(psi_goal, L_baseline, Au_goal, deltaz, c):
        y, err = quad(integrand, 0, psi_goal, args=(Au_goal, deltaz, c))
        return y - L_baseline

    L_baseline, err = quad(integrand, 0, psi_baseline,
                           args=(Au_baseline, deltaz, c_baseline))
    with warnings.catch_warnings():
        warnings.simplefilter("ignore")
        y = fsolve(equation, psi_baseline, args=(L_baseline, Au_goal, deltaz,
                                                 c_goal))
    return y[0] 

def verify_gauss_quad(root_func, eval_func, weight_func, a, b, N,
                      rtol=1e-15, atol=1e-14):
        # this test is copied from numpy's TestGauss in test_hermite.py
        x, w, mu = root_func(N, True)

        n = np.arange(N)
        v = eval_func(n[:,np.newaxis], x)
        vv = np.dot(v*w, v.T)
        vd = 1 / np.sqrt(vv.diagonal())
        vv = vd[:, np.newaxis] * vv * vd
        assert_allclose(vv, np.eye(N), rtol, atol)

        # check that the integral of 1 is correct
        assert_allclose(w.sum(), mu, rtol, atol)

        # compare the results of integrating a function with quad.
        f = lambda x: x**3 - 3*x**2 + x - 2
        resI = integrate.quad(lambda x: f(x)*weight_func(x), a, b)
        resG = np.vdot(f(x), w)
        rtol = 1e-6 if 1e-6 < resI[1] else resI[1] * 10
        assert_allclose(resI[0], resG, rtol=rtol) 

def test_roots_hermitenorm():
    rootf = sc.roots_hermitenorm
    evalf = orth.eval_hermitenorm
    weightf = orth.hermitenorm(5).weight_func

    verify_gauss_quad(rootf, evalf, weightf, -np.inf, np.inf, 5)
    verify_gauss_quad(rootf, evalf, weightf, -np.inf, np.inf, 25, atol=1e-13)
    verify_gauss_quad(rootf, evalf, weightf, -np.inf, np.inf, 100, atol=1e-12)

    x, w = sc.roots_hermitenorm(5, False)
    y, v, m = sc.roots_hermitenorm(5, True)
    assert_allclose(x, y, 1e-14, 1e-14)
    assert_allclose(w, v, 1e-14, 1e-14)

    muI, muI_err = integrate.quad(weightf, -np.inf, np.inf)
    assert_allclose(m, muI, rtol=muI_err)

    assert_raises(ValueError, sc.roots_hermitenorm, 0)
    assert_raises(ValueError, sc.roots_hermitenorm, 3.3) 

def test_roots_chebyt():
    weightf = orth.chebyt(5).weight_func
    verify_gauss_quad(sc.roots_chebyt, orth.eval_chebyt, weightf, -1., 1., 5)
    verify_gauss_quad(sc.roots_chebyt, orth.eval_chebyt, weightf, -1., 1., 25)
    verify_gauss_quad(sc.roots_chebyt, orth.eval_chebyt, weightf, -1., 1., 100, atol=1e-12)

    x, w = sc.roots_chebyt(5, False)
    y, v, m = sc.roots_chebyt(5, True)
    assert_allclose(x, y, 1e-14, 1e-14)
    assert_allclose(w, v, 1e-14, 1e-14)

    muI, muI_err = integrate.quad(weightf, -1, 1)
    assert_allclose(m, muI, rtol=muI_err)

    assert_raises(ValueError, sc.roots_chebyt, 0)
    assert_raises(ValueError, sc.roots_chebyt, 3.3) 

def test_roots_chebyu():
    weightf = orth.chebyu(5).weight_func
    verify_gauss_quad(sc.roots_chebyu, orth.eval_chebyu, weightf, -1., 1., 5)
    verify_gauss_quad(sc.roots_chebyu, orth.eval_chebyu, weightf, -1., 1., 25)
    verify_gauss_quad(sc.roots_chebyu, orth.eval_chebyu, weightf, -1., 1., 100)

    x, w = sc.roots_chebyu(5, False)
    y, v, m = sc.roots_chebyu(5, True)
    assert_allclose(x, y, 1e-14, 1e-14)
    assert_allclose(w, v, 1e-14, 1e-14)

    muI, muI_err = integrate.quad(weightf, -1, 1)
    assert_allclose(m, muI, rtol=muI_err)

    assert_raises(ValueError, sc.roots_chebyu, 0)
    assert_raises(ValueError, sc.roots_chebyu, 3.3) 

def test_roots_chebyc():
    weightf = orth.chebyc(5).weight_func
    verify_gauss_quad(sc.roots_chebyc, orth.eval_chebyc, weightf, -2., 2., 5)
    verify_gauss_quad(sc.roots_chebyc, orth.eval_chebyc, weightf, -2., 2., 25)
    verify_gauss_quad(sc.roots_chebyc, orth.eval_chebyc, weightf, -2., 2., 100, atol=1e-12)

    x, w = sc.roots_chebyc(5, False)
    y, v, m = sc.roots_chebyc(5, True)
    assert_allclose(x, y, 1e-14, 1e-14)
    assert_allclose(w, v, 1e-14, 1e-14)

    muI, muI_err = integrate.quad(weightf, -2, 2)
    assert_allclose(m, muI, rtol=muI_err)

    assert_raises(ValueError, sc.roots_chebyc, 0)
    assert_raises(ValueError, sc.roots_chebyc, 3.3) 

def integral_bounded(fn, lb, ub):
  integral, _ = integrate.quad(fn, lb, ub)
  return integral 

def find_chord(P):
    def _to_minimize(l):
        def _to_integrate(y):
            den = np.sqrt(1-P**2/E**2/I**2*(l*y-y**2/2)**2)
            if np.isnan(den):
                return 100
            else:
                return 1/den
        l = l[0]
        current_L = quad(_to_integrate, 0, l)[0]
        return abs(L-current_L)
        
    return minimize(_to_minimize, L, method = 'Nelder-Mead',).x[0] 

def find_deflection(y, l, P):
    def _to_integrate(y):
        num = P/E/I*(l*y-y**2/2)
        den = np.sqrt(1-P**2/E**2/I**2*(l*y-y**2/2)**2)
        if np.isnan(den):
            return 100
        else:
            return num/den

    x = []
    for y_i in y:
        x_i = current_L = quad(_to_integrate, 0, y_i)[0]
        x.append(x_i)
    return x 

def consumer_surp(self):
        "Compute consumer surplus"
        # == Compute area under inverse demand function == #
        integrand = lambda x: (self.ad / self.bd) - (1 / self.bd) * x
        area, error = quad(integrand, 0, self.quantity())
        return area - self.price() * self.quantity() 

def overlap(ww, eps):
  """ Overlap matrix between a set of Lorentzians using numerical integration """
  from scipy.integrate import quad 
  n = len(ww)
  mat = zeros((n,n), dtype=np.complex128)
  for i,w1 in enumerate(ww):
    for j,w2 in enumerate(ww):
      re = quad(llc_real, -np.inf, np.inf, args=(w1,w2,eps,eps))
      im = quad(llc_imag, -np.inf, np.inf, args=(w1,w2,eps,eps))
      mat[i,j] = re[0]+1j*im[0]
  return mat 

def moment_of_inertia(self):
        """
        Returns the moment of inertia of the layer [kg m^2]
        """
        moment = 0.0
        radii = self.radii
        density = self.evaluate(['density'])
        rhofunc = UnivariateSpline(radii, density)
        moment = np.abs(quad(lambda r: 8.0 / 3.0 * np.pi * rhofunc(r)
                             * r * r * r * r, radii[0], radii[-1])[0])
        return moment 

def producer_surp(self):
        "Compute producer surplus"
        #  == Compute area above inverse supply curve, excluding tax == #
        integrand = lambda x: -(self.az / self.bz) + (1 / self.bz) * x
        area, error = quad(integrand, 0, self.quantity())  
        return (self.price() - self.tax) * self.quantity() - area 

def arclength(self, chord = None):
        def integrand(x1):
            dr = self.r(x1, 'x1')
            # If function is not well defined everywhere, resulting in a nan
            # penalize it
            if np.isnan(np.sqrt(np.inner(dr, dr))):
                return(100)
            else:
                return np.sqrt(np.inner(dr, dr)[0,0])
        if chord is None:
            chord = self.chord
        return integrate.quad(integrand, 0, chord, limit=500) 

def _x(self, s):
        def _to_minimize(l):
            def _to_integrate(x):
                den = np.sqrt(1-self.G(x)**2)
                if np.isnan(den):
                    return 100
                else:
                    return 1/den
            l = l[0]
            current_L = quad(_to_integrate, 0, l)[0]
            return abs(s-current_L)
        return minimize(_to_minimize, s, method = 'Nelder-Mead',).x[0] 

def find_deflection(self):
        def _to_integrate(x):
            G = self.G(x)
            den = np.sqrt(1-G**2)
            if np.isnan(den):
                return 100
            else:
                return G/den

        self.y = []
        for x_i in self.x:
            y_i = current_L = quad(_to_integrate, 0, x_i)[0]
            self.y.append(y_i) 

def __init__(self, smaknee=30, esigma=0.175/np.sqrt(np.pi/2.), 
            prange=[0.083, 0.882], Rprange=[1, 22.6], **specs):
        
        specs['prange'] = prange
        specs['Rprange'] = Rprange
        PlanetPopulation.__init__(self, **specs)
        
        # calculate norm for sma distribution with decay point (knee)
        self.smaknee = float(smaknee)
        ar = self.arange.to('AU').value
        # sma distribution without normalization
        tmp_dist_sma = lambda x,s0=self.smaknee: x**(-0.62)*np.exp(-(x/s0)**2)
        self.smanorm = integrate.quad(tmp_dist_sma, ar[0], ar[1])[0]
        
        # calculate norm for eccentricity Rayleigh distribution 
        self.esigma = float(esigma)
        er = self.erange
        self.enorm = np.exp(-er[0]**2/(2.*self.esigma**2)) \
                - np.exp(-er[1]**2/(2.*self.esigma**2))
        
        # define Kepler radius distribution
        Rs = np.array([1,1.4,2.0,2.8,4.0,5.7,8.0,11.3,16,22.6]) #Earth Radii
        Rvals85 = np.array([0.1555,0.1671,0.1739,0.0609,0.0187,0.0071,0.0102,0.0049,0.0014])
        #sma of 85 days
        a85 = ((85.*u.day/2./np.pi)**2*u.solMass*const.G)**(1./3.)
        # sma of 0.8 days (lower limit of Fressin et al 2012)
        a08 = ((0.8*u.day/2./np.pi)**2*u.solMass*const.G)**(1./3.) 
        fac1 = integrate.quad(tmp_dist_sma, a08.to('AU').value, a85.to('AU').value)[0]
        Rvals = integrate.quad(tmp_dist_sma, ar[0], ar[1])[0]*(Rvals85/fac1)
        Rvals[5:] *= 2.5 #account for longer orbital baseline data
        self.Rs = Rs
        self.Rvals = Rvals
        self.eta = np.sum(Rvals)
        
        # populate outspec with attributes specific to KeplerLike1
        self._outspec['smaknee'] = self.smaknee
        self._outspec['esigma'] = self.esigma
        self._outspec['eta'] = self.eta
        
        self.dist_albedo_built = None 

def value(self, ts, scenario_index):
        a = 0
        if self._lower_parameter is not None:
            a = self._lower_parameter.get_value(scenario_index)
        b = self._value_to_interpolate(ts, scenario_index)

        cost, err = quad(self.interp, a, b)
        return cost 

def calculate_arc_length(psi_initial, psi_final, A_j, deltaz, c_j):
    """Calculate arc length from psi_initial to psi_final for
       shape coefficient A_j, trailing edge thickness deltaz, and
       chord c_j. Output is the dimensional length"""
    def integrand(psi_baseline, A_j, deltaz, c_j):
        return c_j*np.sqrt(1 + dxi_u(psi_baseline, A_j, deltaz/c_j)**2)

    L, err = quad(integrand, psi_initial, psi_final, args=(A_j, deltaz, c_j))
    return L 

def integral_inf(fn):
  integral, _ = integrate.quad(fn, -np.inf, np.inf)
  return integral 

def find_deflection(y, l, P):
    def _to_integrate(y):
        num = P/E/I*(l*y-y**2/2)
        den = np.sqrt(1-P**2/E**2/I**2*(l*y-y**2/2)**2)
        if np.isnan(den):
            return 100
        else:
            return num/den
    
    x = []
    for y_i in y:
        x_i = current_L = quad(_to_integrate, 0, y_i)[0]
        x.append(x_i)
    return x