Python scipy.integrate.dblquad() Examples

def discretize_integrate_2D(model, x_range, y_range):
    """
    Discretize model by integrating the model over the pixel.
    """
    from scipy.integrate import dblquad
    # Set up grid
    x = np.arange(x_range[0] - 0.5, x_range[1] + 0.5)
    y = np.arange(y_range[0] - 0.5, y_range[1] + 0.5)
    values = np.empty((y.size - 1, x.size - 1))

    # Integrate over all pixels
    for i in range(x.size - 1):
        for j in range(y.size - 1):
            values[j, i] = dblquad(lambda y, x: model(x, y), x[i], x[i + 1],
                                   lambda x: y[j], lambda x: y[j + 1])[0]
    return values 

def precise_chern(h,dk=0.01, mode="Wilson",delta=0.0001,operator=None):
  """ Calculates the chern number of a 2d system """
  from scipy import integrate
  err = {"epsabs" : 1.0, "epsrel": 1.0,"limit" : 10}
#  err = [err,err]
  def f(x,y): # function to integrate
    if mode=="Wilson":
      return berry_curvature(h,np.array([x,y]),dk=dk)
    if mode=="Green":
       f2 = h.get_gk_gen(delta=delta) # get generator
       return berry_green(f2,k=[x,y,0.],operator=operator) 
  c = integrate.dblquad(f,0.,1.,lambda x : 0., lambda x: 1.,epsabs=0.01,
                          epsrel=0.01)
  chern = c[0]/(2.*np.pi)
  open("CHERN.OUT","w").write(str(chern)+"\n")
  return chern 

def bivariate_approximate_ground_truth_integral(func, integral_bounds: List[Tuple[float, float]]):
    """
    Estimate the 2D ground truth integral

    :param func: bivariate function
    :param integral_bounds: bounds of integral
    :returns: integral estimate, output of scipy.integrate.dblquad
    """
    def func_dblquad(x, y):
        z = np.array([x, y])
        z = np.reshape(z, [1, 2])
        return func(z)

    lower_bound_x = integral_bounds[0][0]
    upper_bound_x = integral_bounds[0][1]

    lower_bound_y = integral_bounds[1][0]
    upper_bound_y = integral_bounds[1][1]

    return dblquad(func=func_dblquad, a=lower_bound_x, b=upper_bound_x, gfun=lambda x: lower_bound_y,
                   hfun=lambda x: upper_bound_y) 

def test_double_integral3(self):
        def func(x0, x1):
            return x0 + x1 + 1 + 2
        assert_quad(dblquad(func, 1, 2, 1, 2),6.) 

def precise_spin_chern(h,delta=0.00001,tol=0.1):
  """ Calculates the chern number of a 2d system """
  from scipy import integrate
  err = {"epsabs" : 0.01, "epsrel": 0.01,"limit" : 20}
  sz = operators.get_sz(h) # get sz operator
  def f(x,y): # function to integrate
    return operator_berry(h,np.array([x,y]),delta=delta,operator=sz)
  c = integrate.dblquad(f,0.,1.,lambda x : 0., lambda x: 1.,epsabs=tol,
                          epsrel=tol)
  return c[0]/(2.*np.pi) 

def total_energy(h,nk=10,nbands=None,use_kpm=False,random=False,
        kp=None,mode="mesh",tol=1e-1):
  """Return the total energy"""
  h.turn_dense()
  if h.is_sparse and not use_kpm: 
    print("Sparse Hamiltonian but no bands given, taking 20")
    nbands=20
  f = h.get_hk_gen() # get generator
  etot = 0.0 # initialize
  iv = 0
  def enek(k):
    """Compute energy in this kpoint"""
    hk = f(k)  # kdependent hamiltonian
    if use_kpm: # Kernel polynomial method
      return kpm.total_energy(hk,scale=10.,ntries=20,npol=100) # using KPM
    else: # conventional diagonalization
      if nbands is None: vv = lg.eigvalsh(hk) # diagonalize k hamiltonian
      else: 
          vv,aa = slg.eigsh(hk,k=4*nbands,which="LM",sigma=0.0) 
          vv = -np.sort(-(vv[vv<0.0])) # negative eigenvalues
          vv = vv[0:nbands] # get the negative eigenvlaues closest to EF
      return np.sum(vv[vv<0.0]) # sum energies below fermi energy
  # compute energy using different modes
  if mode=="mesh":
    from .klist import kmesh
    kp = kmesh(h.dimensionality,nk=nk)
    etot = np.mean(parallel.pcall(enek,kp)) # compute total energy
  elif mode=="random":
    kp = [np.random.random(3) for i in range(nk)] # random points
    etot = np.mean(parallel.pcall(enek,kp)) # compute total eenrgy
  elif mode=="integrate":
    from scipy import integrate
    if h.dimensionality==1: # one dimensional
        etot = integrate.quad(enek,-1.,1.,epsabs=tol,epsrel=tol)[0]
    elif h.dimensionality==2: # two dimensional
        etot = integrate.dblquad(lambda x,y: enek([x,y]),-1.,1.,-1.,1.,
                epsabs=tol,epsrel=tol)[0]
    else: raise
  else: raise
  return etot 

def bivariateTCdf(yy,xx,rho,nu):    
    t_ans, err = nInt.dblquad(bivariateTDensity, -10, xx,
                   lambda x: -10,
                   lambda x: yy,args=(rho,nu))
    return t_ans 

def quad2d(func=lambda x: 1, lower=(-10,-10), upper=(10,10)):
    def fun(x, y):
        x = np.column_stack((x,y))
        return func(x)
    from scipy.integrate import dblquad
    return dblquad(fun, lower[0], upper[0], lambda y: lower[1],
                   lambda y: upper[1]) 

def kl(self, other):
        '''Kullback-Leibler divergence between this and another distribution

        int f(x) (log f(x) - log g(x)) dx

        where f is the pdf of self, and g is the pdf of other

        uses double integration with scipy.integrate.dblquad

        limits currently hardcoded

        '''
        fun = lambda x : self.logpdf(x) - other.logpdf(x)
        return self.expect(fun) 

def expect(self, func=lambda x: 1, lower=(-10,-10), upper=(10,10)):
        def fun(x, y):
            x = np.column_stack((x,y))
            return func(x) * self.pdf(x)
        from scipy.integrate import dblquad
        return dblquad(fun, lower[0], upper[0], lambda y: lower[1],
                       lambda y: upper[1]) 

def test_matching_dblquad(self):
        def func2d(x0, x1):
            return x0**2 + x1**3 - x0 * x1 + 1

        res, reserr = dblquad(func2d, -2, 2, lambda x: -3, lambda x: 3)
        res2, reserr2 = nquad(func2d, [[-3, 3], (-2, 2)])
        assert_almost_equal(res, res2)
        assert_almost_equal(reserr, reserr2) 

def mutual_information(self):
        def mi_integrand(x, y):
            return np.exp(self.logpdf_xy(x,y)) * \
                (self.logpdf_xy(x,y) - self.logpdf_x(x) - self.logpdf_y(y))
        return dblquad(
            lambda y, x: mi_integrand(x,y), self.D[0], self.D[1],
            self._lower_y, self._upper_y) 

def test_double_integral2(self):
        def func(x0, x1, t0, t1):
            return x0 + x1 + t0 + t1
        g = lambda x: x
        h = lambda x: 2 * x
        args = 1, 2
        assert_quad(dblquad(func, 1, 2, g, h, args=args),35./6 + 9*.5) 

def test_double_integral(self):
        # 8) Double Integral test
        def simpfunc(y, x):       # Note order of arguments.
            return x+y

        a, b = 1.0, 2.0
        assert_quad(dblquad(simpfunc, a, b, lambda x: x, lambda x: 2*x),
                    5/6.0 * (b**3.0-a**3.0)) 

def test_matching_dblquad(self):
        def func2d(x0, x1):
            return x0**2 + x1**3 - x0 * x1 + 1

        res, reserr = dblquad(func2d, -2, 2, lambda x:-3, lambda x:3)
        res2, reserr2 = nquad(func2d, [[-3, 3], (-2, 2)])
        assert_almost_equal(res, res2)
        assert_almost_equal(reserr, reserr2) 

def test_double_integral(self):
        # 8) Double Integral test
        def simpfunc(y,x):       # Note order of arguments.
            return x+y

        a, b = 1.0, 2.0
        assert_quad(dblquad(simpfunc,a,b,lambda x: x, lambda x: 2*x),
                    5/6.0 * (b**3.0-a**3.0)) 

def quad2d(func=lambda x: 1, lower=(-10,-10), upper=(10,10)):
    def fun(x, y):
        x = np.column_stack((x,y))
        return func(x)
    from scipy.integrate import dblquad
    return dblquad(fun, lower[0], upper[0], lambda y: lower[1],
                   lambda y: upper[1]) 

def kl(self, other):
        '''Kullback-Leibler divergence between this and another distribution

        int f(x) (log f(x) - log g(x)) dx

        where f is the pdf of self, and g is the pdf of other

        uses double integration with scipy.integrate.dblquad

        limits currently hardcoded

        '''
        fun = lambda x : self.logpdf(x) - other.logpdf(x)
        return self.expect(fun) 

def expect(self, func=lambda x: 1, lower=(-10,-10), upper=(10,10)):
        def fun(x, y):
            x = np.column_stack((x,y))
            return func(x) * self.pdf(x)
        from scipy.integrate import dblquad
        return dblquad(fun, lower[0], upper[0], lambda y: lower[1],
                       lambda y: upper[1]) 

def _sanity_test(self):
        # Marginal of x integrates to one.
        print quad(lambda x: np.exp(self.logpdf_x(x)), self.D[0], self.D[1])

        # Marginal of y integrates to one.
        print quad(lambda y: np.exp(self.logpdf_y(y)), -1 ,1)

        # Joint of x,y integrates to one; quadrature will fail for small noise.
        print dblquad(
            lambda y,x: np.exp(self.logpdf_xy(x,y)), self.D[0], self.D[1],
            lambda x: -1, lambda x: 1) 

def spherical_integral(C,rho):
    """
    Calculate the integral of a function over a unit sphere.
    """
    # phi - azimuthal angle (angle in xy-plane)
    # theta - polar angle (angle between z and xy-plane)
    #       (  y  , x )
    def func(theta,phi,C,rho):  # Test function. Can I get 4*pi^2????
        x = sp.cos(phi)*sp.sin(theta)
        y = sp.sin(phi)*sp.sin(theta)
        z = sp.cos(theta)
        #dir = sp.array((x,y,z))
        #dc = dir_cosines(dir)
        dc = sp.array((x,y,z))  # Turns out these are direction cosines!
        Gamma = make_gamma(dc,C)
        rho_c_square = linalg.eigvals(Gamma).real  # GPa
        rho_c_square = rho_c_square*1e9  # Pa
        sound_vel = sp.sqrt(rho_c_square/rho)  # m/s
        integrand = 1/(sound_vel[0]**3) + 1/(sound_vel[1]**3) + 1/(sound_vel[2]**3)
        return integrand*sp.sin(theta)
    #        (  y  , x      )
    #def sfunc(theta,phi,args=()):
    #    return func(theta,phi,args)*sp.sin(theta)

    integral,error = dblquad(func,0,2*sp.pi,lambda g: 0,lambda h:
            sp.pi,args=(C,rho))
    return integral


#direction = sp.array((1.0,1.0,1.0))
#dc = dir_cosines(direction)
#C = read_file.read_file(sys.argv[1])
#C.pop(0)
#C = sp.array(C,float)
#Gamma = make_gamma(dc,C)
#density = 7500 #kg/m**3
#density = float(read_file.read_file(sys.argv[2])[0][0])
#rho_c_square = linalg.eigvals(Gamma) #GPa
#rho_c_square = rho_c_square*1e9 #Pa
#sound_vel = sp.sqrt(rho_c_square/density).real
#avg_vel = sp.average(sound_vel)
#print Gamma
#print direction
#print C
#print rho_c_square
#print rho_c_square.real
#print sound_vel," in m/s"
#print avg_vel
#print spherical_integral(C,density)