Python 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}