Interpolation

import numpy as np
from typing import Callable, List, Tuple


class InterpolationToolkit:
    """Comprehensive interpolation toolkit."""

    @staticmethod
    def lagrange(x_data: np.ndarray, y_data: np.ndarray,
                 x: np.ndarray) -> np.ndarray:
        """Lagrange interpolation."""
        n = len(x_data)
        result = np.zeros_like(x, dtype=float)

        for j in range(n):
            L_j = np.ones_like(x, dtype=float)
            for i in range(n):
                if i != j:
                    L_j *= (x - x_data[i]) / (x_data[j] - x_data[i])
            result += y_data[j] * L_j

        return result

    @staticmethod
    def divided_differences(x_data: np.ndarray,
                           y_data: np.ndarray) -> np.ndarray:
        """Compute divided difference table."""
        n = len(x_data)
        dd = np.zeros((n, n))
        dd[:, 0] = y_data

        for j in range(1, n):
            for i in range(n - j):
                dd[i, j] = (dd[i+1, j-1] - dd[i, j-1]) / (x_data[i+j] - x_data[i])

        return dd

    @staticmethod
    def newton(x_data: np.ndarray, y_data: np.ndarray,
               x: np.ndarray) -> np.ndarray:
        """Newton interpolation using Horner's method."""
        dd = InterpolationToolkit.divided_differences(x_data, y_data)
        n = len(x_data)

        result = np.full_like(x, dd[0, n-1], dtype=float)
        for k in range(n-2, -1, -1):
            result = result * (x - x_data[k]) + dd[0, k]

        return result

    @staticmethod
    def natural_cubic_spline(x_data: np.ndarray, y_data: np.ndarray,
                             x_eval: np.ndarray) -> np.ndarray:
        """Natural cubic spline interpolation."""
        n = len(x_data) - 1
        h = np.diff(x_data)

        # Set up tridiagonal system for moments
        A = np.zeros((n+1, n+1))
        b = np.zeros(n+1)
        A[0, 0] = 1
        A[n, n] = 1

        for i in range(1, n):
            A[i, i-1] = h[i-1]
            A[i, i] = 2 * (h[i-1] + h[i])
            A[i, i+1] = h[i]
            b[i] = 6 * ((y_data[i+1] - y_data[i]) / h[i] -
                        (y_data[i] - y_data[i-1]) / h[i-1])

        m = np.linalg.solve(A, b)

        # Evaluate
        results = np.zeros_like(x_eval, dtype=float)
        for j, x in enumerate(x_eval):
            i = np.searchsorted(x_data, x, side='right') - 1
            i = max(0, min(i, n-1))

            dx = x - x_data[i]
            hi = h[i]

            A_c = (m[i+1] - m[i]) / (6 * hi)
            B_c = m[i] / 2
            C_c = (y_data[i+1] - y_data[i]) / hi - (m[i+1] + 2*m[i]) * hi / 6
            D_c = y_data[i]

            results[j] = A_c * dx**3 + B_c * dx**2 + C_c * dx + D_c

        return results

    @staticmethod
    def chebyshev_nodes(a: float, b: float, n: int) -> np.ndarray:
        """Generate Chebyshev nodes on [a, b]."""
        k = np.arange(n)
        nodes = np.cos((2*k + 1) * np.pi / (2*n))
        return 0.5 * (b - a) * nodes + 0.5 * (b + a)

    @staticmethod
    def linear_spline(x_data: np.ndarray, y_data: np.ndarray,
                      x_eval: np.ndarray) -> np.ndarray:
        """Linear spline (piecewise linear interpolation)."""
        return np.interp(x_eval, x_data, y_data)


def demonstrate_interpolation():
    """Comprehensive interpolation demonstration."""
    print("=" * 60)
    print("INTERPOLATION METHOD COMPARISON")
    print("=" * 60)

    toolkit = InterpolationToolkit()

    # Test function
    f = lambda x: np.sin(x) + 0.5 * np.cos(2 * x)
    x_data = np.linspace(0, 2*np.pi, 8)
    y_data = f(x_data)
    x_fine = np.linspace(0, 2*np.pi, 200)
    y_true = f(x_fine)

    print("\n--- Method Comparison (8 data points) ---")
    print(f"{'Method':<25} {'Max Error':<14}")
    print("-" * 39)

    y_lagrange = toolkit.lagrange(x_data, y_data, x_fine)
    print(f"{'Lagrange':<25} {np.max(np.abs(y_lagrange - y_true)):<14.6e}")

    y_newton = toolkit.newton(x_data, y_data, x_fine)
    print(f"{'Newton':<25} {np.max(np.abs(y_newton - y_true)):<14.6e}")

    y_spline = toolkit.natural_cubic_spline(x_data, y_data, x_fine)
    print(f"{'Cubic Spline':<25} {np.max(np.abs(y_spline - y_true)):<14.6e}")

    y_linear = toolkit.linear_spline(x_data, y_data, x_fine)
    print(f"{'Linear Spline':<25} {np.max(np.abs(y_linear - y_true)):<14.6e}")

    # Chebyshev nodes comparison
    print("\n--- Chebyshev vs Equally-Spaced (Runge function) ---")
    runge = lambda x: 1 / (1 + 25*x**2)

    for n in [5, 10, 15]:
        # Equally spaced
        x_eq = np.linspace(-1, 1, n)
        y_eq = runge(x_eq)
        y_eq_interp = toolkit.lagrange(x_eq, y_eq, x_fine)
        err_eq = np.max(np.abs(y_eq_interp - runge(x_fine)))

        # Chebyshev
        x_cheb = toolkit.chebyshev_nodes(-1, 1, n)
        y_cheb = runge(x_cheb)
        y_cheb_interp = toolkit.lagrange(x_cheb, y_cheb, x_fine)
        err_cheb = np.max(np.abs(y_cheb_interp - runge(x_fine)))

        print(f"  n={n:>2}: Equally-spaced error = {err_eq:.6f}, "
              f"Chebyshev error = {err_cheb:.6f}")


demonstrate_interpolation()