Repository: scipy/scipy
Pull request: BUG: interpolate.RegularGridInterpolator: fix linear interp. on 2D grid with a single point
#### Reference issue
Closes GH-23171

#### What does this implement/fix?
`RegularGridInterpolator` (method=linear) works also in the degenerate case of single-points grids.
However in the specific case of a 2D, float64 grid a Nan is returned.

Tests for generic single-point grids are added and the spurious case is fixed.

#### Additional information
Changes are kept to a bare minimum, however unnecessary branches (where `out[point] = NAN` is set) were simply discarded.

Problem statement:
BUG: `interpolate.RegularGridInterpolator` fails on a 2D grid with a single point.
### Describe your issue.

This is a very specific edge case, which I found by running a parametrised code.

In general `RegularGridInterpolator(..., method="linear")` works as expected even if the number of grid points along one ore even all dimensions is one:
```python
In [3]: RegularGridInterpolator(
   ...:     points=([0], [0], [0]),
   ...:     values=[[[0]]],
   ...:     method="linear",
   ...: )([[0,0,0]])
Out[3]: array([0.])
```

This however fails for 2D grids, when the fast path `evaluate_linear_2d` is taken. In this case the interpolator returns `nan` instead of the correct value:
```python
In [10]: rgi = RegularGridInterpolator(
    ...:     points=([0], [0]),
    ...:     values=np.array([[0]]),
    ...:     method="linear",
    ...: )

In [11]: rgi([[0,0]])
Out[11]: array([nan])
```

The fact that the fast path is the culprit is confirmed by this example
```python
In [8]: rgi = RegularGridInterpolator(
   ...:     points=([0], [0]),
   ...:     values=np.array([[0]], dtype=np.float32),
   ...:     method="linear",
   ...: )

In [9]: rgi([[0,0]])
Out[9]: array([0.])
```
where `values` with `np.float32` type forces the generic path.

I understand that this is a very specific edge case, nevertheless the inconsistency between the generic and the fast 2D path should be addressed, in my opinion.

I can open a PR, since the fix appears to be trivial.

### Reproducing Code Example

```python
from scipy.interpolate import RegularGridInterpolator

rgi = RegularGridInterpolator(
    points=([0], [0]),
    values=[[0]],
    method="linear",
)
print(rgi([[0,0]]))
```

### Error message

```shell
No error, output is 
[nan]
```

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

```shell
1.17.0.dev0+git20250617.4c1a4da 2.3.0 sys.version_info(major=3, minor=13, micro=3, releaselevel='final', serial=0)
Build Dependencies:
  blas:
    detection method: extraframeworks
    found: true
    include directory: unknown
    lib directory: unknown
    name: Accelerate
    openblas configuration: unknown
    pc file directory: unknown
    version: unknown
  lapack:
    detection method: extraframeworks
    found: true
    include directory: unknown
    lib directory: unknown
    name: Accelerate
    openblas configuration: unknown
    pc file directory: unknown
    version: unknown
  pybind11:
    detection method: config-tool
    include directory: unknown
    name: pybind11
    version: 2.13.6
Compilers:
  c:
    commands: cc
    linker: ld64
    name: clang
    version: 17.0.0
  c++:
    commands: c++
    linker: ld64
    name: clang
    version: 17.0.0
  cython:
    commands: cython
    linker: cython
    name: cython
    version: 3.1.2
  fortran:
    commands: gfortran
    linker: ld64
    name: gcc
    version: 15.1.0
  pythran:
    include directory: ../.venv/lib/python3.13/site-packages/pythran
    version: 0.18.0
Machine Information:
  build:
    cpu: aarch64
    endian: little
    family: aarch64
    system: darwin
  cross-compiled: false
  host:
    cpu: aarch64
    endian: little
    family: aarch64
    system: darwin
Python Information:
  path: /Users/miccoli/work/miccoli/scipy-fork/.venv/bin/python3
  version: '3.13'
```

Hints from issue discussion:
> I can open a PR, since the fix appears to be trivial.

Superb, please do!

That fast path was meant to be temporary until we find a better general solution. We now have the infra for that general solution, in `method="slinear"`, now only need to find the bandwidth to implement it. No reason to not fix the inconsistency meanwhile. Am looking forward to reviewing a PR. 
> Superb, please do!

OK, I will start crafting a PR.
For reference only, some old issues defining the behaviour for axes with a single point:
- #5890
- #7560

Commit message: BUG:interpolate.RegularGridInterpolator: fix linear interp. on 2D grid with a single point

Explicitly handle the case in which there is a single grid point on both
axes.

Branches in which `out[point] = NAN` is set are deleted because
redundant: out-of-bounds points/nans are handled by caller.

Patch to review:
diff --git a/scipy/interpolate/_rgi_cython.pyx b/scipy/interpolate/_rgi_cython.pyx
index cfbc8a627384..80c17459657a 100644
--- a/scipy/interpolate/_rgi_cython.pyx
+++ b/scipy/interpolate/_rgi_cython.pyx
@@ -27,29 +27,29 @@ def evaluate_linear_2d(const double_or_complex[:, :] values, # cannot declare as
         double_or_complex y0, y1, result
     assert out.shape[0] == num_points
 
-    if grid[1].shape[0] == 1:
-        # linear interpolation along axis=0
+    if grid[0].shape[0] == grid[1].shape[0] == 1:
+        # single grid point along both axes: no interpolation
+        for point in range(num_points):
+            assert indices[0, point] == indices[1, point] == -1
+            out[point] = values[0, 0]
+    elif grid[1].shape[0] == 1:
+        # single grid line along axis=1: linear interpolation along axis=0 alone
         for point in range(num_points):
             i0 = indices[0, point]
             if i0 >= 0:
                 y0 = norm_distances[0, point]
                 result = values[i0, 0]*(1 - y0) + values[i0+1, 0]*y0
                 out[point] = result
-            else:
-                # xi was nan: find_interval returns -1
-                out[point] = NAN
     elif grid[0].shape[0] == 1:
-        # linear interpolation along axis=1
+        # single grid line along axis=0: linear interpolation along axis=1 alone
         for point in range(num_points):
             i1 = indices[1, point]
             if i1 >= 0:
                 y1 = norm_distances[1, point]
                 result = values[0, i1]*(1 - y1) + values[0, i1+1]*y1
                 out[point] = result
-            else:
-                # xi was nan: find_interval returns -1
-                out[point] = NAN
     else:
+        # linear interpolation along both axes
         for point in range(num_points):
             i0, i1 = indices[0, point], indices[1, point]
             if i0 >=0 and i1 >=0:
@@ -61,9 +61,6 @@ def evaluate_linear_2d(const double_or_complex[:, :] values, # cannot declare as
                 result = result + values[i0+1, i1] * y0 * (1 - y1)
                 result = result + values[i0+1, i1+1] * y0 * y1
                 out[point] = result
-            else:
-                # xi was nan
-                out[point] = NAN
 
     return np.asarray(out)
 
diff --git a/scipy/interpolate/tests/test_rgi.py b/scipy/interpolate/tests/test_rgi.py
index 6ec922093244..ef97d1f3f3f6 100644
--- a/scipy/interpolate/tests/test_rgi.py
+++ b/scipy/interpolate/tests/test_rgi.py
@@ -3,7 +3,7 @@
 import pytest
 import numpy as np
 
-from numpy.testing import assert_warns
+from numpy.testing import assert_array_equal, assert_warns
 from numpy.exceptions import ComplexWarning
 
 from scipy._lib._array_api import (
@@ -335,6 +335,41 @@ def test_fillvalue_type(self):
         RegularGridInterpolator(points, values)
         RegularGridInterpolator(points, values, fill_value=0.)
 
+    @pytest.mark.parametrize("dtype", [np.float32, np.float64])
+    @pytest.mark.parametrize("ndim", [1, 2, 3])
+    @pytest.mark.parametrize("method", ["linear", "nearest"])
+    def test_length_one_axes(self, dtype, ndim, method):
+        # gh-23171: length-1 axes in all dimensions are legal
+        # for all methods parametrized above
+        xi = np.array([0], dtype=dtype)
+        points = (xi, ) * ndim
+        val = dtype(1/7)
+        fill = dtype(1/42)
+        nshape = (1, ) * ndim
+        values = np.array(val, dtype=dtype)
+        values.shape = nshape
+        interp = RegularGridInterpolator(
+            points, values, method=method, bounds_error=False, fill_value=fill
+        )
+
+        # grid point
+        x0 = np.concatenate(points)
+        # all NaN point
+        x_nan = np.array([np.nan, ]*ndim, dtype=dtype)
+        # single NaN point
+        x_nan_0 = np.copy(x0)
+        x_nan_0[0] = np.nan
+        # out-of-bounds point
+        x_oob = np.copy(x0)
+        x_oob[0] += 1
+
+        sample = np.array([x0, x_nan, x_nan_0, x_oob], dtype=dtype)
+        wanted = np.array([val, np.nan, np.nan, fill], dtype=dtype)
+        result = interp(sample)
+        assert_array_equal(result, wanted)
+
+
+
     def test_length_one_axis(self):
         # gh-5890, gh-9524 : length-1 axis is legal for method='linear'.
         # Along the axis it's linear interpolation; away from the length-1


Code review comments:
diff_hunk: @@ -27,29 +27,29 @@ def evaluate_linear_2d(const double_or_complex[:, :] values, # cannot declare as
         double_or_complex y0, y1, result
     assert out.shape[0] == num_points
 
-    if grid[1].shape[0] == 1:
-        # linear interpolation along axis=0
+    if grid[0].shape[0] == grid[1].shape[0] == 1:
+        # single grid point along both axes: no interpolation
+        for point in range(num_points):
+            assert indices[0, point] == indices[1, point] == -1
path: scipy/interpolate/_rgi_cython.pyx
text: @user1:
We usually try to avoid plain `assert` in our source code proper (it is "ok" in `pytest` tests) for a few reasons, including its omission in optimized CPython runs.

It probably isn't too much work to i.e., raise a custom `ValueError` if this check is really needed (and possibly add a regression test for the raising of that error as well, if there is a straightforward way to trigger it).

@user2:
Agreed, we should not be using `assert` for checking things. 

I'm also confused by the logic, actually. This piece of code should not raise at all. Previously, the logic was roughly "if the grid size in this direction is 1, return exactly the corresponding value or nan". I'd expect that the doubly degenerate case might need an extra check, but your patch below removes `out[point] = NAN` for extrapolated values.



@author:
The above assert should never raise: it is added only to explicitly state an invariant that is assumed to be true by the current implementation.

I understand that such use of asserts is not common/idiomatic in python... but I felt authorized to add my assert looking at line 28, where there is indeed an assert of this very same nature:

https://github.com/scipy/scipy/blob/4c1a4da8110cc709dc2ff5050d882318fc649e20/scipy/interpolate/_rgi_cython.pyx#L19-L29

When learning C I was taught that asserts should be used only for debugging, and that it is good practice to leave them also in production code (which is typically compiled with `-DNDEBUG`).
The reasons:
- explicitly document the invariants which are assumed to be valid in some point of the code
- help debugging regressions if for some reason the invariants fail in some future release.

So actually, the above assert should never raise and could be simply deleted or commened out.
Please let me know if I should 
- delete this assert (but then also the one at line 28), 
- add a message that clarifies its debugging-only use
- comment it out.


@user2:
> delete this assert (but then also the one at line 28),

This please. 

This whole thing is quite firmly in the "tastes differ" territory and indeed there might be very good reasons to keep asserts for self-documenting reasons, either in python or C. Here though, having both these asserts is more confusing then helpful --- and I'm to blame for the one in line 28, so thank you for cleaning up the mess I left in this code!

My reasoning for leaving the assert on line 28 was probably along the same lines of reasoning you discuss, but then a couple of years later it confuses multiple people, which means it was a bad idea.

diff_hunk: @@ -335,6 +335,41 @@ def test_fillvalue_type(self):
         RegularGridInterpolator(points, values)
         RegularGridInterpolator(points, values, fill_value=0.)
 
+    @pytest.mark.parametrize("dtype", [np.float32, np.float64])
+    @pytest.mark.parametrize("ndim", [1, 2, 3])
+    @pytest.mark.parametrize("method", ["linear", "nearest"])
+    def test_length_one_axes(self, dtype, ndim, method):
+        # gh-23171: length-1 axes in all dimensions are legal
+        # for all methods parametrized above
+        xi = np.array([0], dtype=dtype)
+        points = (xi, ) * ndim
+        val = dtype(1/7)
+        fill = dtype(1/42)
+        nshape = (1, ) * ndim
+        values = np.array(val, dtype=dtype)
+        values.shape = nshape
path: scipy/interpolate/tests/test_rgi.py
text: @user1:
perhaps a bit more compact/idiomatic to avoid manually setting the shape of the array?

<details>

```diff
--- a/scipy/interpolate/tests/test_rgi.py
+++ b/scipy/interpolate/tests/test_rgi.py
@@ -343,11 +343,9 @@ class TestRegularGridInterpolator:
         # for all methods parametrized above
         xi = np.array([0], dtype=dtype)
         points = (xi, ) * ndim
-        val = dtype(1/7)
-        fill = dtype(1/42)
-        nshape = (1, ) * ndim
-        values = np.array(val, dtype=dtype)
-        values.shape = nshape
+        val = np.asarray(1 / 7, dtype=dtype)
+        fill = np.asarray(1 / 42, dtype=dtype)
+        values = np.full((1,) * ndim, val, dtype=dtype)
```

</details>


@author:
Here setting the shape directly stresses the fact that, regardless of the number of dimensions `ndim`, the `values` array should have size=1: so it is first constructed as a 0D array which is promoted to `ndim`D by assigning the shape.

Similarly `val` and `fill` are meant to be numpy scalars and not 0D arrays... 

Currently numpy is a little bit confusing in the treatment of 0D arrays vs. scalars, so I prefer to be 100% explicit on what is meant.

@user2:
numpy scalars vs 0D arrays are confusing indeed! That said, the way it's written now makes it more confusing then just creating the right-shaped array (have I mentioned that we're in the "matters of taste" territory already?)

> val and fill are meant to be numpy scalars and not 0D arrays

why though? You want to reject fill_value and values being numpy scalars? Then please add an explicit assertion that it raises. I wouldn't though, exactly because it's confusing in numpy, and it's not always easy to predict if something is a scalar or a 0D array. (and in the spirit of array API, consider actually just calling `np.asarray` on something that might be a scalar)


@author:
I thought `val` (scalar) as a building block for `values` (array). Similarly with `fill` (scalar), I build the array of desired values in `wanted`.

Why? Actually there is no compelling reason, apart from the fact that I seems more natural to me to build N-dimensional arrays from scalars, and not 0-dimensional arrays.

Again I see here that this is not idiomatic, so I will rewrite the tests along your lines.

diff_hunk: @@ -335,6 +335,41 @@ def test_fillvalue_type(self):
         RegularGridInterpolator(points, values)
         RegularGridInterpolator(points, values, fill_value=0.)
 
+    @pytest.mark.parametrize("dtype", [np.float32, np.float64])
+    @pytest.mark.parametrize("ndim", [1, 2, 3])
+    @pytest.mark.parametrize("method", ["linear", "nearest"])
+    def test_length_one_axes(self, dtype, ndim, method):
+        # gh-23171: length-1 axes in all dimensions are legal
+        # for all methods parametrized above
+        xi = np.array([0], dtype=dtype)
+        points = (xi, ) * ndim
+        val = dtype(1/7)
+        fill = dtype(1/42)
+        nshape = (1, ) * ndim
+        values = np.array(val, dtype=dtype)
+        values.shape = nshape
+        interp = RegularGridInterpolator(
+            points, values, method=method, bounds_error=False, fill_value=fill
+        )
+
+        # grid point
+        x0 = np.concatenate(points)
path: scipy/interpolate/tests/test_rgi.py
text: @user1:
Nit: why is it not `points = np.asarray([0]*ndim, dtype=dtype)`

diff_hunk: @@ -335,6 +335,41 @@ def test_fillvalue_type(self):
         RegularGridInterpolator(points, values)
         RegularGridInterpolator(points, values, fill_value=0.)
 
+    @pytest.mark.parametrize("dtype", [np.float32, np.float64])
+    @pytest.mark.parametrize("ndim", [1, 2, 3])
+    @pytest.mark.parametrize("method", ["linear", "nearest"])
+    def test_length_one_axes(self, dtype, ndim, method):
path: scipy/interpolate/tests/test_rgi.py
text: @user1:
Nit: could you rename the test please (took me a while to register the difference between `length_one_axes` and `length_one_axis` below).

@author:
OK, I will rename the new test... 

@author:
Just a question of style: some projects implement a strict requirement of a single assert per unit test (this to have a single point of failure... instead of a sequnce of `assert A`, some code, `assert B`, some other code `assert C`).

Personally I like the idea of having a single assert, and indeed this was the reason for dumping all the tests in a single vector.

Should I split the different cases as multiple asserts in the same method, or should I write multiple methods with a single assert?

@user1:
This is even further into the "tastes differ" territory indeed. We certainly do not have a strict requirement of a single assertion for a single test. In the opposite limit, having a single test with 25 assertions is not nice either. Where to draw a line is a matter of taste in the end. Here though, you've what three or four assertions, and it seems at least to me that having separate assertions makes the test more readable (and saves a step in debugging if it ever fails).



@author:
OK, clear enough, hope to fix this soon.