Merge pull request #743 from MothNik/dev

Minor improvements for `signalprocessing.Interp` and `signalprocessing.InterpCubicSpline`

Author: mrava87

pylops/signalprocessing/_interp_utils.py

@@ -28,7 +28,7 @@ def _ensure_iava_is_unique(
 def _clip_iava_above_last_sample_index(
     iava: NDArray,
     sample_size: int,
-) -> None:
+) -> NDArray:
     """
     Ensures that elements in ``iava`` do not exceed the last sample index.
     Elements above the penultimate sample are clipped to the next closest float value
@@ -47,8 +47,9 @@ def _clip_iava_above_last_sample_index(
         # NOTE: ``numpy.nextafter(x, -np.inf)`` gives the closest float-value that is
         #       less than ``x``, i.e., this logic clips ``iava`` to the highest possible
         #       value that is still below the last sample
+        iava = iava.copy()  # to avoid silent input mutation
         iava[np.where(outside)] = np.nextafter(last_sample_index, -np.inf)
 
     _ensure_iava_is_unique(iava=iava)
 
-    return
+    return iava

pylops/signalprocessing/interp.py

@@ -47,7 +47,10 @@ def _linearinterp(
 
     # ensure that samples are not beyond the last sample, in that case set to
     # penultimate sample and raise a warning
-    _clip_iava_above_last_sample_index(iava=iava, sample_size=sample_size)
+    iava = _clip_iava_above_last_sample_index(  # type: ignore
+        iava=iava,  # type: ignore
+        sample_size=sample_size,
+    )
 
     # find indices and weights
     iva_l = ncp.floor(iava).astype(int)
@@ -133,7 +136,7 @@ def Interp(
       polynomial fitted between ``np.floor(iava)`` and ``np.floor(iava) + 1``.
       It offers an excellent tradeoff between accuracy and computational complexity
       and its results oscillate less than those obtained from sinc interpolation.
-      It can also be accessed directly via :class:`pylops.singalprocessing.InterpCubicSpline`.
+      It can also be accessed directly via :class:`pylops.signalprocessing.InterpCubicSpline`.
 
     .. note:: The vector ``iava`` should contain unique values. If the same
       index is repeated twice an error will be raised. This also applies when
@@ -179,7 +182,7 @@ def Interp(
     ValueError
         If the vector ``iava`` contains repeated values.
     NotImplementedError
-        If ``kind`` is not ``nearest``, ``linear`` or ``sinc``
+        If ``kind`` is not ``nearest``, ``linear``, ``sinc``, or ``cubic_spline``
 
     See Also
     --------

pylops/signalprocessing/interpspline.py

@@ -210,7 +210,7 @@ class _BandedLUDecomposition:
     Represents the LU decomposition of a general banded matrix as performed by the
     LAPACK routines ``?gbtrf``.
     This class was implemented for spline interpolations between only 2 data points
-    because the class :class:`_BandedLUDecomposition` uses the LAPACK routines
+    because the class :class:`_TridiagonalLUDecomposition` uses the LAPACK routines
     ``?gttrf`` that cannot handle 2 x 2 tridiagonal matrices.
 
     """
@@ -251,7 +251,11 @@ def from_tridiagonal_matrix(
         banded_representation[2, ::] = matrix.main_diagonal
         banded_representation[3, 0:-1] = matrix.sub_diagonal
 
-        (lu_banded, pivot_indices, info,) = lapack_factorizer(
+        (
+            lu_banded,
+            pivot_indices,
+            info,
+        ) = lapack_factorizer(
             ab=banded_representation,
             kl=1,
             ku=1,
@@ -389,7 +393,6 @@ def from_tridiagonal_matrix(
 
         raise np.linalg.LinAlgError(
             f"Could not LU-factorize tridiagonal matrix! Got {info=}."
-            f"Could not LU-factorize tridiagonal matrix! Got {info=}."
         )
 
     def solve(
@@ -818,7 +821,10 @@ def __init__(
             )
 
         iava = np.asarray(iava, dtype=np.float64)
-        _clip_iava_above_last_sample_index(iava=iava, sample_size=num_cols)
+        iava = _clip_iava_above_last_sample_index(  # type: ignore
+            iava=iava,
+            sample_size=num_cols,
+        )
 
         if isinstance(bc_type, str) and bc_type.lower() in {"natural"}:
             self.bc_type = bc_type.lower()
@@ -950,7 +956,7 @@ def _rmatvec(self, x: InexactNDArray) -> InexactNDArray:
             x_mod[0 : self.num_cols]
             + self._rmatmat_difference_method(
                 self._lhs_B_matrix_transposed_lu.solve(
-                    rhs=x_mod[self.num_cols : x_mod.size],
+                    rhs=x_mod[self.num_cols :],
                     lapack_solver=self._tridiag_lu_solve,
                 )
             )

pytests/test_interpolation_spline.py

@@ -1,10 +1,12 @@
+from math import prod
 from typing import Final, Tuple
 
 import numpy as np
 import pytest
 from scipy.interpolate import CubicSpline
 
 from pylops.signalprocessing import InterpCubicSpline
+from pylops.utils import dottest
 
 TEST_ARRAY_SHAPE: Final[Tuple] = (
     20,
@@ -16,6 +18,90 @@
 MIN_NUM_TEST_SAMPLES: Final[int] = 1
 
 
+def test_cubic_spline_raises_on_not_supported_bc_type() -> None:
+    """
+    Tests whether ``pylops.signalprocessing.InterpCubicSpline`` raises a
+    ``NotImplementedError`` for boundary conditions that are not supported.
+
+    """
+
+    with pytest.raises(NotImplementedError):
+        InterpCubicSpline(
+            dims=(5, 2),
+            iava=np.array([0.5, 2.3]),
+            bc_type="erroneous",  # type: ignore
+        )
+
+
+@pytest.mark.parametrize(
+    "with_complex",
+    [
+        pytest.param(False, id="real"),
+        pytest.param(True, id="complex"),
+    ],
+)
+@pytest.mark.parametrize(
+    "axis",
+    [
+        0,
+        1,
+        2,
+        3,
+        -1,
+        -2,
+        -3,
+    ],
+)
+@pytest.mark.parametrize(
+    "subsample_fraction",
+    [
+        pytest.param(0.5, id="decimation"),
+        pytest.param(5.0, id="upsampling"),
+    ],
+)
+def test_natural_cubic_spline_dottest(
+    subsample_fraction: float,
+    axis: int,
+    with_complex: bool,
+) -> None:
+    """
+    Tests ``pylops.signalprocessing.InterpCubicSpline`` with the ``dottest``.
+
+    """
+
+    # Setup
+
+    num_samples = TEST_ARRAY_SHAPE[axis]
+    x_eval_fractions = np.random.rand(
+        max(
+            round(num_samples * subsample_fraction),
+            MIN_NUM_TEST_SAMPLES,
+        )
+    )
+    x_eval_for_pylops = (num_samples - 1) * x_eval_fractions
+
+    shape_list = list(TEST_ARRAY_SHAPE)
+    shape_list[axis] = x_eval_fractions.size  # type: ignore
+    num_rows = prod(shape_list, start=1)
+    num_columns = prod(TEST_ARRAY_SHAPE, start=1)
+
+    # Test
+
+    splinop = InterpCubicSpline(
+        dims=TEST_ARRAY_SHAPE,
+        iava=x_eval_for_pylops,
+        axis=axis,
+        dtype="complex128" if with_complex else "float64",
+    )
+
+    assert dottest(
+        Op=splinop,
+        nr=num_rows,
+        nc=num_columns,
+        complexflag=0 if not with_complex else 3,
+    )
+
+
 @pytest.mark.parametrize(
     "with_complex",
     [