Repository: scipy/scipy
Pull request: BUG: add const for Read-only views in interpnd.pyx
<!-- 
Thanks for contributing a pull request! Please ensure that
your PR satisfies the checklist before submitting:
http://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://docs.scipy.org/doc/numpy/dev/development_workflow.html#writing-the-commit-message

Note that we are a team of volunteers; we appreciate your
patience during the review process.

Again, thanks for contributing!
-->

#### Reference issue
fix #12499 

#### What does this implement/fix?
I added const for Read-only views in interpnd.pyx and add a test for it.


#### Additional information
Ref: https://cython.readthedocs.io/en/latest/src/userguide/memoryviews.html#read-only-views

Problem statement:
LinearNDInterpolator raises ValueError when value array has writeable=False flag.
LinearNDInterpolator and CloughToucher2DInterppolator raise ValueError when value array has writeable=False flag. This is unanticipated and undocumented. Suggest that documentation or code should be changed. Given that they work with writeable
arrays, other bugs might result as well such as side-effects from code executed later on. 

#### Reproducing code example:
```
import numpy as np
from scipy.interpolate import LinearNDInterpolator, CloughTocher2DInterpolator

xs = np.arange(0, 4, dtype=float)
xs, ys = np.meshgrid(xs, xs)
values = np.ravel(xs + ys)
points = np.reshape([xs, ys], (2, -1)).T


print(70*'*')
print('LinearNDInterpolator')
interp_0 = LinearNDInterpolator(points, values)

print('interpolated values = :', interp_0(2, 2))
print('\nWith values.setflags(write=False), evaluation raises a ValueError')
values.setflags(write=False)
interp_1 = LinearNDInterpolator(points, values)
print('interpolated values = :', interp_1(2, 2))
```

#### Output and Error message:
<!-- If any, paste the *full* error message inside a code block
as above (starting from line Traceback)
-->

```
**********************************************************************
LinearNDInterpolator
interpolated values = : 4.0

With values.setflags(write=False), evaluation raises a ValueError
********************************************************************** 

Traceback (most recent call last):

  File "<ipython-input-2-fd2384f87746>", line 19, in <module>
    print('interpolated values = :', interp_1(2, 2))

  File "interpnd.pyx", line 151, in scipy.interpolate.interpnd.NDInterpolatorBase.__call__

  File "interpnd.pyx", line 251, in scipy.interpolate.interpnd.LinearNDInterpolator._evaluate_double

  File "interpnd.pyx", line 259, in scipy.interpolate.interpnd.LinearNDInterpolator._do_evaluate

  File "stringsource", line 658, in View.MemoryView.memoryview_cwrapper

  File "stringsource", line 349, in View.MemoryView.memoryview.__cinit__

ValueError: buffer source array is read-only
```

#### Scipy/Numpy/Python version information:
<!-- You can simply run the following and paste the result in a code block
```
import sys, scipy, numpy; print(scipy.__version__, numpy.__version__, sys.version_info)
```
--> 1.4.1 1.18.1 sys.version_info(major=3, minor=7, micro=7, releaselevel='final', serial=0)

Hints from issue discussion:
A similar issue arises with CloughToucher2DInterpolator, but with instantiation:

```
import numpy as np
from scipy.interpolate import LinearNDInterpolator, CloughTocher2DInterpolator

print(70*'*')
print('CloughTocher2DInterpolator')
xs = np.arange(0, 4, dtype=float)
xs, ys = np.meshgrid(xs, xs)
values = np.ravel(xs + ys)
points = np.reshape([xs, ys], (2, -1)).T
interp_0 = CloughTocher2DInterpolator(points, values)
print('interpolated values = :', interp_0(2, 2))

print('\nWith values.setflags(write=False), instantiation raises a ValueError')
print(70*'*', '\n')

values.setflags(write=False)
interp_1 = CloughTocher2DInterpolator(points, values)
```

#### Output and Error message
```
**********************************************************************
CloughTocher2DInterpolator
interpolated values = : 4.0

With values.setflags(write=False), instantiation raises a ValueError
********************************************************************** 

Traceback (most recent call last):

  File "<ipython-input-1-95325d66f264>", line 17, in <module>
    interp_1 = CloughTocher2DInterpolator(points, values)

  File "interpnd.pyx", line 847, in scipy.interpolate.interpnd.CloughTocher2DInterpolator.__init__

  File "interpnd.pyx", line 520, in scipy.interpolate.interpnd.estimate_gradients_2d_global

  File "stringsource", line 658, in View.MemoryView.memoryview_cwrapper

  File "stringsource", line 349, in View.MemoryView.memoryview.__cinit__

ValueError: buffer source array is read-only
```

Commit message: add const for Read-only views in interpnd.pyx and add a test

Patch to review:
diff --git a/scipy/interpolate/interpnd.pyx b/scipy/interpolate/interpnd.pyx
index 44707e11b93c..dac18d95cff7 100644
--- a/scipy/interpolate/interpnd.pyx
+++ b/scipy/interpolate/interpnd.pyx
@@ -257,7 +257,7 @@ class LinearNDInterpolator(NDInterpolatorBase):
     @cython.boundscheck(False)
     @cython.wraparound(False)
     def _do_evaluate(self, double[:,::1] xi, double_or_complex dummy):
-        cdef double_or_complex[:,::1] values = self.values
+        cdef const double_or_complex[:,::1] values = self.values
         cdef double_or_complex[:,::1] out
         cdef double[:,::1] points = self.points
         cdef int[:,::1] simplices = self.tri.simplices
@@ -491,7 +491,7 @@ cdef int _estimate_gradients_2d_global(qhull.DelaunayInfo_t *d, double *data,
 @cython.boundscheck(False)
 @cython.wraparound(False)
 cpdef estimate_gradients_2d_global(tri, y, int maxiter=400, double tol=1e-6):
-    cdef double[:,::1] data
+    cdef const double[:,::1] data
     cdef double[:,:,::1] grad
     cdef qhull.DelaunayInfo_t info
     cdef int k, ret, nvalues
@@ -857,7 +857,7 @@ class CloughTocher2DInterpolator(NDInterpolatorBase):
     @cython.boundscheck(False)
     @cython.wraparound(False)
     def _do_evaluate(self, double[:,::1] xi, double_or_complex dummy):
-        cdef double_or_complex[:,::1] values = self.values
+        cdef const double_or_complex[:,::1] values = self.values
         cdef double_or_complex[:,:,:] grad = self.grad
         cdef double_or_complex[:,::1] out
         cdef double[:,::1] points = self.points
diff --git a/scipy/interpolate/tests/test_interpolate.py b/scipy/interpolate/tests/test_interpolate.py
index 30497a56e3a5..988d4db294e3 100644
--- a/scipy/interpolate/tests/test_interpolate.py
+++ b/scipy/interpolate/tests/test_interpolate.py
@@ -2613,6 +2613,29 @@ def test_broadcastable_input(self):
                          interp_points3.size ==
                          interp_points4.size, True)
 
+    def test_read_only(self):
+        # input data
+        np.random.seed(0)
+        x = np.random.random(10)
+        y = np.random.random(10)
+        z = np.hypot(x, y)
+        
+        # x-y grid for interpolation
+        X = np.linspace(min(x), max(x))
+        Y = np.linspace(min(y), max(y))
+        X, Y = np.meshgrid(X, Y)
+
+        x.setflags(write=False)
+        y.setflags(write=False)
+        z.setflags(write=False)
+        X.setflags(write=False)
+        Y.setflags(write=False)
+
+        for interpolator in (NearestNDInterpolator, LinearNDInterpolator,
+                             CloughTocher2DInterpolator):
+            interp = interpolator(list(zip(x, y)), z)
+            interp((X, Y))
+
 
 class MyValue(object):
     """


Code review comments:
diff_hunk: @@ -2613,6 +2613,29 @@ def test_broadcastable_input(self):
                          interp_points3.size ==
                          interp_points4.size, True)
 
+    def test_read_only(self):
+        # input data
+        np.random.seed(0)
+        x = np.random.random(10)
+        y = np.random.random(10)
+        z = np.hypot(x, y)
+        
+        # x-y grid for interpolation
+        X = np.linspace(min(x), max(x))
+        Y = np.linspace(min(y), max(y))
+        X, Y = np.meshgrid(X, Y)
+
+        x.setflags(write=False)
+        y.setflags(write=False)
+        z.setflags(write=False)
+        X.setflags(write=False)
+        Y.setflags(write=False)
+
+        for interpolator in (NearestNDInterpolator, LinearNDInterpolator,
+                             CloughTocher2DInterpolator):
+            interp = interpolator(list(zip(x, y)), z)
+            interp((X, Y))
path: scipy/interpolate/tests/test_interpolate.py
text: @user1:
Isn't this going to create a new array internally which will be writeable. I think that's the reason the test is passing even though `xi` isn't marked `const`.

@user1:
Just to be clear, by "internally" I mean the interpolator will effectively call `np.ascontiguousarray((X, Y))` which will make the new array. The same should happen for `list(zip(x, y))` above where it will be turned into a single ndarray.

So, I would say make sure you're passing contiguous ndarrays for all arguments in this test.

@author:
Hi. thank you for reviewing and catching it. you are right, we need more const. I fixed it. PTAL.