From 508c1e7c57469d833bd109ab9af64d2f7c3ee668 Mon Sep 17 00:00:00 2001
From: Wolfgang Kastaun <wolfgang.kastaun@aei.mpg.de>
Date: Fri, 20 Dec 2024 22:31:33 +0100
Subject: [PATCH] Consistently specify lagrange interpolators by order (not
 stencil size)

---
 lisainstrument/dynamic_delay_dask.py    |  6 ++---
 lisainstrument/dynamic_delay_numpy.py   |  6 ++---
 lisainstrument/fixed_shift_dask.py      | 10 +++----
 lisainstrument/fixed_shift_dsp.py       | 36 ++++++++++++-------------
 lisainstrument/fixed_shift_numpy.py     | 25 ++++++++---------
 lisainstrument/regular_interpolators.py | 24 ++++++++---------
 tests/test_dynamic_delay_dask.py        |  5 ++--
 tests/test_dynamic_delay_dsp.py         | 12 ++++-----
 tests/test_fixed_shift_dask.py          | 20 ++++++--------
 tests/test_fixed_shift_numpy.py         |  9 +++----
 10 files changed, 71 insertions(+), 82 deletions(-)

diff --git a/lisainstrument/dynamic_delay_dask.py b/lisainstrument/dynamic_delay_dask.py
index 4c52e8a..aa7d7a0 100644
--- a/lisainstrument/dynamic_delay_dask.py
+++ b/lisainstrument/dynamic_delay_dask.py
@@ -170,7 +170,7 @@ def make_dynamic_shift_linear_dask(
 
 
 def make_dynamic_shift_lagrange_dask(
-    length: int,
+    order: int,
     min_delay: float,
     max_delay: float,
     left_bound: ShiftBC,
@@ -179,7 +179,7 @@ def make_dynamic_shift_lagrange_dask(
     """Set up DynamicShiftDask instance with Lagrange interpolation method.
 
     Arguments:
-        length: Number of lagrange polynomials (of order length - 1)
+        order: Order of the Lagrange polynomials
         min_delay: Assume that any shift > -max_delay
         max_delay: Assume that any shift < -min_delay
         left_bound: Treatment of left boundary
@@ -189,7 +189,7 @@ def make_dynamic_shift_lagrange_dask(
         Interpolation function of type DynamicShiftDask
     """
 
-    interp = make_regular_interpolator_lagrange(length)
+    interp = make_regular_interpolator_lagrange(order)
     cfg = DynShiftCfg(min_delay, max_delay, left_bound, right_bound)
     return DynamicShiftDask(cfg, interp)
 
diff --git a/lisainstrument/dynamic_delay_numpy.py b/lisainstrument/dynamic_delay_numpy.py
index 0c277c9..6dc71f9 100644
--- a/lisainstrument/dynamic_delay_numpy.py
+++ b/lisainstrument/dynamic_delay_numpy.py
@@ -177,7 +177,7 @@ class DynamicShiftNumpy:
 
 
 def make_dynamic_shift_lagrange_numpy(
-    length: int,
+    order: int,
     min_delay: float,
     max_delay: float,
     left_bound: ShiftBC,
@@ -186,7 +186,7 @@ def make_dynamic_shift_lagrange_numpy(
     """Set up DynamicShiftNumpy instance with Lagrange interpolation method.
 
     Arguments:
-        length: Number of lagrange polynomials (of order length - 1)
+        order: Order of the Lagrange polynomials
         min_delay: Assume that any shift > -max_delay
         max_delay: Assume that any shift < -min_delay
         left_bound: Treatment of left boundary
@@ -195,7 +195,7 @@ def make_dynamic_shift_lagrange_numpy(
     Returns:
         Interpolation function of type DynamicShiftNumpy
     """
-    interp = make_regular_interpolator_lagrange(length)
+    interp = make_regular_interpolator_lagrange(order)
     cfg = DynShiftCfg(min_delay, max_delay, left_bound, right_bound)
     return DynamicShiftNumpy(cfg, interp)
 
diff --git a/lisainstrument/fixed_shift_dask.py b/lisainstrument/fixed_shift_dask.py
index 4db20a3..5b26022 100644
--- a/lisainstrument/fixed_shift_dask.py
+++ b/lisainstrument/fixed_shift_dask.py
@@ -142,17 +142,17 @@ class FixedShiftDask:  # pylint: disable=too-few-public-methods
 
 
 def make_fixed_shift_lagrange_dask(
-    left_bound: ShiftBC, right_bound: ShiftBC, length: int
+    left_bound: ShiftBC, right_bound: ShiftBC, order: int
 ) -> FixedShiftDask:
     """Create a FixedShiftDask instance that uses Lagrange interpolator
 
     Arguments:
-        left_bound: boundary treatment on the left
-        right_bound: boundary treatment on the right
-        length: number of Lagrange plolynomials (=order + 1)
+        left_bound: Boundary treatment on the left
+        right_bound: Boundary treatment on the right
+        order: Order of Lagrange plolynomials
 
     Returns:
         Fixed shift interpolator
     """
-    fac = FixedShiftLagrange.factory(length)
+    fac = FixedShiftLagrange.factory(order)
     return FixedShiftDask(left_bound, right_bound, fac)
diff --git a/lisainstrument/fixed_shift_dsp.py b/lisainstrument/fixed_shift_dsp.py
index a85d643..131d3b4 100644
--- a/lisainstrument/fixed_shift_dsp.py
+++ b/lisainstrument/fixed_shift_dsp.py
@@ -36,27 +36,27 @@ class FixedShiftWrapDSP(FixedShiftCore):
     See FixedShiftCore for general properties not specific to the interpolation method.
     """
 
-    def __init__(self, length: int, shift: float):
+    def __init__(self, order: int, shift: float):
         """Set up interpolation parameters.
 
-        The length parameter specifies the number of interpolation polynomials,
-        which have order length - 1. The length is also the number of samples used
-        for each interpolation point. The order of the interpolating polynomials
-        is also the order of plynomials that are interpolated with zero error.
+        The order parameter specifies the order of the interpolation polynomials. The
+        number of samples used for each interpolation point is order + 1. The order of
+        the interpolating polynomials is also the order of plynomials that are interpolated
+        with zero error.
 
         The shift sign is defined such that positive values refer to locations
         to the right of the output sample at hand.
 
         Arguments:
-            length: Size of the interpolation stencil
+            order: The order of the interpolation polynomials
             shift: The constant shift
         """
 
-        if length <= 1:
-            msg = f"FixedShiftWrapDSP: stencil length must be 2 or more, got {length}"
+        if order < 1:
+            msg = f"FixedShiftWrapDSP: order must be >= 1, got {order}"
             raise ValueError(msg)
 
-        if length % 2 != 0:
+        if order % 2 == 0:
             msg = "FixedShiftWrapDSP: even Lagrange polynomial order not supported by dsp.timeshift"
             raise ValueError(msg)
 
@@ -64,13 +64,13 @@ class FixedShiftWrapDSP(FixedShiftCore):
         loc_int = int(np.floor(loc))
         # ~ loc_frac = loc - loc_int
 
-        margin_lagr_left = (length - 1) // 2
-        margin_lagr_right = length - 1 - margin_lagr_left
+        margin_lagr_left = order // 2
+        margin_lagr_right = order - margin_lagr_left
 
         self._margin_left: Final = max(0, margin_lagr_left - loc_int)
         self._margin_right: Final = max(0, margin_lagr_right + loc_int)
         self._shift: Final = shift
-        self._order: Final = length - 1
+        self._order: Final = order
 
     @property
     def margin_left(self) -> int:
@@ -107,11 +107,11 @@ class FixedShiftWrapDSP(FixedShiftCore):
         return make_numpy_array_1d(res)
 
     @staticmethod
-    def factory(length: int) -> FixedShiftFactory:
+    def factory(order: int) -> FixedShiftFactory:
         """Factory for making FixedShiftWrapDSP instances
 
         Arguments:
-            length: number of Lagrange polynomials to use
+            order: Order of the Lagrange polynomials to use
 
         Returns:
             Factory function for making FixedShiftWrapDSP instances from shift
@@ -122,23 +122,23 @@ class FixedShiftWrapDSP(FixedShiftCore):
             Arguments:
                 shift: The fixed shift
             """
-            return FixedShiftWrapDSP(length, shift)
+            return FixedShiftWrapDSP(order, shift)
 
         return factory
 
 
 def make_fixed_shift_dsp_dask(
-    left_bound: ShiftBC, right_bound: ShiftBC, length: int
+    left_bound: ShiftBC, right_bound: ShiftBC, order: int
 ) -> FixedShiftDask:
     """Create a FixedShiftDask instance that uses dsp.timeshift interpolator
 
     Arguments:
         left_bound: boundary treatment on the left
         right_bound: boundary treatment on the right
-        length: number of Lagrange polynomials (=order + 1)
+        order: Order of Lagrange polynomials
 
     Returns:
         Fixed shift interpolator
     """
-    fac = FixedShiftWrapDSP.factory(length)
+    fac = FixedShiftWrapDSP.factory(order)
     return FixedShiftDask(left_bound, right_bound, fac)
diff --git a/lisainstrument/fixed_shift_numpy.py b/lisainstrument/fixed_shift_numpy.py
index 5aec405..53e4756 100644
--- a/lisainstrument/fixed_shift_numpy.py
+++ b/lisainstrument/fixed_shift_numpy.py
@@ -157,22 +157,23 @@ class FixedShiftLagrange(FixedShiftCore):
     See FixedShiftCore for general properties not specific to the interpolation method.
     """
 
-    def __init__(self, length: int, shift: float):
+    def __init__(self, order: int, shift: float):
         """Set up interpolation parameters.
 
-        The length parameter specifies the number of interpolation polynomials,
-        which have order length - 1. The length is also the number of samples used
-        for each interpolation point. The order of the interpolating polynomials
-        is also the order of plynomials that are interpolated with zero error.
+        The order parameter specifies the order of the interpolation polynomials. The
+        number of samples used for each interpolation point is order + 1. The order of
+        the interpolating polynomials is also the order of plynomials that are interpolated
+        with zero error.
 
         The shift sign is defined such that positive values refer to locations
         to the right of the output sample at hand.
 
         Arguments:
-            length: Size of the interpolation stencil
+            order: Order of the interpolation polynomials
             shift: The constant shift
         """
 
+        length = order + 1
         loc = float(shift)
         if length % 2 == 0:
             loc_int = int(np.floor(loc))
@@ -217,11 +218,11 @@ class FixedShiftLagrange(FixedShiftCore):
         return make_numpy_array_1d(self._filt(a))
 
     @staticmethod
-    def factory(length: int) -> FixedShiftFactory:
+    def factory(order: int) -> FixedShiftFactory:
         """Factory for making FixedShiftLagrange instances
 
         Arguments:
-            length: number of Lagrange polynomials to use
+            order: The order of the Lagrange polynomials to use
 
         Returns:
             Factory function for making FixedShiftLagrange instances from shift
@@ -236,7 +237,7 @@ class FixedShiftLagrange(FixedShiftCore):
             Returns:
                 FixedShiftLagrange instance
             """
-            return FixedShiftLagrange(length, shift)
+            return FixedShiftLagrange(order, shift)
 
         return factory
 
@@ -356,17 +357,17 @@ class FixedShiftNumpy:  # pylint: disable=too-few-public-methods
 
 
 def make_fixed_shift_lagrange_numpy(
-    left_bound: ShiftBC, right_bound: ShiftBC, length: int
+    left_bound: ShiftBC, right_bound: ShiftBC, order: int
 ) -> FixedShiftNumpy:
     """Create a FixedShiftNumpy instance that uses Lagrange interpolator
 
     Arguments:
         left_bound: boundary treatment on the left
         right_bound: boundary treatment on the right
-        length: number of Lagrange plolynomials (=order + 1)
+        order: Order of the Lagrange plolynomials
 
     Returns:
         Fixed shift interpolator
     """
-    fac = FixedShiftLagrange.factory(length)
+    fac = FixedShiftLagrange.factory(order)
     return FixedShiftNumpy(left_bound, right_bound, fac)
diff --git a/lisainstrument/regular_interpolators.py b/lisainstrument/regular_interpolators.py
index cc54d10..e5691f5 100644
--- a/lisainstrument/regular_interpolators.py
+++ b/lisainstrument/regular_interpolators.py
@@ -166,22 +166,20 @@ class RegularInterpLagrange(RegularInterpCore):
             filts.append(filt)
         return filts
 
-    def __init__(self, length: int):
+    def __init__(self, order: int):
         """Set up interpolation parameters.
 
-        The length parameter specifies the number of interpolation polynomials,
-        which have order length - 1. The length is also the number of samples used
-        for each interpolation point. The order of the interpolating polynomials
-        is also the order of plynomials that are interpolated with zero error.
+        The order parameter specifies the order of the interpolation polynomials. The
+        number of samples used for each interpolation point is order + 1. The order of
+        the interpolating polynomials is also the order of plynomials that are interpolated
+        with zero error.
 
         Arguments:
-            length: Size of the interpolation stencil
+            order: order of the interpolation polynomials
         """
-
+        length = order + 1
         if length <= 1:
-            msg = (
-                f"RegularInterpLagrange: stencil length must be 2 or more, got {length}"
-            )
+            msg = f"RegularInterpLagrange: order must be >= 1, got {length}"
             raise ValueError(msg)
 
         offset = -((length - 1) // 2)
@@ -478,17 +476,17 @@ class RegularInterpolator:
         return make_numpy_array_1d(res)
 
 
-def make_regular_interpolator_lagrange(length: int) -> RegularInterpolator:
+def make_regular_interpolator_lagrange(order: int) -> RegularInterpolator:
     """Create an interpolator using Lagrange interpolation
 
     See RegularInterpLagrange for details of the method.
 
     Arguments:
-        length: size of interpolation stencil
+        order: order of the interpolating polynomials
     Returns:
         Interpolation function
     """
-    return RegularInterpolator(RegularInterpLagrange(length))
+    return RegularInterpolator(RegularInterpLagrange(order))
 
 
 def make_regular_interpolator_linear() -> RegularInterpolator:
diff --git a/tests/test_dynamic_delay_dask.py b/tests/test_dynamic_delay_dask.py
index ce444d3..3b7e543 100644
--- a/tests/test_dynamic_delay_dask.py
+++ b/tests/test_dynamic_delay_dask.py
@@ -38,7 +38,6 @@ def test_dynamic_shift_linear_dask() -> None:
 
 def test_dynamic_shift_lagrange_dask() -> None:
     order = 5
-    length = order + 1
 
     t, dt = np.linspace(-5.345, 10.345, 1003, retstep=True)
 
@@ -57,7 +56,7 @@ def test_dynamic_shift_lagrange_dask() -> None:
     d = (0.93456456 + 0.0235345 * np.cos(4.3354 * t)) / dt
 
     op_da = make_dynamic_shift_lagrange_dask(
-        length, d.min(), d.max(), ShiftBC.FLAT, ShiftBC.EXCEPTION
+        order, d.min(), d.max(), ShiftBC.FLAT, ShiftBC.EXCEPTION
     )
     op_na = numpyfy_dask_multi(op_da, chunks=19)
 
@@ -69,7 +68,7 @@ def test_dynamic_shift_lagrange_dask() -> None:
     )
 
     op_np = make_dynamic_shift_lagrange_numpy(
-        length, d.min(), d.max(), ShiftBC.FLAT, ShiftBC.EXCEPTION
+        order, d.min(), d.max(), ShiftBC.FLAT, ShiftBC.EXCEPTION
     )
 
     s_np = op_np(y, -d)
diff --git a/tests/test_dynamic_delay_dsp.py b/tests/test_dynamic_delay_dsp.py
index ad24b06..59da8fc 100644
--- a/tests/test_dynamic_delay_dsp.py
+++ b/tests/test_dynamic_delay_dsp.py
@@ -71,7 +71,6 @@ def test_dynamic_shift_dsp_dask_orders() -> None:
     Compare against other lagrange interpolator"""
 
     for order in (5, 7, 13, 31):
-        length = order + 1
         t, dt = np.linspace(-5.345, 10.345, 1003, retstep=True)
         y = g(t)
 
@@ -85,7 +84,7 @@ def test_dynamic_shift_dsp_dask_orders() -> None:
         s1_na = op1_na(y, -d)
 
         op2_da = make_dynamic_shift_lagrange_dask(
-            length, d.min(), d.max(), ShiftBC.FLAT, ShiftBC.EXCEPTION
+            order, d.min(), d.max(), ShiftBC.FLAT, ShiftBC.EXCEPTION
         )
         op2_na = numpyfy_dask_multi(op2_da, chunks=19)
 
@@ -103,17 +102,16 @@ def test_fixed_shift_dsp_dask() -> None:
     they should
     """
     order = 5
-    length = order + 1
 
     t, dt = np.linspace(-5.345, 10.345, 1003, retstep=True)
 
     y = g(t)
 
     for d in (0.93456456 / dt, 3.1, 3.5, 3.9, 4.0, 4.1, 4.5, 4.9, 116.0):
-        op1_da = make_fixed_shift_lagrange_dask(ShiftBC.FLAT, ShiftBC.EXCEPTION, length)
+        op1_da = make_fixed_shift_lagrange_dask(ShiftBC.FLAT, ShiftBC.EXCEPTION, order)
         op1_na = numpyfy_dask_multi(op1_da, chunks=113)
 
-        op2_da = make_fixed_shift_dsp_dask(ShiftBC.FLAT, ShiftBC.EXCEPTION, length)
+        op2_da = make_fixed_shift_dsp_dask(ShiftBC.FLAT, ShiftBC.EXCEPTION, order)
         op2_na = numpyfy_dask_multi(op2_da, chunks=113)
 
         s1_na = op1_na(y, -d)
@@ -124,10 +122,10 @@ def test_fixed_shift_dsp_dask() -> None:
         with pytest.raises(RuntimeError):
             op2_na(y, d)
 
-        op3_da = make_fixed_shift_lagrange_dask(ShiftBC.EXCEPTION, ShiftBC.FLAT, length)
+        op3_da = make_fixed_shift_lagrange_dask(ShiftBC.EXCEPTION, ShiftBC.FLAT, order)
         op3_na = numpyfy_dask_multi(op3_da, chunks=113)
 
-        op4_da = make_fixed_shift_dsp_dask(ShiftBC.EXCEPTION, ShiftBC.FLAT, length)
+        op4_da = make_fixed_shift_dsp_dask(ShiftBC.EXCEPTION, ShiftBC.FLAT, order)
         op4_na = numpyfy_dask_multi(op4_da, chunks=113)
 
         s3_na = op3_na(y, d)
diff --git a/tests/test_fixed_shift_dask.py b/tests/test_fixed_shift_dask.py
index c74037e..f8c20b0 100644
--- a/tests/test_fixed_shift_dask.py
+++ b/tests/test_fixed_shift_dask.py
@@ -33,15 +33,14 @@ def test_fixed_shift_lagrange_dask() -> None:
     expected result (excluding points affected by boundary treatment).
     """
     order = 5
-    length = order + 1
 
     t, dt = np.linspace(-5.345, 10.345, 1003, retstep=True, dtype=float)
     y = g(t)
 
     for d in (0.93456456 / dt, 3.1, 3.5, 3.9, 4.0, 4.1, 4.5, 4.9, 116.0):
 
-        op_np = make_fixed_shift_lagrange_numpy(ShiftBC.FLAT, ShiftBC.EXCEPTION, length)
-        op_da = make_fixed_shift_lagrange_dask(ShiftBC.FLAT, ShiftBC.EXCEPTION, length)
+        op_np = make_fixed_shift_lagrange_numpy(ShiftBC.FLAT, ShiftBC.EXCEPTION, order)
+        op_da = make_fixed_shift_lagrange_dask(ShiftBC.FLAT, ShiftBC.EXCEPTION, order)
         op_na = numpyfy_dask_multi(op_da, chunks=113)
 
         s_np = op_np(y, -d)
@@ -56,17 +55,15 @@ def test_fixed_shift_lagrange_dask() -> None:
         with pytest.raises(RuntimeError):
             op_na(y, d)
 
-        op2_np = make_fixed_shift_lagrange_numpy(
-            ShiftBC.EXCEPTION, ShiftBC.FLAT, length
-        )
-        op2_da = make_fixed_shift_lagrange_dask(ShiftBC.EXCEPTION, ShiftBC.FLAT, length)
+        op2_np = make_fixed_shift_lagrange_numpy(ShiftBC.EXCEPTION, ShiftBC.FLAT, order)
+        op2_da = make_fixed_shift_lagrange_dask(ShiftBC.EXCEPTION, ShiftBC.FLAT, order)
         op2_na = numpyfy_dask_multi(op2_da, chunks=113)
 
         s2_np = op2_np(y, d)
         s2_na = op2_na(y, d)
         s2_ex = g(t + d * dt)
 
-        margin2_ex = int(np.floor(d)) + (length - 1 - order // 2)
+        margin2_ex = int(np.floor(d)) + (order - order // 2)
 
         assert s2_na == pytest.approx(s2_np, rel=1e-14, abs=0)
         assert s2_ex[:-margin2_ex] == pytest.approx(
@@ -76,7 +73,7 @@ def test_fixed_shift_lagrange_dask() -> None:
         with pytest.raises(RuntimeError):
             op2_na(y, -d)
 
-        op3_da = make_fixed_shift_lagrange_dask(ShiftBC.FLAT, ShiftBC.FLAT, length)
+        op3_da = make_fixed_shift_lagrange_dask(ShiftBC.FLAT, ShiftBC.FLAT, order)
         op3_na = numpyfy_dask_multi(op3_da, chunks=113)
 
         s3_na = op3_na(y, -d)
@@ -96,18 +93,17 @@ def test_shift_lagrange_dask_fixed_vs_dyn() -> None:
     y = g(t)
 
     for order in (5, 6, 31):
-        length = order + 1
 
         for d in (0.93456456 / dt, 23.1, 23.5, 23.9, 24.0, 116.0):
             d1d = np.ones_like(y) * d
 
             op1_da = make_dynamic_shift_lagrange_dask(
-                length, d, d, ShiftBC.FLAT, ShiftBC.EXCEPTION
+                order, d, d, ShiftBC.FLAT, ShiftBC.EXCEPTION
             )
             op1_na = numpyfy_dask_multi(op1_da, chunks=113)
 
             op2_da = make_fixed_shift_lagrange_dask(
-                ShiftBC.FLAT, ShiftBC.EXCEPTION, length
+                ShiftBC.FLAT, ShiftBC.EXCEPTION, order
             )
             op2_na = numpyfy_dask_multi(op2_da, chunks=113)
 
diff --git a/tests/test_fixed_shift_numpy.py b/tests/test_fixed_shift_numpy.py
index 26277d0..907cbd2 100644
--- a/tests/test_fixed_shift_numpy.py
+++ b/tests/test_fixed_shift_numpy.py
@@ -9,7 +9,6 @@ from lisainstrument.fixed_shift_numpy import make_fixed_shift_lagrange_numpy
 
 def test_fixed_shift_lagrange_dask() -> None:
     order = 5
-    length = order + 1
 
     t, dt = np.linspace(-5.345, 10.345, 1003, retstep=True)
 
@@ -27,7 +26,7 @@ def test_fixed_shift_lagrange_dask() -> None:
 
     for d in (0.93456456 / dt, 3.1, 3.5, 3.9, 4.0, 4.1, 4.5, 4.9):
 
-        op_np = make_fixed_shift_lagrange_numpy(ShiftBC.FLAT, ShiftBC.EXCEPTION, length)
+        op_np = make_fixed_shift_lagrange_numpy(ShiftBC.FLAT, ShiftBC.EXCEPTION, order)
 
         s_np = op_np(y, -d)
         s_ex = g(t - d * dt)
@@ -41,14 +40,12 @@ def test_fixed_shift_lagrange_dask() -> None:
         with pytest.raises(RuntimeError):
             op_np(y, d)
 
-        op2_np = make_fixed_shift_lagrange_numpy(
-            ShiftBC.EXCEPTION, ShiftBC.FLAT, length
-        )
+        op2_np = make_fixed_shift_lagrange_numpy(ShiftBC.EXCEPTION, ShiftBC.FLAT, order)
 
         s2_np = op2_np(y, d)
         s2_ex = g(t + d * dt)
 
-        margin2_ex = int(np.floor(d)) + (length - 1 - order // 2)
+        margin2_ex = int(np.floor(d)) + (order - order // 2)
 
         assert s2_ex[:-margin2_ex] == pytest.approx(
             s2_np[:-margin2_ex], abs=1e-15, rel=5e-13
-- 
GitLab