我们从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, b, 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 feature_transform(X, mode='polynomial', degree=1): poly = PolynomialFeatures(degree) process_X = poly.fit_transform(X) if mode == 'legendre': lege = legendre(degree) process_X = lege(process_X) return process_X
def test_eval(t, ref, tol=1.0e-14): n = 5 p0, a, b, c = orthopy.line.recurrence_coefficients.legendre( n, 'monic', symbolic=True ) value = orthopy.line.evaluate_orthogonal_polynomial(t, p0, a, b, 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, ref, tol=1.0e-14): n = 5 p0, a, b, c = orthopy.line.recurrence_coefficients.legendre( n, 'monic', symbolic=True ) value = orthopy.line.evaluate_orthogonal_polynomial(t, p0, a, b, c) assert (value == ref).all() # 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).all() return
def test_clenshaw(tol=1.0e-14): n = 5 _, _, alpha, beta = \ orthopy.line.recurrence_coefficients.legendre(n, 'monic') t = 1.0 a = numpy.ones(n+1) value = orthopy.line.clenshaw(a, alpha, beta, t) ref = math.fsum([ numpy.polyval(legendre(i, monic=True), t) for i in range(n+1)]) assert abs(value - ref) < tol 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): _, _, b, c = orthopy.line.recurrence_coefficients.legendre( 2*n, standardization='p(1)=1' ) alpha, beta = \ orthopy.line.chebyshev_modified(moments[:2*n], b, c) orthopy.line.plot(1, len(alpha)*[1], alpha, beta, -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), 0, sympy.Rational(2, 5)] ).all() moments = orthopy.line.compute_moments( lambda x: 1, -1, +1, 5, polynomial_class=orthopy.line.legendre ) assert (moments == [2, 0, 0, 0, 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, -1, +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, sympy.oo, # 5 # ) # print(alpha) # print(beta) # return
def test_xk(k): n = 10 moments = orthopy.line.compute_moments(lambda x: x**k, -1, +1, 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], dtype=float), 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, "ml":smoothml, "refmags":refmags}