Source code for boundlab.ops.legendre

r'''Legendre polynomial bases and basis-to-monomial conversion.

Legendre polynomials :math:`P_i` are orthogonal on :math:`[-1, 1]`
(:math:`\int P_i P_j = \frac{2}{2i+1}\delta_{ij}`) and satisfy Bonnet's
recurrence :math:`(i + 1) P_{i+1}(x) = (2i + 1) x\, P_i(x) - i P_{i-1}(x)`,
which :func:`generate_legendre_polynomials` unrolls.  A truncated Legendre
series is the :math:`L^2`-optimal polynomial approximation on the interval,
which is why the :mod:`boundlab.polysp.legendre` relaxations expand in this
basis before converting to monomial coefficients.
'''

import torch

from boundlab import utils
from boundlab.utils import Polynomial



def generate_legendre_polynomials(count: int) -> list[Polynomial]:
    """Build the first ``count`` standard Legendre polynomials."""
    if count == 1:
        return [Polynomial([1.0])]

    polynomials = [[1.0], [0.0, 1.0]]
    for degree in range(1, count - 1):
        denominator = degree + 1
        current = polynomials[-1]
        previous = polynomials[-2]
        following = [0.0] * (degree + 2)

        x_scale = (2 * degree + 1) / denominator
        for power, coefficient in enumerate(current):
            following[power + 1] += x_scale * coefficient

        previous_scale = degree / denominator
        for power, coefficient in enumerate(previous):
            following[power] -= previous_scale * coefficient

        polynomials.append(following)

    return [Polynomial(coeffs) for coeffs in polynomials]





def legendre_polynomial(
    li: list[torch.Tensor], hw: torch.Tensor
) -> utils.Polynomial:
    """Convert ``sum_i li[i] * hw**i * P_i(x/hw)`` to monomial coefficients."""
    basis = generate_legendre_polynomials(len(li))
    polyc = utils.sum(
        legendre.apply_ax_reversed(hw) * coefficient
        for coefficient, legendre in zip(li, basis)
    )

    return polyc


__all__ = [
    "generate_legendre_polynomials",
    "legendre_polynomial",
]