<jupyter_start><jupyter_text>This notebook is part of https://github.com/AudioSceneDescriptionFormat/splines, see also https://splines.readthedocs.io/.  Lagrange InterpolationBefore diving into splines,let's have a look at an arguably simpler interpolation method using polynomials:[Lagrange interpolation](https://en.wikipedia.org/wiki/Lagrange_polynomial).This is easy to implement, but as we will see,it has quite severe limitations,which will motivate us to look into splines later.<jupyter_code>import matplotlib.pyplot as plt
import numpy as np<jupyter_output><empty_output><jupyter_text>One-dimensional ExampleAssume we have $N$ time instants $t_i$, with $0 \le i < N$:<jupyter_code>ts = -1.5, 0.5, 1.7, 3, 4<jupyter_output><empty_output><jupyter_text>... and for each time instant we are given an associated value $x_i$:<jupyter_code>xs = 2, -1, 1.3, 3.14, 1<jupyter_output><empty_output><jupyter_text>Our task is now to find a function that yieldsthe given $x_i$ values for the given times $t_i$ andsome "reasonable" interpolated values when evaluated at time values in-between.The idea of Lagrange interpolation is to create a separate polynomial$\ell_j(t)$for each of the $N$ given time instants,which will be weighted by the associated $x_j$.The final interpolation function is the weighted sum of these $N$ polynomials:\begin{equation*}L(t) = \sum_{j=0}^{N-1} x_j \ell_j(t)\end{equation*} In order for this to actually work,the polynomials must fulfill the following requirements:* Each polynomial must yield $1$ when evaluated at its associated time $t_i$.* Each polynomial must yield $0$ at all other instances in the set of given times. To satisfy the second point,let's create a product with a term for each of the relevant timesand make each of those factors vanish when evaluated at their associated time.As an example we look at the basis for $t_3 = 3$:<jupyter_code>def maybe_polynomial_3(t):
    t = np.asarray(t)
    return (t - (-1.5)) * (t - 0.5) * (t - 1.7) * (t - 4)
maybe_polynomial_3(ts)<jupyter_output><empty_output><jupyter_text>As we can see, this indeed fulfills the second requirement.Note that we were given 5 time instants,but we need only 4 product terms(corresponding to the 4 roots of the polynomial).Mathematically, this could be written as\begin{equation*}\ell_j(t) \overset{?}{=} \prod_{\substack{i=0\\i \ne j}}^{N-1} t - t_i.\end{equation*} Now, for the first requirement,we can divide each term to yield $1$ when evaluated at $t = 3$(luckily, this will not violate the second requirement).If each term is $1$, the whole product will also be $1$:<jupyter_code>def polynomial_3(t):
    t = np.asarray(t)
    return (
        (t - (-1.5)) / (3 - (-1.5)) *
        (t - 0.5) / (3 - 0.5) *
        (t - 1.7) / (3 - 1.7) *
        (t - 4) / (3 - 4))
polynomial_3(ts)<jupyter_output><empty_output><jupyter_text>That's it! Generally,Lagrange basis polynomials can be written as\begin{equation*}\ell_j(t) = \prod_{\substack{i=0\\i \ne j}}^{N-1}\frac{t - t_i}{t_j - t_i}.\end{equation*} To get a better idea what's going on between the given time instances $t_i$,let's plot this polynomial(with a little help from [helper.py](helper.py)):<jupyter_code>from helper import grid_lines
plot_times = np.linspace(ts[0], ts[-1], 100)
plt.plot(plot_times, polynomial_3(plot_times))
grid_lines(ts, [0, 1])<jupyter_output><empty_output><jupyter_text>We can see from its shape that this is a polynomial of degree 4,which makes sense because the product we are using has 4 termscontaining one $t$ each.We can also see that it has the value $0$ at each of the initially providedtime instances $t_i$, except for $t_3 = 3$, where it has the value $1$.The above calculation can be easily generalized to be able to getany one of the set of polynomials defined by an arbitrary list of time instants:<jupyter_code>def lagrange_polynomial(times, i, t):
    """i-th Lagrange polynomial for the given time values, evaluated at t."""
    t = np.asarray(t)
    product = np.multiply.reduce
    return product([
        (t - times[j]) / (times[i] - times[j])
        for j in range(len(times))
        if i != j
    ])<jupyter_output><empty_output><jupyter_text>Now we can calculate and visualize all 5 polynomialsfor our 5 given time instants:<jupyter_code>polys = np.column_stack([lagrange_polynomial(ts, i, plot_times)
                         for i in range(len(ts))])
plt.plot(plot_times, polys)
grid_lines(ts, [0, 1])<jupyter_output><empty_output><jupyter_text>Finally, the interpolated values can be obtainedby applying the given $x_i$ values as weights to the polynomialsand summing everything together:<jupyter_code>weighted_polys = polys * xs
interpolated = np.sum(weighted_polys, axis=-1)
plt.plot(plot_times, weighted_polys)
plt.plot(plot_times, interpolated, color='black', linestyle='dashed')
plt.scatter(ts, xs, color='black')
grid_lines(ts)<jupyter_output><empty_output><jupyter_text>Neville's AlgorithmAn alternative way to calculate interpolated values is[Neville's algorithm](https://en.wikipedia.org/wiki/Neville%27s_algorithm)(see alsoBarry and Goldman (1988),figure 2).We mention this algorithm mainly because it is referenced in the[derivation of non-uniform Catmull--Rom splines](catmull-rom-non-uniform.ipynb)and the[description of the Barry--Goldman algorithm](catmull-rom-barry-goldman.ipynb). As main building block, we need a linear interpolationbetween two values in a given time interval:<jupyter_code>def lerp(xs, ts, t):
    """Linear intERPolation.

    Returns the interpolated value(s) at time(s) *t*,
    given two values/vertices *xs* at times *ts*.

    The two x-values can be scalars or vectors,
    or even higher-dimensional arrays
    (as long as the shape of *t* is compatible).

    """
    x_begin, x_end = map(np.asarray, xs)
    t_begin, t_end = ts
    if not np.isscalar(t):
        # This allows using an array of *t* values:
        t = np.expand_dims(t, axis=-1)
    return (x_begin * (t_end - t) + x_end * (t - t_begin)) / (t_end - t_begin)<jupyter_output><empty_output><jupyter_text>In each stage of the algorithm,linear interpolation is used to interpolate between adjacent values,leading to one fewer value than in the stage before.The new values are used as input to the next stage and so on.When there is only one value left, this value is the result.The only tricky part isto choose the appropriate time interval for each interpolation.In the first stage,the intervals between the given time values are used.In the second stage,each time interval is combined with the following one,leading to one fewer time intervals in total.In the third stage,each time interval is combined with the following two intervals,and so on until the last stage,where all time intervals are combined into a single large interval.<jupyter_code>def neville(xs, ts, t):
    """Lagrange interpolation using Neville's algorithm.
    
    Returns the interpolated value(s) at time(s) *t*,
    given the values *xs* at times *ts*.
    
    """
    assert len(xs) == len(ts)
    while len(xs) > 1:
        step = len(ts) - len(xs) + 1
        xs = [
            lerp(*args, t)
            for args in zip(zip(xs, xs[1:]), zip(ts, ts[step:]))]
    return xs[0]
plt.plot(plot_times, neville(xs, ts, plot_times))
plt.scatter(ts, xs)
grid_lines(ts)<jupyter_output><empty_output><jupyter_text>Two-dimensional ExampleLagrange interpolation can of course also be used in higher-dimensional spaces.To show this, let's create a little class:<jupyter_code>class Lagrange:
    
    def __init__(self, vertices, grid):
        assert len(vertices) == len(grid)
        self.vertices = vertices
        self.grid = grid
    
    def evaluate(self, t):
        return neville(self.vertices, self.grid, t)<jupyter_output><empty_output><jupyter_text>Since this class has the same interface as the splinesthat will be discussed in the following sections,we can use a spline helper function from [helper.py](helper.py)for plotting:<jupyter_code>from helper import plot_spline_2d<jupyter_output><empty_output><jupyter_text>This time, we have a list of two-dimensional vectorsand the same list of associated times as before:<jupyter_code>l1 = Lagrange([(2, -2), (-1, 0), (1.3, 0.5), (3.14, 0), (1, -1)], ts)
plot_spline_2d(l1)<jupyter_output><empty_output><jupyter_text>Runge's PhenomenonThis seems to work reasonably well,but as indicated above,Lagrange implementation has a severe limitation.This limitation gets more apparent when using more vertices,which leads to a higher-degree polynomial.<jupyter_code>vertices = [
    (1, 0),
    (1, 2),
    (3, 0),
    (2, -1),
    (2.5, 1.5),
    (5, 2),
    (6, 1),
    (5, 0),
    (6, -1),
    (7, 2),
    (4, 3),
]
times = range(len(vertices))
l2 = Lagrange(vertices, times)
plot_spline_2d(l2)<jupyter_output><empty_output><jupyter_text>Here we see a severe overshooting effect,most pronounced at the beginning and the end of the curve.This effect is called[Runge's phenomenon](https://en.wikipedia.org/wiki/Runge's_phenomenon).Long story short,Lagrange interpolation is typically not usable for drawing curves.For comparison, let's use the same positions and time valuesand create a [Catmull--Rom spline](catmull-rom.ipynb):<jupyter_code>import splines
cr_spline = splines.CatmullRom(vertices, times)
plot_spline_2d(cr_spline)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10659517426273458, 'fraction_numerical': 0.015013404825737266, 'mean_word_length': 4.561121049493143, 'pattern_counts': {'":': 0, '<': 66, '<?xml version=': 0, '>': 66, 'https://': 5, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '16', 'dir': 'jupyter-structured-clean-dedup', 'id': '225380', 'license': 'MIT', 'n_tokens_mistral': 3105, 'n_tokens_neox': 2938, 'n_words': 1039, 'path': 'doc/euclidean/lagrange.ipynb', 'repo_name': 'AudioSceneDescriptionFormat/splines', 'size': '9325'}