feat(interpolation): add Lagrange interpolation helpers

docs/misc.rst

@@ -22,6 +22,11 @@ Finite Difference Approximations
 
 .. automodule:: pycaputo.finite_difference
 
+Interpolation
+-------------
+
+.. automodule:: pycaputo.interpolation
+
 Generating Functions
 --------------------
 

pycaputo/interpolation.py

@@ -0,0 +1,146 @@
+# SPDX-FileCopyrightText: 2023 Alexandru Fikl <alexfikl@gmail.com>
+# SPDX-License-Identifier: MIT
+
+from __future__ import annotations
+
+import math
+from dataclasses import dataclass
+from functools import cached_property
+from typing import Any, Iterator
+
+import numpy as np
+
+from pycaputo.utils import Array
+
+
+@dataclass(frozen=True)
+class InterpStencil:
+    r"""Approximation of a function by Lagrange poylnomials on a uniform grid.
+
+    .. math::
+
+        f(x_m) = \sum_{k \in \text{offsets}} f_{m + k} \ell_k(x_m)
+
+    where :math:`\ell_k` is the :math:`k`-th Lagrange polynomial. The approximation
+    is of the :attr:`order` of the Lagrange polynomials.
+
+    Note that Lagrange polynomials are notoriously ill-conditioned on uniform
+    grids (Runge phenomenon), so they should not be used with high-order.
+    """
+
+    #: Coefficients used in the stencil.
+    coeffs: Array
+    #: Offsets around the centered :math:`0` used in the stencil.
+    offsets: Array
+    #: Point at which the interpolation was evaluated.
+    x: float
+
+    @property
+    def order(self) -> int:
+        """Order of the Lagrange polynomial."""
+        return self.coeffs.size
+
+    @cached_property
+    def padded_coeffs(self) -> Array:
+        """Padded coefficients that are symmetric around the :math:`0`
+        index and can be easily applied as a convolution.
+        """
+        n = np.max(np.abs(self.offsets))
+        coeffs = np.zeros(2 * n + 1, dtype=self.coeffs.dtype)
+        coeffs[n + self.offsets] = self.coeffs
+
+        return coeffs
+
+    @cached_property
+    def trunc(self) -> float:
+        """Truncation error of the interpolation."""
+        return determine_truncation_error(self.offsets, self.x)
+
+
+def apply_interpolation(s: InterpStencil, f: Array) -> Array:
+    """Apply the stencil to a function *f*.
+
+    Note that only interior points are correctly computed. Any boundary
+    points will contain invalid values.
+
+    :returns: the stencil applied to the function *f*.
+    """
+    a = s.padded_coeffs.astype(f.dtype)
+    return np.convolve(f, a, mode="same")
+
+
+def determine_truncation_error(offsets: Array, x: float, h: float = 1.0) -> float:
+    r"""Approximate the truncation error of the Lagrange interpolation.
+
+    .. math::
+
+        f(x) - \sum_{k \in \text{offsets}} \ell_k(x) f_{m + k} =
+        c \frac{\mathrm{d}^n f}{\mathrm{d} x^n}(\xi),
+
+    where the constant :math:`c` is approxiated as the truncation error. The
+    Error itself also depends on the derivatives of :math:`f` at a undetermined
+    point :math:`\xi` in the interpolation interval.
+
+    :arg h: grid size on the physical grid, which is necessary for an accurate
+        estimate.
+    """
+
+    c = h**offsets.size / math.factorial(offsets.size) * np.prod(x - offsets)
+    return float(c)
+
+
+def wandering(
+    n: int, wanderer: int | bool | float = 1.0, landscape: int | bool | float = 0.0
+) -> Iterator[Array]:
+    for i in range(n):
+        yield np.array([landscape] * i + [wanderer] + [landscape] * (n - i - 1))
+
+
+def make_lagrange_approximation(
+    bounds: tuple[int, int],
+    x: int | float = 0.5,
+    *,
+    dtype: np.dtype[Any] | None = None,
+) -> InterpStencil:
+    r"""Construct interpolation coefficients at *x* on a uniform grid.
+
+    :arg bounds: inclusive left and right bounds on the stencil around a point
+        :math:`x_i`. For example, ``(-1, 2)`` defines the 4 point stencil
+        :math:`\{x_{i - 1}, x_i, x_{i + 1}, x_{i + 2}\}`.
+    :arg x: point (in index space) at which to evaluate the function, i.e. to
+        evaluate at :math:`x_{i + 1/2}` this should be :math:`1/2`, to evaluate at
+        :math:`x_{i - 3/2}` this should be :math:`-3/2`, etc.
+    """
+    if len(bounds) != 2:
+        raise ValueError(f"Stencil bounds are invalid: {bounds}")
+
+    if bounds[0] > bounds[1]:
+        bounds = (bounds[1], bounds[0])
+
+    if bounds[0] > 0 or bounds[1] < 0:
+        raise ValueError(f"Bounds must be (smaller <= 0, bigger >= 0): {bounds}")
+
+    if dtype is None:
+        dtype = np.dtype(np.float64)
+    dtype = np.dtype(dtype)
+
+    # evaluate Lagrange polynomials
+    offsets = np.arange(bounds[0], bounds[1] + 1)
+
+    x = float(x)
+    if x.is_integer() and bounds[0] <= x <= bounds[1]:
+        coeffs = np.zeros(offsets.size, dtype=dtype)
+        coeffs[abs(bounds[0]) + int(x)] = 1.0
+    else:
+        # NOTE: this evaluates l_i(x) = prod((x - x_m) / (x_i - x_m), i != m)
+        coeffs = np.array(
+            [
+                np.prod((x - offsets[not_n]) / (offsets[n] - offsets[not_n]))
+                for n, not_n in enumerate(
+                    wandering(offsets.size, wanderer=False, landscape=True)
+                )
+            ],
+            dtype=dtype,
+        )
+
+    return InterpStencil(coeffs=coeffs, offsets=offsets, x=x)

tests/test_finite_difference.py

@@ -15,15 +15,17 @@
     modified_wavenumber,
 )
 from pycaputo.logging import get_logger
-from pycaputo.utils import EOCRecorder, savefig, set_recommended_matplotlib
+from pycaputo.utils import savefig, set_recommended_matplotlib
 
 logger = get_logger("pycaputo.test_finite_difference")
 set_recommended_matplotlib()
 
 # {{{ test_finite_difference_taylor
 
 
-def finite_difference_convergence(d: DiffStencil) -> EOCRecorder:
+def finite_difference_convergence(d: DiffStencil) -> float:
+    from pycaputo.utils import EOCRecorder
+
     eoc = EOCRecorder()
 
     s = np.s_[abs(d.offsets[0]) + 1 : -abs(d.offsets[-1]) - 1]
@@ -40,7 +42,8 @@ def finite_difference_convergence(d: DiffStencil) -> EOCRecorder:
         error = np.linalg.norm(df_dx[s] - num_df_dx[s]) / np.linalg.norm(df_dx[s])
         eoc.add_data_point(h, error)
 
-    return eoc
+    logger.info("\n%s", eoc)
+    return eoc.estimated_order
 
 
 def test_finite_difference_taylor_stencil(*, visualize: bool = False) -> None:
@@ -102,9 +105,8 @@ def test_finite_difference_taylor_stencil(*, visualize: bool = False) -> None:
         assert np.allclose(s.trunc.error, coefficient)
         assert s.trunc.order == order
 
-        eoc = finite_difference_convergence(s)
-        logger.info("\n%s", eoc)
-        assert eoc.estimated_order >= order - 0.25
+        estimated_order = finite_difference_convergence(s)
+        assert estimated_order >= order - 0.25
 
         if visualize:
             part = np.real if s.derivative % 2 == 0 else np.imag

tests/test_interpolation.py

@@ -0,0 +1,75 @@
+# SPDX-FileCopyrightText: 2023 Alexandru Fikl <alexfikl@gmail.com>
+# SPDX-License-Identifier: MIT
+
+from __future__ import annotations
+
+import numpy as np
+
+from pycaputo.interpolation import (
+    InterpStencil,
+    apply_interpolation,
+    make_lagrange_approximation,
+)
+from pycaputo.logging import get_logger
+from pycaputo.utils import set_recommended_matplotlib
+
+logger = get_logger("pycaputo.test_interpolation")
+set_recommended_matplotlib()
+
+
+# {{{ test_interpolation_lagrange
+
+
+def interpolation_convergence(s: InterpStencil) -> float:
+    from pycaputo.utils import EOCRecorder
+
+    eoc = EOCRecorder()
+
+    k = np.s_[abs(s.offsets[0]) + 1 : -abs(s.offsets[-1]) - 1]
+    for n in [32, 64, 128, 256, 512]:
+        theta = np.linspace(0.0, 2.0 * np.pi, n, dtype=s.coeffs.dtype)
+        theta_m = (theta[1:] + theta[:-1]) / 2.0
+        h = theta[1] - theta[0]
+
+        f = np.sin(theta)
+        fhat = apply_interpolation(s, f)[:-1]
+        f_ref = np.sin(theta_m)
+
+        error = np.linalg.norm(fhat[k] - f_ref[k]) / np.linalg.norm(f_ref[k])
+        eoc.add_data_point(h, error)
+
+    logger.info("\n%s", eoc)
+    return eoc.estimated_order
+
+
+def test_interpolation_lagrange(*, visualize: bool = False) -> None:
+    stencils = [
+        # (
+        #     make_lagrange_approximation((0, 1), 0.5),
+        #     np.array([1 / 2, 1 / 2]),
+        #     2,
+        # ),
+        # (
+        #     make_lagrange_approximation((-1, 1), 0.5),
+        #     np.array([-1 / 8, 6 / 8, 3 / 8]),
+        #     3,
+        # ),
+        (
+            make_lagrange_approximation((-1, 2), -0.5),
+            np.array([5 / 16, 15 / 16, -5 / 16, 1 / 16]),
+            4,
+        )
+    ]
+
+    for s, a, order in stencils:
+        logger.info("stencil:\n%r", s)
+
+        assert np.allclose(np.sum(s.coeffs), 1.0)
+        assert np.allclose(s.coeffs, np.array(a, dtype=s.coeffs.dtype))
+        assert s.order == order
+
+        estimated_order = interpolation_convergence(s)
+        assert estimated_order >= order - 0.25
+
+
+# }}}