Repository: scipy/scipy
Pull request: BUG: Fix `Akima1DInterpolator` by returning linear interpolant for `y.shape[0] == 2`
<!-- 
Thanks for contributing a pull request! Please ensure that
your PR satisfies the checklist before submitting:
https://scipy.github.io/devdocs/dev/contributor/development_workflow.html#checklist-before-submitting-a-pr

Also, please name and describe your PR as you would write a
commit message:
https://scipy.github.io/devdocs/dev/contributor/development_workflow.html#writing-the-commit-message.
However, please only include an issue number in the description, not the title,
and please ensure that any code names containing underscores are enclosed in backticks.

Depending on your changes, you can skip CI operations and save time and energy: 
https://scipy.github.io/devdocs/dev/contributor/continuous_integration.html#skipping

Note that we are a team of volunteers; we appreciate your
patience during the review process.

Again, thanks for contributing!
-->

#### Reference issue
<!--Example: Closes gh-WXYZ.-->

Closes https://github.com/scipy/scipy/issues/22011

#### What does this implement/fix?
<!--Please explain your changes.-->

Does the fix as described in https://github.com/scipy/scipy/issues/22011#issuecomment-2525067847

cc: @ev-br 

#### Additional information
<!--Any additional information you think is important.-->


Problem statement:
BUG: interpolate.Akima1DInterpolator: different values on subsequent calls
### Describe your issue.

The first and third calls to Akima1DInterpolator pass the same parameters but the results are different. The result of the first and third call is the same if the second call is removed:




### Reproducing Code Example

```python
import numpy as np
import scipy

x = np.array([0.0, 1.0], dtype=float)
y = np.array([0.55, 0.55])
akima = scipy.interpolate.Akima1DInterpolator(x, y, axis=0, extrapolate=None)
print(akima(0.45))

x = np.array([0.0, 0.5, 1.0])
y = np.array([0.15, 0.67, 0.15])
akima = scipy.interpolate.Akima1DInterpolator(x, y, axis=0, extrapolate=None)
print(akima(0.45))

x = np.array([0.0, 1.0])
y = np.array([0.55,  0.55])
akima = scipy.interpolate.Akima1DInterpolator(x, y, axis=0, extrapolate=None)
print(akima(0.45))
```


### Error message

```shell
Output:

0.55
0.6648000000000001
0.2926
```


### SciPy/NumPy/Python version and system information

```shell
1.14.1 1.26.4 sys.version_info(major=3, minor=10, micro=7, releaselevel='final', serial=0)
Build Dependencies:
  blas:
    detection method: pkgconfig
    found: true
    include directory: C:/Users/runneradmin/AppData/Local/Temp/cibw-run-_mbvtoth/cp310-win_amd64/build/venv/Lib/site-packages/scipy_openblas32/include
    lib directory: C:/Users/runneradmin/AppData/Local/Temp/cibw-run-_mbvtoth/cp310-win_amd64/build/venv/Lib/site-packages/scipy_openblas32/lib
    name: scipy-openblas
    openblas configuration: OpenBLAS 0.3.27.dev DYNAMIC_ARCH NO_AFFINITY Zen MAX_THREADS=24
    pc file directory: D:/a/scipy/scipy
    version: 0.3.27.dev
  lapack:
    detection method: pkgconfig
    found: true
    include directory: C:/Users/runneradmin/AppData/Local/Temp/cibw-run-_mbvtoth/cp310-win_amd64/build/venv/Lib/site-packages/scipy_openblas32/include
    lib directory: C:/Users/runneradmin/AppData/Local/Temp/cibw-run-_mbvtoth/cp310-win_amd64/build/venv/Lib/site-packages/scipy_openblas32/lib
    name: scipy-openblas
    openblas configuration: OpenBLAS 0.3.27.dev DYNAMIC_ARCH NO_AFFINITY Zen MAX_THREADS=24
    pc file directory: D:/a/scipy/scipy
    version: 0.3.27.dev
  pybind11:
    detection method: config-tool
    include directory: unknown
    name: pybind11
    version: 2.12.0
Compilers:
  c:
    commands: cc
    linker: ld.bfd
    name: gcc
    version: 10.3.0
  c++:
    commands: c++
    linker: ld.bfd
    name: gcc
    version: 10.3.0
  cython:
    commands: cython
    linker: cython
    name: cython
    version: 3.0.11
  fortran:
    commands: gfortran
    linker: ld.bfd
    name: gcc
    version: 10.3.0
  pythran:
    include directory: C:\Users\runneradmin\AppData\Local\Temp\pip-build-env-w8i5lw_j\overlay\Lib\site-packages/pythran
    version: 0.16.1
Machine Information:
  build:
    cpu: x86_64
    endian: little
    family: x86_64
    system: windows
  cross-compiled: false
  host:
    cpu: x86_64
    endian: little
    family: x86_64
    system: windows
Python Information:
  path: C:\Users\runneradmin\AppData\Local\Temp\cibw-run-_mbvtoth\cp310-win_amd64\build\venv\Scripts\python.exe
  version: '3.10'
```

Hints from issue discussion:
Confirmed. Wow. The last result is nonsense. And it depends on previous calls, so it's some unitialized memory. And indeed, changing `np.empty` to `np.zeros` in https://github.com/scipy/scipy/blob/main/scipy/interpolate/_cubic.py#L511 fixes the immediate problem.
I'm not immediately convinced that np.zeros is the right fix though. 

Clearly, the Akima1D code is not expecting `len(x) == 2`. 
Note that `PchipInterpolator` bails out early in this case: https://github.com/scipy/scipy/blob/main/scipy/interpolate/_cubic.py#L282

So I believe this would be the right fix: detect that `x.shape[0] == 2` and return a linear interpolant. Would you be interested in sending a fix @DougGard ?
I've never looked into SciPy code, not sure I would be the best person for this.
No worries, we all started somewhere :-). 
The fix should go around this line, I'd think: https://github.com/scipy/scipy/blob/main/scipy/interpolate/_cubic.py#L510
When the inputs have length 2, `m` has length 5, and the interior points of `m` are computed in

https://github.com/scipy/scipy/blob/846ecd4648c3da29bfc1e05f7487c3fa047c94ca/scipy/interpolate/_cubic.py#L513-L513

which fills in just the single point `m[2]`.  But in
 
https://github.com/scipy/scipy/blob/846ecd4648c3da29bfc1e05f7487c3fa047c94ca/scipy/interpolate/_cubic.py#L516

the value at `m[3]` is referenced.

So it looks like the case where the length is 2 needs to be handled separately, as @ev-br suggested above.
I am able to reproduce the behaviour (the numbers are different for me),

```
-0.7803125
0.6648000000000001
0.2926
```

Commit message: TEST: Add 2D and 3D case

Patch to review:
diff --git a/scipy/interpolate/_cubic.py b/scipy/interpolate/_cubic.py
index fed25e254870..30771aa22e88 100644
--- a/scipy/interpolate/_cubic.py
+++ b/scipy/interpolate/_cubic.py
@@ -405,10 +405,10 @@ class Akima1DInterpolator(CubicHermiteSpline):
         .. versionadded:: 1.13.0
 
     extrapolate : {bool, None}, optional
-        If bool, determines whether to extrapolate to out-of-bounds points 
-        based on first and last intervals, or to return NaNs. If None, 
+        If bool, determines whether to extrapolate to out-of-bounds points
+        based on first and last intervals, or to return NaNs. If None,
         ``extrapolate`` is set to False.
-        
+
     Methods
     -------
     __call__
@@ -491,7 +491,7 @@ class Akima1DInterpolator(CubicHermiteSpline):
 
     """
 
-    def __init__(self, x, y, axis=0, *, method: Literal["akima", "makima"]="akima", 
+    def __init__(self, x, y, axis=0, *, method: Literal["akima", "makima"]="akima",
                  extrapolate:bool | None = None):
         if method not in {"akima", "makima"}:
             raise NotImplementedError(f"`method`={method} is unsupported.")
@@ -509,37 +509,45 @@ def __init__(self, x, y, axis=0, *, method: Literal["akima", "makima"]="akima",
         # Akima extrapolation historically False; parent class defaults to True.
         extrapolate = False if extrapolate is None else extrapolate
 
-        # determine slopes between breakpoints
-        m = np.empty((x.size + 3, ) + y.shape[1:])
-        dx = dx[(slice(None), ) + (None, ) * (y.ndim - 1)]
-        m[2:-2] = np.diff(y, axis=0) / dx
-
-        # add two additional points on the left ...
-        m[1] = 2. * m[2] - m[3]
-        m[0] = 2. * m[1] - m[2]
-        # ... and on the right
-        m[-2] = 2. * m[-3] - m[-4]
-        m[-1] = 2. * m[-2] - m[-3]
-
-        # if m1 == m2 != m3 == m4, the slope at the breakpoint is not
-        # defined. This is the fill value:
-        t = .5 * (m[3:] + m[:-3])
-        # get the denominator of the slope t
-        dm = np.abs(np.diff(m, axis=0))
-        if method == "makima":
-            pm = np.abs(m[1:] + m[:-1])
-            f1 = dm[2:] + 0.5 * pm[2:]
-            f2 = dm[:-2] + 0.5 * pm[:-2]
+        if y.shape[0] == 2:
+            # edge case: only have two points, use linear interpolation
+            xp = x.reshape((x.shape[0],) + (1,)*(y.ndim-1))
+            hk = xp[1:] - xp[:-1]
+            mk = (y[1:] - y[:-1]) / hk
+            t = np.zeros_like(y)
+            t[...] = mk
         else:
-            f1 = dm[2:]
-            f2 = dm[:-2]
-        f12 = f1 + f2
-        # These are the mask of where the slope at breakpoint is defined:
-        ind = np.nonzero(f12 > 1e-9 * np.max(f12, initial=-np.inf))
-        x_ind, y_ind = ind[0], ind[1:]
-        # Set the slope at breakpoint
-        t[ind] = (f1[ind] * m[(x_ind + 1,) + y_ind] +
-                  f2[ind] * m[(x_ind + 2,) + y_ind]) / f12[ind]
+            # determine slopes between breakpoints
+            m = np.empty((x.size + 3, ) + y.shape[1:])
+            dx = dx[(slice(None), ) + (None, ) * (y.ndim - 1)]
+            m[2:-2] = np.diff(y, axis=0) / dx
+
+            # add two additional points on the left ...
+            m[1] = 2. * m[2] - m[3]
+            m[0] = 2. * m[1] - m[2]
+            # ... and on the right
+            m[-2] = 2. * m[-3] - m[-4]
+            m[-1] = 2. * m[-2] - m[-3]
+
+            # if m1 == m2 != m3 == m4, the slope at the breakpoint is not
+            # defined. This is the fill value:
+            t = .5 * (m[3:] + m[:-3])
+            # get the denominator of the slope t
+            dm = np.abs(np.diff(m, axis=0))
+            if method == "makima":
+                pm = np.abs(m[1:] + m[:-1])
+                f1 = dm[2:] + 0.5 * pm[2:]
+                f2 = dm[:-2] + 0.5 * pm[:-2]
+            else:
+                f1 = dm[2:]
+                f2 = dm[:-2]
+            f12 = f1 + f2
+            # These are the mask of where the slope at breakpoint is defined:
+            ind = np.nonzero(f12 > 1e-9 * np.max(f12, initial=-np.inf))
+            x_ind, y_ind = ind[0], ind[1:]
+            # Set the slope at breakpoint
+            t[ind] = (f1[ind] * m[(x_ind + 1,) + y_ind] +
+                    f2[ind] * m[(x_ind + 2,) + y_ind]) / f12[ind]
 
         super().__init__(x, y, t, axis=0, extrapolate=extrapolate)
         self.axis = axis
diff --git a/scipy/interpolate/tests/test_interpolate.py b/scipy/interpolate/tests/test_interpolate.py
index 24a6907b7b05..17031f82d223 100644
--- a/scipy/interpolate/tests/test_interpolate.py
+++ b/scipy/interpolate/tests/test_interpolate.py
@@ -906,6 +906,60 @@ def test_eval_3d(self):
         yi[:, 1, 1] = 4. * yi_
         xp_assert_close(ak(xi), yi)
 
+    def test_linear_interpolant_edge_case(self):
+        # Reference - https://github.com/scipy/scipy/issues/22011
+        x = np.array([0.0, 1.0], dtype=float)
+        y = np.array([0.55, 0.55])
+        akima = Akima1DInterpolator(x, y, axis=0, extrapolate=None)
+        xp_assert_close(akima(0.45), np.array(0.55))
+
+        x = np.array([0.0, 0.5, 1.0])
+        y = np.array([0.15, 0.67, 0.15])
+        akima = Akima1DInterpolator(x, y, axis=0, extrapolate=None)
+        xp_assert_close(akima(0.45), np.array(0.6648))
+
+        x = np.array([0.0, 1.0])
+        y = np.array([0.55,  0.55])
+        akima = Akima1DInterpolator(x, y, axis=0, extrapolate=None)
+        xp_assert_close(akima(0.45), np.array(0.55))
+
+        # 2D case
+        x = np.array([0., 1.])
+        y = np.array([0., 1.])
+        y = np.column_stack((y, 2. * y, 3. * y, 4. * y))
+
+        ak = Akima1DInterpolator(x, y)
+        xi = np.array([0.5, 1.])
+        yi = np.array([[0.5, 1., 1.5, 2. ],
+                       [1., 2., 3., 4.]])
+        xp_assert_close(ak(xi), yi)
+
+        ak = Akima1DInterpolator(x, y.T, axis=1)
+        xp_assert_close(ak(xi), yi.T)
+
+        x = np.arange(0., 2.)
+        y_ = np.array([0., 1.])
+        y = np.empty((2, 2, 2))
+        y[:, 0, 0] = y_
+        y[:, 1, 0] = 2. * y_
+        y[:, 0, 1] = 3. * y_
+        y[:, 1, 1] = 4. * y_
+        ak = Akima1DInterpolator(x, y)
+        yi_ = np.array([0.5, 1.])
+        yi = np.empty((2, 2, 2))
+        yi[:, 0, 0] = yi_
+        yi[:, 1, 0] = 2. * yi_
+        yi[:, 0, 1] = 3. * yi_
+        yi[:, 1, 1] = 4. * yi_
+        xp_assert_close(ak(xi), yi)
+
+        ak = Akima1DInterpolator(x, y.transpose(1, 0, 2), axis=1)
+        xp_assert_close(ak(xi), yi.transpose(1, 0, 2))
+
+        ak = Akima1DInterpolator(x, y.transpose(2, 1, 0), axis=2)
+        xp_assert_close(ak(xi), yi.transpose(2, 1, 0))
+
+
     def test_degenerate_case_multidimensional(self):
         # This test is for issue #5683.
         x = np.array([0, 1, 2])


Code review comments:
diff_hunk: @@ -906,6 +906,60 @@ def test_eval_3d(self):
         yi[:, 1, 1] = 4. * yi_
         xp_assert_close(ak(xi), yi)
 
+    def test_linear_interpolant_edge_case(self):
+        # Reference - https://github.com/scipy/scipy/issues/22011
+        x = np.array([0.0, 1.0], dtype=float)
+        y = np.array([0.55, 0.55])
+        akima = Akima1DInterpolator(x, y, axis=0, extrapolate=None)
+        xp_assert_close(akima(0.45), np.array(0.55))
+
+        x = np.array([0.0, 0.5, 1.0])
+        y = np.array([0.15, 0.67, 0.15])
+        akima = Akima1DInterpolator(x, y, axis=0, extrapolate=None)
+        xp_assert_close(akima(0.45), np.array(0.6648))
+
+        x = np.array([0.0, 1.0])
+        y = np.array([0.55,  0.55])
+        akima = Akima1DInterpolator(x, y, axis=0, extrapolate=None)
+        xp_assert_close(akima(0.45), np.array(0.55))
+
+        # 2D case
+        x = np.array([0., 1.])
+        y = np.array([0., 1.])
+        y = np.column_stack((y, 2. * y, 3. * y, 4. * y))
path: scipy/interpolate/tests/test_interpolate.py
text: @user1:
```suggestion
        y = np.column_stack((x, 2. * x, 3. * x, 4. * x))
```

and remove `y = np.array([0., 1.0])` above

@author:
Ah! Yes. Nice. Thanks. I don't know why I didn't see it before. Thanks a lot. :-)