Python scipy.special.lpmv() Examples

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
#------------------------------------------------------------------------------ 

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)]) 

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
# ------------------------------------------------------------------------------ 

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) 

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)) 

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) 

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 

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)))) 

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)))) 

def _naive_csh_unnormalized(l, m, theta, phi):
    """
    Compute unnormalized SH
    """
    return lpmv(m, l, np.cos(theta)) * np.exp(1j * m * phi) 

def test_lpmv(self):
        lp = special.lpmv(0,2,.5)
        assert_almost_equal(lp,-0.125,7)
        lp = special.lpmv(0,40,.001)
        assert_almost_equal(lp,0.1252678976534484,7)

        # XXX: this is outside the domain of the current implementation,
        #      so ensure it returns a NaN rather than a wrong answer.
        olderr = np.seterr(all='ignore')
        try:
            lp = special.lpmv(-1,-1,.001)
        finally:
            np.seterr(**olderr)
        assert_(lp != 0 or np.isnan(lp)) 

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) 

def test_lpmv(self):
        assert_equal(cephes.lpmv(0,0,1),1.0) 

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 

def test_lpmv(self):
        lp = special.lpmv(0,2,.5)
        assert_almost_equal(lp,-0.125,7)
        lp = special.lpmv(0,40,.001)
        assert_almost_equal(lp,0.1252678976534484,7)

        # XXX: this is outside the domain of the current implementation,
        #      so ensure it returns a NaN rather than a wrong answer.
        olderr = np.seterr(all='ignore')
        try:
            lp = special.lpmv(-1,-1,.001)
        finally:
            np.seterr(**olderr)
        assert_(lp != 0 or np.isnan(lp)) 

def test_lpmv(self):
        assert_equal(cephes.lpmv(0,0,1),1.0) 

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() 

def test_legenp(self):
        def lpnm(n, m, z):
            try:
                v = sc.lpmn(m, n, z)[0][-1,-1]
            except ValueError:
                return np.nan
            if abs(v) > 1e306:
                # harmonize overflow to inf
                v = np.inf * np.sign(v.real)
            return v

        def lpnm_2(n, m, z):
            v = sc.lpmv(m, n, z)
            if abs(v) > 1e306:
                # harmonize overflow to inf
                v = np.inf * np.sign(v.real)
            return v

        def legenp(n, m, z):
            if (z == 1 or z == -1) and int(n) == n:
                # Special case (mpmath may give inf, we take the limit by
                # continuity)
                if m == 0:
                    if n < 0:
                        n = -n - 1
                    return mpmath.power(mpmath.sign(z), n)
                else:
                    return 0

            if abs(z) < 1e-15:
                # mpmath has bad performance here
                return np.nan

            typ = 2 if abs(z) < 1 else 3
            v = exception_to_nan(mpmath.legenp)(n, m, z, type=typ)

            if abs(v) > 1e306:
                # harmonize overflow to inf
                v = mpmath.inf * mpmath.sign(v.real)

            return v

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

        assert_mpmath_equal(lpnm_2,
                            legenp,
                            [IntArg(-100, 100), Arg(-100, 100), Arg(-1, 1)],
                            atol=1e-10) 

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
#------------------------------------------------------------------------------