Cubic Spline Interpolator

[MothNik]

This pull request aims to close #632 by adding the kind="cubic_spline"-option to pylops.signalprocessing.Interp.
It is implemented in a way that can handle 1D-Arrays and 2D-Arrays along a single axis (0 or 1) in a vectorized fashion that heavily reuses pre-computations (optimised to match speed of kind="linear" while providing similar accuracy as kind="sinc" with less oscillations).

Currently, only a numpy-backend is available.

[MothNik]

Pull Request Checklist

• implement CubicSplineInterpolator InterpCubicSpline basic logic
• handle real and complex data types
• optimize runtime to be comparable with Interp(..., kind="linear")
• handle 1D-vectors
• handle 2D-Arrays
• handle ND-Arrays with axis-argument (thanks to @mrava87 🚀 )
• protect against wrong document availability of backends (tridiagonal solver seems unavailable for cupy)
• check if 2D-Array memory layout (C/F-order) matters
• integrate CubicSplineInterpolator InterpCubicSpline into existing tests
• validate CubicSplineInterpolator InterpCubicSpline against scipy.interpolate.CubicSpline(..., bc_type="natural")
• add tiny documentation section to Interp
• add mathematical documentation Notes section to CubicSplineInterpolator
• move CubicSplineInterpolator InterpCubicSpline to its own linear operator file
• fix floating point rounding error in iava-clipping for very large samples
• fix typing of iava to SamplingLike
• integrate InterpCubicSpline into pylops documentation

¹⁾ What's not obvious to me right now is how dims and dimsd work and why the @ (__matmul__) operator decides to flatten 2D-input before passing it down to the _matmat and _rmatmat methods.

[mrava87]

@MothNik thanks a lot for taking the time to turn the long stand Github issue into a PR ❤️

I am going to try to reply to all the points you have raised here and in your last comment in the Github issue, I will later do a proper review of the code and provide inline comments:

• handling of 2D/ND arrays: I recall we discussed this briefly in the Github issue. Basically in PyLops we historically allowed two type of inputs to an operator (I discuss for simplicity the case of forward, the same applies to adjoint)

  1. 1d arrays with size equal to the number of columns of the operator. Applying the operator is equivalent to a matvec if the operator was a dense matrix, and the output is a 1d array of size equal to the number of rows of the operator.
  2. 2d arrays (or matrices) where the number of rows equals to the number of columns of the operator and an arbitrary number of columns. Here the operator is applied to each column with a naive for...loop -

     pylops/pylops/linearoperator.py

     Line 451 in eeb8c47

     def _matmat(self, X: NDArray) -> NDArray:

  In PyLops v2 we also started allowing one to pass inputs that are ND-arrays in the natural shape an operator would work (so case 1 but with the input not being flattened)... this is the scenario we implement with the help of dims/dimsd as they are the unflattened version of shape[1]/shape[0] (basically np.prod(dims)=shape[1] and np.prod(dims)=shape[0]), and this is what I think you did not really implement - hence the tests failing. I have played a bit and found out that the current code wasn't that far for allowing this, see my commit for the changes that allow ANY ND-array (so not just 2d but also 3/4/ND) to work with any choice of axis.

• matmat: I see you implemented some custom _matmat/_rmatmat... as I explained above we do not do it as matmat isn't a very common use case for PyLops (it is there more for historical reasons as SciPy supports it)... I wonder if you actually did it thinking it will handle the 2D/ND array scenario I was asking?
  Also I tried this:

np.random.seed(1)
x = np.random.normal(0, 1, (10, 50))

perc_subsampling = 0.4
Nsub = int(np.round(10 * perc_subsampling))
iava = np.sort(np.random.permutation(np.arange(10))[:Nsub])

Iop, _ = Interp(
    10,
    iava,
    axis=0,
    kind="cubic_spline",
    dtype=par["dtype"],
)

Iop._matmat(x)

and seems to fail ValueError: all the input array dimensions except for the concatenation axis must match exactly, but along dimension 1, the array at index 0 has size 50 and the array at index 1 has size 52 🤔

• cupy: it is ok to have an operator that has only a numpy (aka cpu) backend, and your scenario (ie lack of equivalent in cupy to numpy/scipy method used in the operator) is a common one... for those operators we however do not feel it is important to have checks if one passes cupy arrays, but we simply state the capabilities of each operator in https://github.com/PyLops/pylops/blob/dev/docs/source/gpu.rst. In this case, since Interp will be partially CuPy/JAX compatible, we can use the warning emoji and had a warning at the bottom specifying that for spline the only backend is numpy... see Convolve1D for an example

• Notes: I think this is definitely important as this is not a trivial operator, and you had a nice description in the Github issue. And here comes a suggestion: since this operator is effectively implemented in CubicSplineInterpolator, I would make this a public class (calling the file interpspline and the operator InterpCubicSpline). We can still also have it as an option of Interp but having the class public we can more easily expose it in the doc and add a good documentation of what it does, I feel it would be too much to be part of the Notes of Interp (basically we can treat it the same way we treat nearest. I also suggest to add the spline interpolator to this example https://github.com/PyLops/pylops/blob/dev/examples/plot_restriction.py, and feel free to add more (like the example you had in the Github issue)

• Minor Finding No. 1 - Typing of iava: you are right, there is an inconsistency here, which I think originated from this massive PR (d765af9#diff-db0feeec6c30154190addf833c81fba3c7a670f4fce6579966fcc942307094a4) when we added a lot of new type hints. The correct type for iava should be NDArray as we expect floats as element of the array. In fact, if one passed just ints we will have already the correct values in the input and we could just fallback to the Restriction operator (and so the operator would work with int values, but it will be an overkill). So if you don't mind, could you fix the types to SamplingLike since in the doc we say we also allow lists 😄

• Minor Finding No. 2 - Clipping of iava - Interesting, I never noticed this (probably because I never had such large inputs), but worth fixing it now that you spotted a problem. I am just not sure why you have added the 1e6; this still seems to work and gives a number close to the actual last value:

>>> n * (1 - np.finfo(np.float64).eps)
np.float64(1999999.9999999995)

probably nothing that will be ever perceived, but just curious 😉

¹⁾ not sure what you mean, if I run the code I gave you above with Iop @ x instead of directly Iop._matmat(x) and I print the shape of x in your _matmat method, I get (10, 50), which is what I expect...

[MothNik]

@mrava87
Thanks a lot for taking the time for such detailed feedback, the explanations, and providing a fix!
I misunderstood the 2D-Array-handling then, and I must say that the handling of ND-Arrays is really a genius addition 🤩
I will go through your feedback and adapt the PR 👍

[mrava87]

Perfect 😄 Do not worry about some of the CI failing, these are problems unrelated to your PR (as long as the doc and main Github action CI passes it is all good...)

[MothNik]

@mrava87
I was adapting the tests in test_interpolation and I uncovered and fixed a lingering test bug (commit) in the tests for Bilinear (the tests, not the operator itself).

The lines

    with pytest.raises(ValueError):
        iava_rep = iava.copy()
        iava_rep[-2] = [0, 0]
        iava_rep[-1] = [0, 0]
        _, _ = Bilinear(
            iava_rep, dims=(par.y_num, par.x_num, par.t_num), dtype=par.dtype
        )

passed the assertion, but never hit the iava-duplicate check in Bilinear.
Changing this using regex-matching with a relevant part of the error message like

    with pytest.raises(ValueError, match="repeated"):
        iava_rep = iava.copy()
        iava_rep[-2] = [0, 0]
        iava_rep[-1] = [0, 0]
        _, _ = Bilinear(
            iava_rep, dims=(par.y_num, par.x_num, par.t_num), dtype=par.dtype
        )

revealed the following

>       with pytest.raises(ValueError, match="repeated"):
E       AssertionError: Regex pattern did not match.
E         Expected regex: 'repeated'
E         Actual message: 'could not broadcast input array from shape (2,) into shape (10,)'

The error happened already in iava_rep[-2] = [0, 0] due to a shape mismatch, which also raises a ValueError.

Checking against common exceptions like TypeError or ValueError with pytest.raises alone is relatively risky because they are frequently raised in other contexts. Then, pytest.raises just swallows the context and nobody can infer what happened unless closely inspecting the coverage reports. I recommed to use the match=... with a meaningful regex pattern to avoid tests passing without hitting the targeted branch.

Altneratively, and probably more readable, one can use something more verbose like

    with pytest.raises(ValueError) as exception_info:
        iava_rep = iava.copy()
        iava_rep[-2] = [0, 0]
        iava_rep[-1] = [0, 0]
        _, _ = Bilinear(
            iava_rep, dims=(par.y_num, par.x_num, par.t_num), dtype=par.dtype
        )
        
        error_message = str(exception_info.value)
        assert "iava" in error_message
        assert "repeated" in error message

[mrava87]

@mrava87

I was adapting the tests in test_interpolation and I uncovered and fixed a lingering test bug (commit) in the tests for Bilinear (the tests, not the operator itself).

The lines

    with pytest.raises(ValueError):

        iava_rep = iava.copy()

        iava_rep[-2] = [0, 0]

        iava_rep[-1] = [0, 0]

        _, _ = Bilinear(

            iava_rep, dims=(par.y_num, par.x_num, par.t_num), dtype=par.dtype

        )

passed the assertion, but never hit the iava-duplicate check in Bilinear.

Changing this using regex-matching with a relevant part of the error message like

    with pytest.raises(ValueError, match="repeated"):

        iava_rep = iava.copy()

        iava_rep[-2] = [0, 0]

        iava_rep[-1] = [0, 0]

        _, _ = Bilinear(

            iava_rep, dims=(par.y_num, par.x_num, par.t_num), dtype=par.dtype

        )

revealed the following

>       with pytest.raises(ValueError, match="repeated"):

E       AssertionError: Regex pattern did not match.

E         Expected regex: 'repeated'

E         Actual message: 'could not broadcast input array from shape (2,) into shape (10,)'

The error happened already in iava_rep[-2] = [0, 0] due to a shape mismatch, which also raises a ValueError.

Checking against common exceptions like TypeError or ValueError with pytest.raises alone is relatively risky because they are frequently raised in other contexts. Then, pytest.raises just swallows the context and nobody can infer what happened unless closely inspecting the coverage reports. I recommed to use the match=... with a meaningful regex pattern to avoid tests passing without hitting the targeted branch.

Altneratively, and probably more readable, one can use something more verbose like

    with pytest.raises(ValueError) as exception_info:

        iava_rep = iava.copy()

        iava_rep[-2] = [0, 0]

        iava_rep[-1] = [0, 0]

        _, _ = Bilinear(

            iava_rep, dims=(par.y_num, par.x_num, par.t_num), dtype=par.dtype

        )

        

        error_message = str(exception_info.value)

        assert "iava" in error_message

        assert "repeated" in error message

Oh yeah, good catch! I agree, sometimes our with... tests have too generic assertments leading to unexpected behaviors like the one you spotted. I like the second one so we can check multiple words or even partial sentences in the error message 😉

[mrava87]

@MothNik great job!

I went through the PR once again and left a few minor comments. Once those are addressed I think this PR is finally ready to be merged 🥳

[mrava87]

@MothNik let me know if you have time to make these final changes or if you want me to go ahead and do them before I can merge 😄

[MothNik]

@mrava87 Thanks for the thorough review!
I will go over it when I'm back from vacation after this week.
Once I'm done, I will drop a note an convert this draft to a PR.

[MothNik]

@mrava87
I implemented your changes and also increased the version in the last commit to 2.7.0 (I guess?).
I made a changelog entry and converted this draft to an actual PR.

[MothNik]

There is definitely also some follow-up work to do. The current $LU$-decomposition logic for the spline interpolation is borrowed from this scipy issue, but I think this should at one point belong to scipy and not be part of pylops.

[mrava87]

Fantastic, I will do one final review but I think this is likely ready to be merged 😄 and yes, if the logic of the LU-decomposition were to move to scipy that would be even better.. for now, happy to host it as part of the spline interpolator