Python special.lpmv方法代碼示例

# 需要導入模塊: from scipy import special [as 別名]
# 或者: from scipy.special import lpmv [as 別名]
def test_erf_complex():
    # need to increase mpmath precision for this test
    old_dps, old_prec = mpmath.mp.dps, mpmath.mp.prec
    try:
        mpmath.mp.dps = 70
        x1, y1 = np.meshgrid(np.linspace(-10, 1, 31), np.linspace(-10, 1, 11))
        x2, y2 = np.meshgrid(np.logspace(-80, .8, 31), np.logspace(-80, .8, 11))
        points = np.r_[x1.ravel(),x2.ravel()] + 1j*np.r_[y1.ravel(),y2.ravel()]

        assert_func_equal(sc.erf, lambda x: complex(mpmath.erf(x)), points,
                          vectorized=False, rtol=1e-13)
        assert_func_equal(sc.erfc, lambda x: complex(mpmath.erfc(x)), points,
                          vectorized=False, rtol=1e-13)
    finally:
        mpmath.mp.dps, mpmath.mp.prec = old_dps, old_prec


#------------------------------------------------------------------------------
# lpmv
#------------------------------------------------------------------------------ 

# 需要導入模塊: from scipy import special [as 別名]
# 或者: from scipy.special import lpmv [as 別名]
def test_legenp(self):
        def lpnm(n, m, z):
            if m > n:
                return 0.0
            return sc.lpmn(m, n, z)[0][-1,-1]

        def lpnm_2(n, m, z):
            if m > n:
                return 0.0
            return sc.lpmv(m, n, z)

        def legenp(n, m, z):
            if abs(z) < 1e-15:
                # mpmath has bad performance here
                return np.nan
            return _exception_to_nan(mpmath.legenp)(n, m, z, **HYPERKW)

        assert_mpmath_equal(lpnm,
                            legenp,
                            [IntArg(0, 100), IntArg(0, 100), Arg()])

        assert_mpmath_equal(lpnm_2,
                            legenp,
                            [IntArg(0, 100), IntArg(0, 100), Arg(-1, 1)]) 

# 需要導入模塊: from scipy import special [as 別名]
# 或者: from scipy.special import lpmv [as 別名]
def test_erf_complex():
    # need to increase mpmath precision for this test
    old_dps, old_prec = mpmath.mp.dps, mpmath.mp.prec
    try:
        mpmath.mp.dps = 70
        x1, y1 = np.meshgrid(np.linspace(-10, 1, 31), np.linspace(-10, 1, 11))
        x2, y2 = np.meshgrid(np.logspace(-80, .8, 31), np.logspace(-80, .8, 11))
        points = np.r_[x1.ravel(),x2.ravel()] + 1j*np.r_[y1.ravel(), y2.ravel()]

        assert_func_equal(sc.erf, lambda x: complex(mpmath.erf(x)), points,
                          vectorized=False, rtol=1e-13)
        assert_func_equal(sc.erfc, lambda x: complex(mpmath.erfc(x)), points,
                          vectorized=False, rtol=1e-13)
    finally:
        mpmath.mp.dps, mpmath.mp.prec = old_dps, old_prec


# ------------------------------------------------------------------------------
# lpmv
# ------------------------------------------------------------------------------ 

# 需要導入模塊: from scipy import special [as 別名]
# 或者: from scipy.special import lpmv [as 別名]
def spherical_harmonics(m, n, theta, phi):
    """
    An implementation of spherical harmonics that overcomes conda compilation
    issues. See: https://github.com/nipy/dipy/issues/852
    """
    x = np.cos(phi)
    val = lpmv(m, n, x).astype(complex)
    val *= np.sqrt((2 * n + 1) / 4.0 / np.pi)
    val *= np.exp(0.5 * (gammaln(n - m + 1) - gammaln(n + m + 1)))
    val = val * np.exp(1j * m * theta)
    return val 

# 需要導入模塊: from scipy import special [as 別名]
# 或者: from scipy.special import lpmv [as 別名]
def test_lpmv(self):
        assert_equal(cephes.lpmv(0,0,1),1.0) 

# 需要導入模塊: from scipy import special [as 別名]
# 或者: from scipy.special import lpmv [as 別名]
def test_clpmn_close_to_real_2(self):
        eps = 1e-10
        m = 1
        n = 3
        x = 0.5
        clp_plus = special.clpmn(m, n, x+1j*eps, 2)[0][m, n]
        clp_minus = special.clpmn(m, n, x-1j*eps, 2)[0][m, n]
        assert_array_almost_equal(array([clp_plus, clp_minus]),
                                  array([special.lpmv(m, n, x),
                                         special.lpmv(m, n, x)]),
                                  7) 

# 需要導入模塊: from scipy import special [as 別名]
# 或者: from scipy.special import lpmv [as 別名]
def test_clpmn_close_to_real_3(self):
        eps = 1e-10
        m = 1
        n = 3
        x = 0.5
        clp_plus = special.clpmn(m, n, x+1j*eps, 3)[0][m, n]
        clp_minus = special.clpmn(m, n, x-1j*eps, 3)[0][m, n]
        assert_array_almost_equal(array([clp_plus, clp_minus]),
                                  array([special.lpmv(m, n, x)*np.exp(-0.5j*m*np.pi),
                                         special.lpmv(m, n, x)*np.exp(0.5j*m*np.pi)]),
                                  7) 

# 需要導入模塊: from scipy import special [as 別名]
# 或者: from scipy.special import lpmv [as 別名]
def _naive_csh_unnormalized(l, m, theta, phi):
    """
    Compute unnormalized SH
    """
    return lpmv(m, l, np.cos(theta)) * np.exp(1j * m * phi) 

# 需要導入模塊: from scipy import special [as 別名]
# 或者: from scipy.special import lpmv [as 別名]
def _naive_csh_quantum(l, m, theta, phi):
    """
    Compute orthonormalized spherical harmonics in a naive way.
    """
    return (((-1.) ** m) * lpmv(m, l, np.cos(theta)) * np.exp(1j * m * phi) *
            np.sqrt(((2 * l + 1) * factorial(l - m))
                    /
                    (4 * np.pi * factorial(l + m)))) 

# 需要導入模塊: from scipy import special [as 別名]
# 或者: from scipy.special import lpmv [as 別名]
def _naive_csh_seismology(l, m, theta, phi):
    """
    Compute the spherical harmonics according to the seismology convention, in a naive way.
    This appears to be equal to the sph_harm function in scipy.special.
    """
    return (lpmv(m, l, np.cos(theta)) * np.exp(1j * m * phi) *
            np.sqrt(((2 * l + 1) * factorial(l - m))
                    /
                    (4 * np.pi * factorial(l + m)))) 

# 需要導入模塊: from scipy import special [as 別名]
# 或者: from scipy.special import lpmv [as 別名]
def _naive_csh_ph(l, m, theta, phi):
    """
    CSH as defined by Pinchon-Hoggan. Same as wikipedia's quantum-normalized SH = naive_Y_quantum()
    """
    if l == 0 and m == 0:
        return 1. / np.sqrt(4 * np.pi)
    else:
        phase = ((1j) ** (m + np.abs(m)))
        normalizer = np.sqrt(((2 * l + 1.) * factorial(l - np.abs(m)))
                             /
                             (4 * np.pi * factorial(l + np.abs(m))))
        P = lpmv(np.abs(m), l, np.cos(theta))
        e = np.exp(1j * m * phi)
        return phase * normalizer * P * e 

# 需要導入模塊: from scipy import special [as 別名]
# 或者: from scipy.special import lpmv [as 別名]
def _naive_rsh_ph(l, m, theta, phi):

    if m == 0:
        return np.sqrt((2 * l + 1.) / (4 * np.pi)) * lpmv(m, l, np.cos(theta))
    elif m < 0:
        return np.sqrt(((2 * l + 1.) * factorial(l + m)) /
                       (2 * np.pi * factorial(l - m))) * lpmv(-m, l, np.cos(theta)) * np.sin(-m * phi)
    elif m > 0:
        return np.sqrt(((2 * l + 1.) * factorial(l - m)) /
                       (2 * np.pi * factorial(l + m))) * lpmv(m, l, np.cos(theta)) * np.cos(m * phi) 

# 需要導入模塊: from scipy import special [as 別名]
# 或者: from scipy.special import lpmv [as 別名]
def _rspherical_harm(m, l, theta, cos_phi):
    r""" Calculates the real spherical harmonics using :math:`Y_l^m(\theta, \varphi)` with :math:`\mathbf R\to \{r, \theta, \varphi\}`.

    These real spherical harmonics are via these equations:

    .. math::
        Y^m_l(\theta,\varphi) &= -(-1)^m\sqrt{2\frac{2l+1}{4\pi} \frac{(l-m)!}{(l+m)!}}
           P^{m}_l (\cos(\varphi)) \sin(m \theta) & m < 0\\
        Y^m_l(\theta,\varphi) &= \sqrt{\frac{2l+1}{4\pi}} P^{m}_l (\cos(\varphi)) & m = 0\\
        Y^m_l(\theta,\varphi) &= \sqrt{2\frac{2l+1}{4\pi} \frac{(l-m)!}{(l+m)!}}
           P^{m}_l (\cos(\varphi)) \cos(m \theta) & m > 0

    Parameters
    ----------
    m : int
       order of the spherical harmonics
    l : int
       degree of the spherical harmonics
    theta : array_like
       angle in :math:`x-y` plane (azimuthal)
    cos_phi : array_like
       cos(phi) to angle from :math:`z` axis (polar)
    """
    # Calculate the associated Legendre polynomial
    # Since the real spherical harmonics has slight differences
    # for positive and negative m, we have to implement them individually.
    # Currently this is a re-write of what Inelastica does and a combination of
    # learned lessons from Denchar.
    # As such the choice of these real spherical harmonics is that of Siesta.
    if m == 0:
        return _rspher_harm_fact[l][m] * lpmv(m, l, cos_phi)
    elif m < 0:
        return _rspher_harm_fact[l][m] * (lpmv(m, l, cos_phi) * sin(m*theta))
    return _rspher_harm_fact[l][m] * (lpmv(m, l, cos_phi) * cos(m*theta)) 

# 需要導入模塊: from scipy import special [as 別名]
# 或者: from scipy.special import lpmv [as 別名]
def gscontrol_raw(dtrank=4):
	"""
	This function uses the spatial global signal estimation approach to modify catd (global variable) to 
	removal global signal out of individual echo time series datasets. The spatial global signal is estimated
	from the optimally combined data after detrending with a Legendre polynomial basis of order=0 and degree=dtrank.
	"""

	print "++ Applying amplitude-based T1 equilibration correction"
    
	#Legendre polynomial basis for denoising
	from scipy.special import lpmv
	Lmix = np.array([lpmv(0,vv,np.linspace(-1,1,OCcatd.shape[-1])) for vv in range(dtrank)]).T

	#Compute mean, std, mask local to this function - inefficient, but makes this function a bit more modular
	Gmu = OCcatd.mean(-1)
	Gstd = OCcatd.std(-1)
	Gmask = Gmu!=0

	#Find spatial global signal
	dat = OCcatd[Gmask] - Gmu[Gmask][:,np.newaxis]
	sol = np.linalg.lstsq(Lmix,dat.T) #Legendre basis for detrending
	detr = dat - np.dot(sol[0].T,Lmix.T)[0]
	sphis = (detr).min(1)
	sphis -= sphis.mean()
	niwrite(unmask(sphis,Gmask),aff,'T1gs.nii',head)
    
	#Find time course of the spatial global signal, make basis with the Legendre basis
	glsig = np.linalg.lstsq(np.atleast_2d(sphis).T,dat)[0]
	glsig = (glsig-glsig.mean() ) / glsig.std()
	np.savetxt('glsig.1D',glsig)
	glbase = np.hstack([Lmix,glsig.T])

	#Project global signal out of optimally combined data
	sol = np.linalg.lstsq(np.atleast_2d(glbase),dat.T)
	tsoc_nogs = dat - np.dot(np.atleast_2d(sol[0][dtrank]).T,np.atleast_2d(glbase.T[dtrank])) + Gmu[Gmask][:,np.newaxis]

	global OCcatd, sphis_raw
	sphis_raw = sphis
	niwrite(OCcatd,aff,'tsoc_orig.nii',head)
	OCcatd = unmask(tsoc_nogs,Gmask)
	niwrite(OCcatd,aff,'tsoc_nogs.nii',head)
	
	#Project glbase out of each echo
	for ii in range(Ne):
		dat = catd[:,:,:,ii,:][Gmask]
		sol = np.linalg.lstsq(np.atleast_2d(glbase),dat.T)
		e_nogs = dat - np.dot(np.atleast_2d(sol[0][dtrank]).T,np.atleast_2d(glbase.T[dtrank]))
		catd[:,:,:,ii,:] = unmask(e_nogs,Gmask)
		if 'DEBUG' in sys.argv:
			ipdb.set_trace() 

# 需要導入模塊: from scipy import special [as 別名]
# 或者: from scipy.special import lpmv [as 別名]
def test_lpmv():
    pts = []
    for x in [-0.99, -0.557, 1e-6, 0.132, 1]:
        pts.extend([
            (1, 1, x),
            (1, -1, x),
            (-1, 1, x),
            (-1, -2, x),
            (1, 1.7, x),
            (1, -1.7, x),
            (-1, 1.7, x),
            (-1, -2.7, x),
            (1, 10, x),
            (1, 11, x),
            (3, 8, x),
            (5, 11, x),
            (-3, 8, x),
            (-5, 11, x),
            (3, -8, x),
            (5, -11, x),
            (-3, -8, x),
            (-5, -11, x),
            (3, 8.3, x),
            (5, 11.3, x),
            (-3, 8.3, x),
            (-5, 11.3, x),
            (3, -8.3, x),
            (5, -11.3, x),
            (-3, -8.3, x),
            (-5, -11.3, x),
        ])

    def mplegenp(nu, mu, x):
        if mu == int(mu) and x == 1:
            # mpmath 0.17 gets this wrong
            if mu == 0:
                return 1
            else:
                return 0
        return mpmath.legenp(nu, mu, x)

    dataset = [p + (mplegenp(p[1], p[0], p[2]),) for p in pts]
    dataset = np.array(dataset, dtype=np.float_)

    def evf(mu, nu, x):
        return sc.lpmv(mu.astype(int), nu, x)

    olderr = np.seterr(invalid='ignore')
    try:
        FuncData(evf, dataset, (0,1,2), 3, rtol=1e-10, atol=1e-14).check()
    finally:
        np.seterr(**olderr)


#------------------------------------------------------------------------------
# beta
#------------------------------------------------------------------------------