Python scipy.special 模块，legendre() 实例源码
我们从Python开源项目中，提取了以下14个代码示例，用于说明如何使用scipy.special.legendre()。
def test_chebyshev_modified(tol=1.0e-14):
alpha = 2.0
# Get the moments corresponding to the Legendre polynomials and the weight
# function omega(x) = |x^alpha|:
#
# / 2/3 if k == 0,
# int_{-1}^{+1} |x^alpha| P_k(x) dx ={ 8/45 if k == 2,
# \ 0 otherwise.
#
n = 5
moments = numpy.zeros(2*n)
moments[0] = 2.0/3.0
moments[2] = 8.0/45.0
_, _, b, c = \
orthopy.line.recurrence_coefficients.legendre(2*n, 'monic')
alpha, beta = orthopy.line.chebyshev_modified(moments, c)
assert numpy.all(abs(alpha) < tol)
assert numpy.all(
abs(beta - [2.0/3.0, 3.0/5.0, 4.0/35.0, 25.0/63.0, 16.0/99.0])
< tol
)
return
def calcmlmags(self, lightcurve):
"""
Returns a "lc.mags"-like array made using the ml-parameters.
It has the same size as lc.mags,and contains the microlensing to be added to them.
The lightcurve object is not changed !
For normal use,call getmags() from the lightcurve.
Idea : think about only returning the seasons mags to speed it up ? Not sure if reasonable,as no seasons defined outside ?
"""
jds = lightcurve.jds[self.season.indices] # Is this already a copy ? It seems so. So no need for an explicit copy().
# We do not need to apply shifts (i.e. getjds()),as anyway we "center" the jds.
# Old method :
if self.mltype == "poly":
refjd = np.mean(jds)
jds -= refjd # This is apparently safe,it does not shifts the lightcurves jds.
allmags = np.zeros(len(lightcurve.jds))
allmags[self.season.indices] = np.polyval(self.params, jds) # probably faster then +=
return allmags
# Legendre polynomials :
if self.mltype == "leg":
rjd = (np.max(jds) - np.min(jds))/2.0
cjd = (np.max(jds) + np.min(jds))/2.0
jds = (jds - cjd)/rjd
allmags = np.zeros(len(lightcurve.jds))
for (n, p) in enumerate(self.params):
allmags[self.season.indices] += p * ss.legendre(n)(jds)
return allmags
def factory(seasons, nparams, mltype="poly"):
"""
A factory function to create a microlensings object filled by seasonfct objects.
seasons is a list of season objects
nparams is an array or list of ints. "default" = one constant per season.
All parameters will be set to 0.0,and free to be optimized.
mltype = "poly" : simple polynomial microlensing,very stupid but fast,ok for degree <= 3
default type.
mltype = "leg" : legendre polynomials,very cLever but slow :-) These are fine for degree <= 25
Above deg 25,some numerical errors set in. Could perhaps be rewritten to make this faster.
Or implement in C !!!!
"""
#if nparams == "default":
# nparams = ones(len(seasons),dtype="int")
if len(nparams) != len(seasons):
raise RuntimeError, "Give as many nparams as they are seasons !"
mllist = []
for (season, n) in zip(seasons, nparams):
if n != 0:
p = np.zeros(n, dtype="float")
mask = p > -1.0
sfct = seasonfct(season, mltype, p, mask)
mllist.append(sfct)
return microlensing(mllist)
def test_eval(t, ref, tol=1.0e-14):
n = 5
p0, a, c = orthopy.line.recurrence_coefficients.legendre(
n, 'monic', symbolic=True
)
value = orthopy.line.evaluate_orthogonal_polynomial(t, p0, c)
assert value == ref
# Evaluating the Legendre polynomial in this way is rather unstable,so
# don't go too far with n.
approx_ref = numpy.polyval(legendre(n, monic=True), t)
assert abs(value - approx_ref) < tol
return
def test_eval_vec(t, c)
assert (value == ref).all()
# Evaluating the Legendre polynomial in this way is rather unstable, t)
assert (abs(value - approx_ref) < tol).all()
return
def test_logo():
import matplotlib.pyplot as plt
max_n = 6
moments = numpy.zeros(2*max_n)
moments[0] = 2.0 / 3.0
moments[2] = 8.0 / 45.0
for n in range(max_n):
_, c = orthopy.line.recurrence_coefficients.legendre(
2*n, standardization='p(1)=1'
)
alpha, beta = \
orthopy.line.chebyshev_modified(moments[:2*n], c)
orthopy.line.plot(1, len(alpha)*[1], -1.0, +1.0)
plt.xlim(-1, +1)
plt.ylim(-2, +2)
plt.grid()
plt.tick_params(
axis='both',
which='both',
left='off',
labelleft='off',
bottom='off',
labelbottom='off',
)
plt.gca().set_aspect(0.25)
plt.show()
# plt.savefig('logo.png',transparent=True)
return
def test_compute_moments():
moments = orthopy.line.compute_moments(lambda x: 1, -1, +1, 5)
assert (
moments == [2, 0, sympy.Rational(2, 3), 5)]
).all()
moments = orthopy.line.compute_moments(
lambda x: 1, 5,
polynomial_class=orthopy.line.legendre
)
assert (moments == [2, 0]).all()
# Example from Gautschi's "How to and how not to" article
moments = orthopy.line.compute_moments(
lambda x: sympy.exp(-x**3/3),
0, sympy.oo,
5
)
reference = [
3**sympy.Rational(n-2, 3) * sympy.gamma(sympy.Rational(n+1, 3))
for n in range(5)
]
assert numpy.all([
sympy.simplify(m - r) == 0
for m, r in zip(moments, reference)
])
return
def test_stieltjes():
alpha0, beta0 = orthopy.line.stieltjes(lambda t: 1, 5)
_, alpha1, beta1 = \
orthopy.line.recurrence_coefficients.legendre(5, 'monic')
assert (alpha0 == alpha1).all()
assert (beta0 == beta1).all()
return
# def test_expt3():
# '''Full example from Gautschi's "How to and how not to" article.
# '''
# # moments = orthopy.line.compute_moments(
# # lambda x: sympy.exp(-x**3/3),
# # 0,sympy.oo,
# # 31
# # )
# # print(moments)
# # alpha,beta = orthopy.line.chebyshev(moments)
#
# alpha,beta = orthopy.line.stieltjes(
# lambda x: sympy.exp(-x**3/3),
# 0,
# 5
# )
# print(alpha)
# print(beta)
# return
def test_xk(k):
n = 10
moments = orthopy.line.compute_moments(lambda x: x**k, 2*n)
alpha, beta = orthopy.line.chebyshev(moments)
assert (alpha == 0).all()
assert beta[0] == moments[0]
assert beta[1] == sympy.Rational(k+1, k+3)
assert beta[2] == sympy.Rational(4, (k+5) * (k+3))
orthopy.line.schemes.custom(
numpy.array([sympy.N(a) for a in alpha], dtype=float),
numpy.array([sympy.N(b) for b in beta],
mode='numpy'
)
# a,b = \
# orthopy.line.recurrence_coefficients.legendre(
# 2*n,mode='sympy'
# )
# moments = orthopy.line.compute_moments(
# lambda x: x**2,-1,+1,2*n,
# polynomial_class=orthopy.line.legendre
# )
# alpha,beta = orthopy.line.chebyshev_modified(moments,a,b)
# points,weights = orthopy.line.schemes.custom(
# numpy.array([sympy.N(a) for a in alpha],dtype=float),
# numpy.array([sympy.N(b) for b in beta],dtype=float)
# )
return
def Legendre_matrix(N, ctheta):
r"""Matrix of weighted Legendre polynominals.
Computes a matrix of weighted Legendre polynominals up to order N for
the given angles
.. math::
L_n(\theta) = \frac{2n+1}{4 \pi} P_n(\theta)
Parameters
----------
N : int
Maximum order.
ctheta : (Q,) array_like
Angles.
Returns
-------
Lmn : ((N+1),Q) numpy.ndarray
Matrix containing Legendre polynominals.
"""
if ctheta.ndim == 0:
M = 1
else:
M = len(ctheta)
Lmn = np.zeros([N+1, M], dtype=complex)
for n in range(N+1):
Lmn[n, :] = (2*n+1)/(4*np.pi) * np.polyval(special.legendre(n), ctheta)
return Lmn
def grid_gauss(n):
""" Gauss-Legendre sampling points on sphere.
According to (cf. Rafaely book,sec.3.3)
Parameters
----------
n : int
Maximum order.
Returns
-------
azi : array_like
Azimuth.
elev : array_like
Elevation.
weights : array_like
Quadrature weights.
"""
azi = np.linspace(0, 2*np.pi, 2*n+2, endpoint=False)
x, weights = np.polynomial.legendre.leggauss(n+1)
elev = np.arccos(x)
azi = np.tile(azi, n+1)
elev = np.repeat(elev, 2*n+2)
weights = np.repeat(weights, 2*n+2)
weights *= np.pi / (n+1) # sum(weights) == 4pi
return azi, elev, weights
def smooth(self, lightcurve):
"""
Only for plotting purposes : returns jds,mlmagshifts,and refmags with a tight and regular sampling,
over the range given by the season.
Note that this time we are interested in the acutal shifted jds,so that it looks right when plotted !
We return arrays that can directly be plotted to illustrate the microlensing.
Todo : return refmag,not refmags !
"""
jds = lightcurve.getjds()[self.season.indices]
# Old method :
if self.mltype == "poly":
refjd = np.mean(jds)
jds -= refjd
refmag = np.median(lightcurve.getmags()) # So the reference magnitude is evaluated from the entire lightcurve.
# Independent on seasons.
smoothtime = np.linspace(jds[0], jds[-1], 50)
smoothml = np.polyval(self.params, smoothtime)
smoothtime += refjd # important,to get the time back at the right place.
refmags = np.zeros(50) + refmag
return {"jds":smoothtime, "ml":smoothml, "refmags":refmags}
# Legendre polynomials :
if self.mltype == "leg":
rjd = (np.max(jds) - np.min(jds))/2.0
cjd = (np.max(jds) + np.min(jds))/2.0
jds = (jds - cjd)/rjd
refmag = np.median(lightcurve.getmags()) # So the reference magnitude is evaluated from the entire lightcurve.
# Independent on seasons.
smoothtime = np.linspace(-1, 1, 300)
smoothml = np.zeros(len(smoothtime))
for (n, p) in enumerate(self.params):
smoothml += p * ss.legendre(n)(smoothtime)
smoothtime = smoothtime * rjd + cjd # important,to get the time back at the right place.
refmags = np.zeros(300) + refmag
return {"jds":smoothtime, "refmags":refmags}