Quadratic interpolate.interpo1d() function does not work with read-only array input

[will133]

The interpolate.interpo1d() function with kind='quadratic' does not work with read-only array input. If kind is not given (default) or if the array is not read-only, the function would work properly. For instance:

Reproducing code example:

In [31]: x = np.arange(0, 10)
In [32]: y = np.exp(-x/3.0)
In [33]: f = interpolate.interp1d(x, y, kind='quadratic')
In [34]: xnew = np.arange(0, 9, 0.1)

# this would work fine:
In[35]: ok = f(xnew)

# this is what would cause problem
In [36]: xnew.flags.writeable = False
In [37]: f(xnew)

Error message:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-37-c7d3ba7c04f3> in <module>()
----> 1 f(xnew)

../lib/python2.7/site-packages/scipy/interpolate/polyint.pyc in __call__(self, x)
     77         """
     78         x, x_shape = self._prepare_x(x)
---> 79         y = self._evaluate(x)
     80         return self._finish_y(y, x_shape)
     81

../lib/python2.7/site-packages/scipy/interpolate/interpolate.pyc in _evaluate(self, x_new)
    660         #    The behavior is set by the bounds_error variable.
    661         x_new = asarray(x_new)
--> 662         y_new = self._call(self, x_new)
    663         if not self._extrapolate:
    664             below_bounds, above_bounds = self._check_bounds(x_new)

../lib/python2.7/site-packages/scipy/interpolate/interpolate.pyc in _call_spline(self, x_new)
    648
    649     def _call_spline(self, x_new):
--> 650         return self._spline(x_new)
    651
    652     def _call_nan_spline(self, x_new):

../lib/python2.7/site-packages/scipy/interpolate/_bsplines.pyc in __call__(self, x, nu, extrapolate)
    348         out = np.empty((len(x), prod(self.c.shape[1:])), dtype=self.c.dtype)
    349         self._ensure_c_contiguous()
--> 350         self._evaluate(x, nu, extrapolate, out)
    351         out = out.reshape(x_shape + self.c.shape[1:])
    352         if self.axis != 0:

../lib/python2.7/site-packages/scipy/interpolate/_bsplines.pyc in _evaluate(self, xp, nu, extrapolate, out)
    359     def _evaluate(self, xp, nu, extrapolate, out):
    360         _bspl.evaluate_spline(self.t, self.c.reshape(self.c.shape[0], -1),
--> 361                 self.k, xp, nu, extrapolate, out)
    362
    363     def _ensure_c_contiguous(self):

_bspl.pyx in scipy.interpolate._bspl.evaluate_spline()

../lib/python2.7/site-packages/scipy/interpolate/_bspl.so in View.MemoryView.memoryview_cwrapper()

../lib/python2.7/site-packages/scipy/interpolate/_bspl.so in View.MemoryView.memoryview.__cinit__()

ValueError: buffer source array is read-only

Scipy/Numpy/Python version information:

In [42]: import sys, scipy, numpy; print(scipy.__version__, numpy.__version__, sys.version_info)
('1.1.0', '1.11.3', sys.version_info(major=2, minor=7, micro=12, releaselevel='final', serial=0))

[ev-br]

Cython memoryviews seem to have issues with read-only numpy arrays. Here's a repro without scipy:

In [14]: %load_ext cython

In [15]: %%cython
    ...: def cyfunc(double[:] m):
    ...:      cdef double x = m[0]
    ...:      return x
    ...: 

In [16]: arr = np.arange(5, dtype=float)

In [17]: cyfunc(arr)
Out[17]: 0.0

In [18]: arr.flags.writeable = False

In [19]: cyfunc(arr)
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-19-5fcf616d0a4f> in <module>()
----> 1 cyfunc(arr)

_cython_magic_19a8015f59531c8d68f20d9468335fbf.pyx in _cython_magic_19a8015f59531c8d68f20d9468335fbf.cyfunc()

~/.cache/ipython/cython/_cython_magic_19a8015f59531c8d68f20d9468335fbf.cpython-35m-x86_64-linux-gnu.so in View.MemoryView.memoryview_cwrapper()

~/.cache/ipython/cython/_cython_magic_19a8015f59531c8d68f20d9468335fbf.cpython-35m-x86_64-linux-gnu.so in View.MemoryView.memoryview.__cinit__()

ValueError: buffer source array is read-only

While arguably this is an implementation detail (typed memoryviews) leaking out of scipy.interpolate, this is best fixed on the cython side. No idea if it's easy to fix in cython or not.
Going away from typed memoryviews in scipy.interpolate is certainly not desirable.
Best report it to Cython?

[will133]

I think there's some discussion on the cython github:

cython/cython#1605

It looks like this patch is already merged in?

cython/cython#1869

I'm still going through what exactly is that we need to do though. Using a const keyword maybe?

[will133]

I tried the following, which seems to work:

In [9]: import cython

In [10]: cython.__version__
'0.28.5'

In [11]: %%cython
    ...: def cyfunc(const double[:] m):
    ...:     cdef double x = m[0]
    ...:     return x
    ...:



In [12]: cyfunc(arr)
0.0

In [13]: arr.flags.writeable
False

[ev-br]

Then it's the matter of bumping up the Cython version to 0.28.x

Would you be interested in sending a PR?
I'd suggest:

• add the OP example as a test. It's going to fail currently
• bump the Cython version is CI: .travis.yml, .appveyor.yml and Circle CI.
• if all of the CI is green, we're done.

[will133]

I think the trick is to add the const keyword to various functions in _bspl.pyx file where the array is passed in as not modifiable. I started doing that for things like "double[::1]" and turn that into "const double[::1]". I did have problem doing "const double_or_complex[:,::1]" though. Where would the test go?

[ev-br]

Tests are in scipy/interpolate/tests : test_bsplines.py for b-splines, test_interpolate.py for interp1d and PPoly (const memoryviews are best added simultaneously to b-splines and PPoly).

[ev-br]

cython/cython#1772 mentions const memoryviews and fused types. So it might be that we'd want to wait until this is ironed out in cython (which has happened previously).

[will133]

I think I can get my use case to work when I add a couple of const to the declaration in _bspl.pyx (ignoring the fused type problem for double_or_complex problem for now). It certainly does not cover all the cases where this can happen in scipy though. I think for some of those functions the proposed workaround does not work. For instance, you have double_or_complex as writable and readable in the same function signature in _norm_eq_lsq, so we would probably hit the corner case where cython would fail.

I can look into adding a few tests. Do you still want a pull request on that?