Python numpy.polynomial.legendre.leggauss() Examples

def leggauss_ab(n=96, a=-1.0, b=1.0):
  assert(n>0)
  x,w = np.polynomial.legendre.leggauss(n)
  x = (b-a) * 0.5 * x+(b+a) * 0.5
  w = w * (b-a) * 0.5
  return x,w 

def __init__(self, lm=None, jmx=7, ngl=96, lbdmx=14):
    """
    Expansion of the products of two atomic orbitals placed at given locations and around a center between these locations 
    [1] Talman JD. Multipole Expansions for Numerical Orbital Products, Int. J. Quant. Chem. 107, 1578--1584 (2007)
      ngl : order of Gauss-Legendre quadrature
      log_mesh : instance of log_mesh_c defining the logarithmic mesh (rr and pp arrays)
      jmx : maximal angular momentum quantum number of each atomic orbital in a product
      lbdmx : maximal angular momentum quantum number used for the expansion of the product phia*phib
    """
    from numpy.polynomial.legendre import leggauss
    from pyscf.nao.m_log_interp import log_interp_c
    from pyscf.nao.m_csphar_talman_libnao import csphar_talman_libnao as csphar_jt
    assert ngl>2 
    assert jmx>-1
    assert hasattr(lm, 'rr') 
    assert hasattr(lm, 'pp')
    
    self.ngl = ngl
    self.lbdmx = lbdmx
    self.xx, self.ww = leggauss(ngl)
    log_mesh.__init__(self, rr=lm.rr, pp=lm.pp)

    self.plval=np.zeros([2*(self.lbdmx+jmx+1), ngl])
    self.plval[0,:] = 1.0
    self.plval[1,:] = self.xx
    for kappa in range(1,2*(self.lbdmx+jmx)+1):
      self.plval[kappa+1, :]= ((2*kappa+1)*self.xx*self.plval[kappa, :]-kappa*self.plval[kappa-1, :])/(kappa+1)

    self.log_interp = log_interp_c(self.rr)
    self.ylm_cr = csphar_jt([0.0,0.0,1.0], self.lbdmx+2*jmx)

    return 

def test_100(self):
        x, w = leg.leggauss(100)

        # test orthogonality. Note that the results need to be normalized,
        # otherwise the huge values that can arise from fast growing
        # functions like Laguerre can be very confusing.
        v = leg.legvander(x, 99)
        vv = np.dot(v.T * w, v)
        vd = 1/np.sqrt(vv.diagonal())
        vv = vd[:, None] * vv * vd
        assert_almost_equal(vv, np.eye(100))

        # check that the integral of 1 is correct
        tgt = 2.0
        assert_almost_equal(w.sum(), tgt) 

def test_100(self):
        x, w = leg.leggauss(100)

        # test orthogonality. Note that the results need to be normalized,
        # otherwise the huge values that can arise from fast growing
        # functions like Laguerre can be very confusing.
        v = leg.legvander(x, 99)
        vv = np.dot(v.T * w, v)
        vd = 1/np.sqrt(vv.diagonal())
        vv = vd[:, None] * vv * vd
        assert_almost_equal(vv, np.eye(100))

        # check that the integral of 1 is correct
        tgt = 2.0
        assert_almost_equal(w.sum(), tgt) 

def test_100(self):
        x, w = leg.leggauss(100)

        # test orthogonality. Note that the results need to be normalized,
        # otherwise the huge values that can arise from fast growing
        # functions like Laguerre can be very confusing.
        v = leg.legvander(x, 99)
        vv = np.dot(v.T * w, v)
        vd = 1/np.sqrt(vv.diagonal())
        vv = vd[:, None] * vv * vd
        assert_almost_equal(vv, np.eye(100))

        # check that the integral of 1 is correct
        tgt = 2.0
        assert_almost_equal(w.sum(), tgt) 

def glq_nodes_weights(glq_degrees):
    """
    Calculate GLQ unscaled nodes, weights and total number of nodes

    Parameters
    ----------
    glq_degrees : list
        List of GLQ degrees for each direction: ``longitude``, ``latitude``,
        ``radius``.

    Returns
    -------
    n_nodes : int
        Total number of nodes computed as the product of the GLQ degrees.
    glq_nodes : list
        Unscaled GLQ nodes for each direction: ``longitude``, ``latitude``,
        ``radius``.
    glq_weights : list
        GLQ weights for each node on each direction: ``longitude``,
        ``latitude``, ``radius``.
    """
    # Unpack GLQ degrees
    lon_degree, lat_degree, rad_degree = glq_degrees[:]
    # Get number of point masses
    n_nodes = np.prod(glq_degrees)
    # Get nodes coordinates and weights
    lon_node, lon_weights = leggauss(lon_degree)
    lat_node, lat_weights = leggauss(lat_degree)
    rad_node, rad_weights = leggauss(rad_degree)
    # Reorder nodes and weights
    glq_nodes = (lon_node, lat_node, rad_node)
    glq_weights = (lon_weights, lat_weights, rad_weights)
    return n_nodes, glq_nodes, glq_weights 

def test_100(self):
        x, w = leg.leggauss(100)

        # test orthogonality. Note that the results need to be normalized,
        # otherwise the huge values that can arise from fast growing
        # functions like Laguerre can be very confusing.
        v = leg.legvander(x, 99)
        vv = np.dot(v.T * w, v)
        vd = 1/np.sqrt(vv.diagonal())
        vv = vd[:, None] * vv * vd
        assert_almost_equal(vv, np.eye(100))

        # check that the integral of 1 is correct
        tgt = 2.0
        assert_almost_equal(w.sum(), tgt) 

def test_100(self):
        x, w = leg.leggauss(100)

        # test orthogonality. Note that the results need to be normalized,
        # otherwise the huge values that can arise from fast growing
        # functions like Laguerre can be very confusing.
        v = leg.legvander(x, 99)
        vv = np.dot(v.T * w, v)
        vd = 1/np.sqrt(vv.diagonal())
        vv = vd[:, None] * vv * vd
        assert_almost_equal(vv, np.eye(100))

        # check that the integral of 1 is correct
        tgt = 2.0
        assert_almost_equal(w.sum(), tgt) 

def test_100(self):
        x, w = leg.leggauss(100)

        # test orthogonality. Note that the results need to be normalized,
        # otherwise the huge values that can arise from fast growing
        # functions like Laguerre can be very confusing.
        v = leg.legvander(x, 99)
        vv = np.dot(v.T * w, v)
        vd = 1/np.sqrt(vv.diagonal())
        vv = vd[:, None] * vv * vd
        assert_almost_equal(vv, np.eye(100))

        # check that the integral of 1 is correct
        tgt = 2.0
        assert_almost_equal(w.sum(), tgt) 

def test_100(self):
        x, w = leg.leggauss(100)

        # test orthogonality. Note that the results need to be normalized,
        # otherwise the huge values that can arise from fast growing
        # functions like Laguerre can be very confusing.
        v = leg.legvander(x, 99)
        vv = np.dot(v.T * w, v)
        vd = 1/np.sqrt(vv.diagonal())
        vv = vd[:, None] * vv * vd
        assert_almost_equal(vv, np.eye(100))

        # check that the integral of 1 is correct
        tgt = 2.0
        assert_almost_equal(w.sum(), tgt) 

def test_100(self):
        x, w = leg.leggauss(100)

        # test orthogonality. Note that the results need to be normalized,
        # otherwise the huge values that can arise from fast growing
        # functions like Laguerre can be very confusing.
        v = leg.legvander(x, 99)
        vv = np.dot(v.T * w, v)
        vd = 1/np.sqrt(vv.diagonal())
        vv = vd[:, None] * vv * vd
        assert_almost_equal(vv, np.eye(100))

        # check that the integral of 1 is correct
        tgt = 2.0
        assert_almost_equal(w.sum(), tgt) 

def test_100(self):
        x, w = leg.leggauss(100)

        # test orthogonality. Note that the results need to be normalized,
        # otherwise the huge values that can arise from fast growing
        # functions like Laguerre can be very confusing.
        v = leg.legvander(x, 99)
        vv = np.dot(v.T * w, v)
        vd = 1/np.sqrt(vv.diagonal())
        vv = vd[:, None] * vv * vd
        assert_almost_equal(vv, np.eye(100))

        # check that the integral of 1 is correct
        tgt = 2.0
        assert_almost_equal(w.sum(), tgt) 

def linspace(b, grid_type='Driscoll-Healy'):
    if grid_type == 'Driscoll-Healy':
        beta = np.arange(2 * b) * np.pi / (2. * b)
        alpha = np.arange(2 * b) * np.pi / b
    elif grid_type == 'SOFT':
        beta = np.pi * (2 * np.arange(2 * b) + 1) / (4. * b)
        alpha = np.arange(2 * b) * np.pi / b
    elif grid_type == 'Clenshaw-Curtis':
        # beta = np.arange(2 * b + 1) * np.pi / (2 * b)
        # alpha = np.arange(2 * b + 2) * np.pi / (b + 1)
        # Must use np.linspace to prevent numerical errors that cause beta > pi
        beta = np.linspace(0, np.pi, 2 * b + 1)
        alpha = np.linspace(0, 2 * np.pi, 2 * b + 2, endpoint=False)
    elif grid_type == 'Gauss-Legendre':
        x, _ = leggauss(b + 1)  # TODO: leggauss docs state that this may not be only stable for orders > 100
        beta = np.arccos(x)
        alpha = np.arange(2 * b + 2) * np.pi / (b + 1)
    elif grid_type == 'HEALPix':
        #TODO: implement this here so that we don't need the dependency on healpy / healpix_compat
        from healpix_compat import healpy_sphere_meshgrid
        return healpy_sphere_meshgrid(b)
    elif grid_type == 'equidistribution':
        raise NotImplementedError('Not implemented yet; see Fast evaluation of quadrature formulae on the sphere.')
    else:
        raise ValueError('Unknown grid_type:' + grid_type)
    return beta, alpha 

def test_100(self):
        x, w = leg.leggauss(100)

        # test orthogonality. Note that the results need to be normalized,
        # otherwise the huge values that can arise from fast growing
        # functions like Laguerre can be very confusing.
        v = leg.legvander(x, 99)
        vv = np.dot(v.T * w, v)
        vd = 1/np.sqrt(vv.diagonal())
        vv = vd[:, None] * vv * vd
        assert_almost_equal(vv, np.eye(100))

        # check that the integral of 1 is correct
        tgt = 2.0
        assert_almost_equal(w.sum(), tgt) 

def test_100(self):
        x, w = leg.leggauss(100)

        # test orthogonality. Note that the results need to be normalized,
        # otherwise the huge values that can arise from fast growing
        # functions like Laguerre can be very confusing.
        v = leg.legvander(x, 99)
        vv = np.dot(v.T * w, v)
        vd = 1/np.sqrt(vv.diagonal())
        vv = vd[:, None] * vv * vd
        assert_almost_equal(vv, np.eye(100))

        # check that the integral of 1 is correct
        tgt = 2.0
        assert_almost_equal(w.sum(), tgt) 

def test_100(self):
        x, w = leg.leggauss(100)

        # test orthogonality. Note that the results need to be normalized,
        # otherwise the huge values that can arise from fast growing
        # functions like Laguerre can be very confusing.
        v = leg.legvander(x, 99)
        vv = np.dot(v.T * w, v)
        vd = 1/np.sqrt(vv.diagonal())
        vv = vd[:, None] * vv * vd
        assert_almost_equal(vv, np.eye(100))

        # check that the integral of 1 is correct
        tgt = 2.0
        assert_almost_equal(w.sum(), tgt) 

def test_100(self):
        x, w = leg.leggauss(100)

        # test orthogonality. Note that the results need to be normalized,
        # otherwise the huge values that can arise from fast growing
        # functions like Laguerre can be very confusing.
        v = leg.legvander(x, 99)
        vv = np.dot(v.T * w, v)
        vd = 1/np.sqrt(vv.diagonal())
        vv = vd[:, None] * vv * vd
        assert_almost_equal(vv, np.eye(100))

        # check that the integral of 1 is correct
        tgt = 2.0
        assert_almost_equal(w.sum(), tgt) 

def test_100(self):
        x, w = leg.leggauss(100)

        # test orthogonality. Note that the results need to be normalized,
        # otherwise the huge values that can arise from fast growing
        # functions like Laguerre can be very confusing.
        v = leg.legvander(x, 99)
        vv = np.dot(v.T * w, v)
        vd = 1/np.sqrt(vv.diagonal())
        vv = vd[:, None] * vv * vd
        assert_almost_equal(vv, np.eye(100))

        # check that the integral of 1 is correct
        tgt = 2.0
        assert_almost_equal(w.sum(), tgt) 

def test_100(self):
        x, w = leg.leggauss(100)

        # test orthogonality. Note that the results need to be normalized,
        # otherwise the huge values that can arise from fast growing
        # functions like Laguerre can be very confusing.
        v = leg.legvander(x, 99)
        vv = np.dot(v.T * w, v)
        vd = 1/np.sqrt(vv.diagonal())
        vv = vd[:, None] * vv * vd
        assert_almost_equal(vv, np.eye(100))

        # check that the integral of 1 is correct
        tgt = 2.0
        assert_almost_equal(w.sum(), tgt) 

def quadrature_weights(b, grid_type='Gauss-Legendre'):
    """
    Compute quadrature weights for a given grid-type.
    The function S2.meshgrid generates the points that correspond to the weights generated by this function.

    if convention == 'Gauss-Legendre':
    The quadrature formula is exact for polynomials up to degree M less than or equal to 2b + 1,
    so that we can compute exact Fourier coefficients for f a polynomial of degree at most b.

    if convention == 'Clenshaw-Curtis':
    The quadrature formula is exact for polynomials up to degree M less than or equal to 2b,
    so that we can compute exact Fourier coefficients for f a polynomial of degree at most b.

    :param b: the grid resolution. See S2.meshgrid
    :param grid_type:
    :return:
    """
    if grid_type == 'Clenshaw-Curtis':
        # There is a faster fft based method to compute these weights
        # see "Fast evaluation of quadrature formulae on the sphere"
        # W = np.empty((2 * b + 2, 2 * b + 1))
        # for j in range(2 * b + 1):
        #    eps_j_2b = 0.5 if j == 0 or j == 2 * b else 1.
        #    for k in range(2 * b + 2):  # Doesn't seem to depend on k..
        #        W[k, j] = (4 * np.pi * eps_j_2b) / (b * (2 * b + 2))
        #        sum = 0.
        #        for l in range(b + 1):
        #            eps_l_b = 0.5 if l == 0 or l == b else 1.
        #            sum += eps_l_b / (1 - 4 * l ** 2) * np.cos(j * l * np.pi / b)
        #        W[k, j] *= sum
        w = _clenshaw_curtis_weights(n=2 * b)
        W = np.empty((2 * b + 1, 2 * b + 2))
        W[:] = w[:, None]
    elif grid_type == 'Gauss-Legendre':
        # We found this formula in:
        # "A Fast Algorithm for Spherical Grid Rotations and its Application to Singular Quadrature"
        # eq. 10
        _, w = leggauss(b + 1)
        W = w[:, None] * (2 * np.pi / (2 * b + 2) * np.ones(2 * b + 2)[None, :])
    elif grid_type == 'SOFT':
        print("WARNING: SOFT quadrature weights don't work yet")
        k = np.arange(0, b)
        w = np.array([(2. / b) * np.sin(np.pi * (2. * j + 1.) / (4. * b)) *
                      (np.sum((1. / (2 * k + 1))
                              * np.sin((2 * j + 1) * (2 * k + 1)
                                       * np.pi / (4. * b))))
                      for j in range(2 * b)])
        W = w[:, None] * np.ones(2 * b)[None, :]
    else:
        raise ValueError('Unknown grid_type:' + str(grid_type))

    return W 

def GaussQuadrature(N,a=-1,b=1):
    if a==-1 and b==1:
        return leggauss(N)
    # The following is for historical purposes and when the range is different from [-1,1]
    N0=N-1
    N1 = N0+1
    N2 = N0+2
    xu = np.linspace(-1.,1.,N1)
    
    # Legendre-Gauss-Vandermonde Matrix
    L = 1.0*np.zeros((N1,N2))
    # Derivative of Legendre-Gauss-Vandermonde Matrix
    Lp = 1.0*np.zeros(N1)

    dum = np.linspace(0,N0,N1)
    y=np.cos((2*dum+1)*np.pi/(2*N0+2))+(0.27/N1)*np.sin(np.pi*xu*N0/N2)
    # PI = np.pi
    # y=ne.evaluate("cos((2*dum+1)*PI/(2*N0+2))+(0.27/N1)*sin(PI*xu*N0/N2)")
    deps = np.finfo(np.float64).eps

    # Initial Guess
    y0 = 2.0*np.ones(N1)

    while np.max(np.abs(y-y0)) > deps:
        L[:,0] = np.ones(N1)
        L[:,1] = y
        Lp  = np.zeros(N1)

        for k in range(1,N1):
            L[:,k+1] = ((2*k+1)*L[:,k]*y - k*L[:,k-1])/(k+1)

        Lp = N2*(L[:,N0]-L[:,N1]*y)/(1-y**2)

        y0 = y
        y=y0-L[:,N1]/Lp

    z = (a*(1-y)+b*(1+y))/2.0
    w = (b-a)/((1-y**2)*Lp**2)*pow((np.float64(N2)/N1),2)


    z = np.fliplr(z.reshape(1,z.shape[0])).reshape(z.shape[0]) 
    w = np.fliplr(w.reshape(1,w.shape[0])).reshape(w.shape[0]) 

    return (z,w)