readonly_interpolate_PR17717

===================================================================

Gbp-Pq: Name readonly_interpolate_PR17717.patch

diff --git a/scipy/interpolate/_rgi.py b/scipy/interpolate/_rgi.py
index 37c76bb0b710a5e9b9270ac68f61a50eca741fb9..ec9f55f76aad233fde0ef79f5340763ee7fa6fab 100644 (file)
--- a/scipy/interpolate/_rgi.py
+++ b/scipy/interpolate/_rgi.py
@@ -23,11 +23,12 @@ def _check_points(points):
                 # input is descending, so make it ascending
                 descending_dimensions.append(i)
                 p = np.flip(p)
-                p = np.ascontiguousarray(p)
             else:
                 raise ValueError(
                     "The points in dimension %d must be strictly "
                     "ascending or descending" % i)
+        # see https://github.com/scipy/scipy/issues/17716
+        p = np.ascontiguousarray(p)
         grid.append(p)
     return tuple(grid), tuple(descending_dimensions)
 
@@ -330,7 +331,10 @@ class RegularGridInterpolator:
         if method == "linear":
             indices, norm_distances = self._find_indices(xi.T)
             if (ndim == 2 and hasattr(self.values, 'dtype') and
-                    self.values.ndim == 2):
+                    self.values.ndim == 2 and self.values.flags.writeable and
+                    self.values.dtype.byteorder == '='):
+                # until cython supports const fused types, the fast path
+                # cannot support non-writeable values
                 # a fast path
                 out = np.empty(indices.shape[1], dtype=self.values.dtype)
                 result = evaluate_linear_2d(self.values,

diff --git a/scipy/interpolate/_rgi_cython.pyx b/scipy/interpolate/_rgi_cython.pyx
index 22d4a2eda7afe8439cb280349779cc663b5dafc4..71e030664b75196bbf59a2f661acdaafb4585aec 100644 (file)
--- a/scipy/interpolate/_rgi_cython.pyx
+++ b/scipy/interpolate/_rgi_cython.pyx
@@ -17,7 +17,7 @@ np.import_array()
 @cython.boundscheck(False)
 @cython.initializedcheck(False)
 def evaluate_linear_2d(double_or_complex[:, :] values, # cannot declare as ::1
-                       long[:, :] indices,             # unless prior
+                       const long[:, :] indices,       # unless prior
                        double[:, :] norm_distances,    # np.ascontiguousarray
                        tuple grid not None,
                        double_or_complex[:] out):
@@ -72,11 +72,13 @@ def evaluate_linear_2d(double_or_complex[:, :] values, # cannot declare as ::1
 @cython.boundscheck(False)
 @cython.cdivision(True)
 @cython.initializedcheck(False)
-def find_indices(tuple grid not None, double[:, :] xi):
+def find_indices(tuple grid not None, const double[:, :] xi):
+    # const is required for xi above in case xi is read-only
     cdef:
         long i, j, grid_i_size
         double denom, value
-        double[::1] grid_i
+        # const is required in case grid is read-only
+        const double[::1] grid_i
 
         # Axes to iterate over
         long I = xi.shape[0]

diff --git a/scipy/interpolate/tests/test_rgi.py b/scipy/interpolate/tests/test_rgi.py
index fc255d946f2f2cab19debc98589f392e2514e719..cb2c54cf3701f41e989104c48a98ba19974efd97 100644 (file)
--- a/scipy/interpolate/tests/test_rgi.py
+++ b/scipy/interpolate/tests/test_rgi.py
@@ -933,3 +933,58 @@ class TestInterpN:
                "RegularGridInterpolator has dimension 1")
         with assert_raises(ValueError, match=msg):
             interpn(points, values, xi)
+
+    def test_readonly_grid(self):
+        # https://github.com/scipy/scipy/issues/17716
+        x = np.linspace(0, 4, 5)
+        y = np.linspace(0, 5, 6)
+        z = np.linspace(0, 6, 7)
+        points = (x, y, z)
+        values = np.ones((5, 6, 7))
+        point = np.array([2.21, 3.12, 1.15])
+        for d in points:
+            d.flags.writeable = False
+        values.flags.writeable = False
+        point.flags.writeable = False
+        interpn(points, values, point)
+        RegularGridInterpolator(points, values)(point)
+
+    def test_2d_readonly_grid(self):
+        # https://github.com/scipy/scipy/issues/17716
+        # test special 2d case
+        x = np.linspace(0, 4, 5)
+        y = np.linspace(0, 5, 6)
+        points = (x, y)
+        values = np.ones((5, 6))
+        point = np.array([2.21, 3.12])
+        for d in points:
+            d.flags.writeable = False
+        values.flags.writeable = False
+        point.flags.writeable = False
+        interpn(points, values, point)
+        RegularGridInterpolator(points, values)(point)
+
+    def test_non_c_contiguous_grid(self):
+        # https://github.com/scipy/scipy/issues/17716
+        x = np.linspace(0, 4, 5)
+        x = np.vstack((x, np.empty_like(x))).T.copy()[:, 0]
+        assert not x.flags.c_contiguous
+        y = np.linspace(0, 5, 6)
+        z = np.linspace(0, 6, 7)
+        points = (x, y, z)
+        values = np.ones((5, 6, 7))
+        point = np.array([2.21, 3.12, 1.15])
+        interpn(points, values, point)
+        RegularGridInterpolator(points, values)(point)
+
+    @pytest.mark.parametrize("dtype", ['>f8', '<f8'])
+    def test_endianness(self, dtype):
+        # https://github.com/scipy/scipy/issues/17716
+        # test special 2d case
+        x = np.linspace(0, 4, 5, dtype=dtype)
+        y = np.linspace(0, 5, 6, dtype=dtype)
+        points = (x, y)
+        values = np.ones((5, 6), dtype=dtype)
+        point = np.array([2.21, 3.12], dtype=dtype)
+        interpn(points, values, point)
+        RegularGridInterpolator(points, values)(point)