Merge pull request #134 from jorenham/reworked-bounds

Reworked `diagnostic.l_ratio_bounds`

lmo/_poly.py

@@ -9,15 +9,15 @@
 
 __all__ = (
     'eval_sh_jacobi',
+    'peaks_jacobi',
+    'arg_extrema_jacobi',
+    'extrema_jacobi',
     'jacobi',
     'jacobi_series',
     'roots',
-    'extrema',
-    'minima',
-    'maxima',
 )
 
-from typing import Any, TypeVar, cast, overload
+from typing import TypeVar, cast, overload
 
 import numpy as np
 import numpy.polynomial as npp
@@ -101,6 +101,190 @@ def eval_sh_jacobi(
     return scs.eval_jacobi(n, a, b, u)  # type: ignore
 
 
+def peaks_jacobi(n: int, a: float, b: float) -> npt.NDArray[np.float64]:
+    r"""
+    Finds the \( x \in [-1, 1] \) s.t.
+    \( /frac{\dd{\shjacobi{n}{a}{b}{x}}}{\dd{x}} = 0 \) of a Jacobi polynomial,
+    which includes the endpoints \( x \in \{-1, 1\} \). I.e. the locations of
+    the peaks.
+
+    The Jacobi polynomials with order \( n \) have \( n + 1 \) peaks.
+
+    Examples:
+        For \( n = 0 \) there is only one "peak", since
+        \( \jacobi{0}{a}{b}{x} = 1 \):
+
+        >>> peaks_jacobi(0, 0, 0)
+        array([0.])
+
+        The `0` is arbitrary; all \( x \in [0, 1] \) evaluate to the same
+        constant \( 1 \).
+
+        For \( n = 1 \), it is a positive linear function, so the peaks
+        are exactly the endpoints, and do not depend on \( a \) or \( b \):
+
+        >>> peaks_jacobi(1, 0, 0)
+        array([-1., 1.])
+        >>> peaks_jacobi(1, 3.14, -1 / 12)
+        array([-1., 1.])
+
+        For \( n > 1 \), the effects of the choices for \( a \) and \( b \)
+        become apparent, e.g. for \( n = 4 \):
+
+        >>> peaks_jacobi(4, 0, 0).round(5)
+        array([-1.     , -0.65465,  0.     ,  0.65465,  1.     ])
+        >>> peaks_jacobi(4, 0, 1).round(5)
+        array([-1.     , -0.50779,  0.1323 ,  0.70882,  1.     ])
+        >>> peaks_jacobi(4, 1, 0).round(5)
+        array([-1.     , -0.70882, -0.1323 ,  0.50779,  1.     ])
+        >>> peaks_jacobi(4, 1, 1).round(5)
+        array([-1.     , -0.57735,  0.     ,  0.57735,  1.     ])
+        >>> peaks_jacobi(4, 2.5, 2.5)
+        array([-1. , -0.5,  0. ,  0.5,  1. ])
+        >>> peaks_jacobi(4, 10, 10).round(5)
+        array([-1.     , -0.33333,  0.     ,  0.33333,  1.     ])
+    """
+    if n == 0:
+        # constant; any x is a "peak"; so take the "middle ground"
+        return np.array([0.])
+    if n == 1:
+        # linear; the peaks are only at the ends
+        return np.array([-1., 1.])
+
+    # otherwise, peaks are at the ends, and at the roots of the derivative
+    x = np.empty(n + 1)
+    x[0] = -1
+    x[1:-1] = scs.roots_jacobi(n - 1, a + 1, b + 1)[0]  # type: ignore
+    x[-1] = 1
+
+    return np.round(x, 15) + 0.0  # cleanup of numerical noise
+
+def arg_extrema_jacobi(n: int, a: float, b: float) -> tuple[float, float]:
+    r"""
+    Find the \( x \) of the minimum and maximum values of a Jacobi polynomial
+    on \( [-1, 1] \).
+
+    Note:
+        There can be multiple \( x \) that share the same extremum, but only
+        one of them is returned, which for \( n > 0 \) is the smallest (first)
+        one.
+
+    Examples:
+        For \( n = 1 \), the Jacobi polynomials are positive linear function
+        (i.e. a straight line), so the minimum and maximum are the left and
+        right endpoints of the domain.
+
+        >>> arg_extrema_jacobi(1, 0, 0)
+        (-1.0, 1.0)
+        >>> arg_extrema_jacobi(1, 3.14, -1 / 12)
+        (-1.0, 1.0)
+
+        The 2nd degree Jacobi polynomial is a positive parabola, with one
+        unique minimum, and maxima at \( -1 \) and/or \( 1 \).
+        When \( a == b \), the parabola is centered within the domain, and
+        has maxima at both \( x = -1 \) and \( x=1 \). For the sake of
+        simplicity, only one (the first) value is returned in such cases:
+
+        >>> arg_extrema_jacobi(2, 0, 0)
+        (0.0, -1.0)
+        >>> arg_extrema_jacobi(2, 42, 42)
+        (0.0, -1.0)
+
+        Conversely, when \( a \neq b \), the parabola is "shifted" so that
+        there is only one global maximum:
+
+        >>> arg_extrema_jacobi(2, 0, 1)
+        (0.2, -1.0)
+        >>> arg_extrema_jacobi(2, 1, 0)
+        (-0.2, 1.0)
+        >>> arg_extrema_jacobi(2, 10, 2)
+        (-0.5, 1.0)
+
+    """
+    x = peaks_jacobi(n, a, b)
+    p = eval_sh_jacobi(n, a, b, (x + 1) / 2)
+
+    return x[np.argmin(p)], x[np.argmax(p)]
+
+def extrema_jacobi(n: int, a: float, b: float) -> tuple[float, float]:
+    r"""
+    Find the global minimum and maximum values of a (shifted) Jacobi
+    polynomial on \( [-1, 1] \) (or equivalently \( [0, 1] \) if shifted).
+
+    Examples:
+        With \( n \), \( \jacobi{0}{a}{b}{x} = 1 \), so there is only one
+        "extremum":
+
+        >>> extrema_jacobi(0, 0, 0)
+        (1, 1)
+        >>> extrema_jacobi(0, 3.14, -1 / 12)
+        (1, 1)
+
+        With \( n = 1 \), the extrema are always at \( -1 \) and \( 1 \),
+        but their values depend on \( a \) and \( b \):
+
+        >>> extrema_jacobi(1, 0, 0)
+        (-1.0, 1.0)
+        >>> extrema_jacobi(1, 0, 1)
+        (-2.0, 1.0)
+        >>> extrema_jacobi(1, 1, 0)
+        (-1.0, 2.0)
+        >>> extrema_jacobi(1, 1, 1)
+        (-2.0, 2.0)
+
+        For \( n = 2 \) (a parabola), the relation between \( a, b \)
+        and the extrema isn't as obvious:
+
+        >>> extrema_jacobi(2, 0, 0)
+        (-0.5, 1.0)
+        >>> extrema_jacobi(2, 0, 4)
+        (-0.75, 15.0)
+        >>> extrema_jacobi(2, 4, 0)
+        (-0.75, 15.0)
+        >>> extrema_jacobi(2, 4, 4)
+        (-1.5, 15.0)
+
+        With \( n = 3 \), the extrema appear to behave very predictable:
+
+        >>> extrema_jacobi(3, 0, 0)
+        (-1.0, 1.0)
+        >>> extrema_jacobi(3, 0, 1)
+        (-4.0, 1.0)
+        >>> extrema_jacobi(3, 1, 0)
+        (-1.0, 4.0)
+        >>> extrema_jacobi(3, 1, 1)
+        (-4.0, 4.0)
+
+        However, if we keep \( a \) fixed at \( 0 \), and increase \( b \),
+        the plot-twist emerges:
+
+        >>> extrema_jacobi(3, 0, 2)
+        (-10.0, 1.0)
+        >>> extrema_jacobi(3, 0, 3)
+        (-20.0, 1.0)
+        >>> extrema_jacobi(3, 0, 4)
+        (-35.0, 1.13541...)
+        >>> extrema_jacobi(3, 0, 5)
+        (-56.0, 1.25241...)
+
+        Looking at the corresponding \( x \) can help to understand the
+        "movement" of the maximum.
+
+        >>> arg_extrema_jacobi(3, 0, 2)
+        (-1.0, 1.0)
+        >>> arg_extrema_jacobi(3, 0, 3)
+        (-1.0, 0.0)
+        >>> arg_extrema_jacobi(3, 0, 4)
+        (-1.0, 0.09449...)
+        >>> arg_extrema_jacobi(3, 0, 5)
+        (-1.0, 0.17287...)
+
+    """
+    x = peaks_jacobi(n, a, b)
+    p = eval_sh_jacobi(n, a, b, (x + 1) / 2)
+    return cast(float, np.min(p)), cast(float, np.max(p))
+
+
 def _jacobi_coefs(n: int, a: float, b: float) -> npt.NDArray[np.float64]:
     p_n: np.poly1d
     p_n = scs.jacobi(n, a, b)  # type: ignore [reportUnknownMemberType]
@@ -182,41 +366,3 @@ def roots(
         return x[(x >= a) & (x <= b)]
 
     return x
-
-
-def integrate(p: PolySeries, /, a: float | None = None) -> PolySeries:
-    r"""Calculate the anti-derivative: $P(x) = \int_a^x p(u) \, du$."""
-    return p.integ(lbnd=p.domain[0] if a is None else a)
-
-
-def extrema(
-    p: PolySeries,
-    /,
-    outside: bool = False,
-) -> npt.NDArray[np.inexact[Any]]:
-    """Return the $x$ in the domain of $p$, where $p'(x) = 0$."""
-    return roots(p.deriv(), outside=outside)
-
-
-def minima(
-    p: PolySeries,
-    /,
-    outside: bool = False,
-) -> npt.NDArray[np.inexact[Any]]:
-    """
-    Return the $x$ in the domain of $p$, where $p'(x) = 0$ and $p''(x) > 0$.
-    """  # noqa: D200
-    x = extrema(p, outside=outside)
-    return x[p.deriv(2)(x) > 0] if len(x) else x
-
-
-def maxima(
-    p: PolySeries,
-    /,
-    outside: bool = False,
-) -> npt.NDArray[np.inexact[Any]]:
-    """
-    Return the $x$ in the domain of $p$, where $p'(x) = 0$ and $p''(x) < 0$.
-    """  # noqa: D200
-    x = extrema(p, outside=outside)
-    return x[p.deriv(2)(x) < 0] if len(x) else x