Is this an error in scipy.interpolate.interp1d?

Question

When using zero order interpolation I found that the last Y value in the input array is not returned for the last value in the X array:

from scipy.interpolate import interp1d

xx = [0.0, 1.0, 2.0]
xi = interp1d(xx, xx, kind='zero')
print(xi(xx))

seems like it should return [0., 1., 2.] but it returns [0., 1., 1.]. The last value in xx is considered to be in the interpolation range but is not returned as the value of the last point. The documentation does not provide details for 'zero' but I would expect that it would either:

a) raise a ValueError because the input values are considered to define values on the semi-closed ranges [0., 1.) and [1., 2.), thereby leaving 2.0 undefined, or

b) return 2.0 because the ranges are [0., 1.), [1., 2.) and [2., 2.)

The interp1d function appears to consider the correct answer to be:

c) return 1.0 because the last range is a special case, defined as a closed interval [1., 2.]

Is there a correct choice? If so, is it what is implemented by interp1d?

Answer 1

The zero order spline is piecewise constant, and has discontinuities at the knots, which here are the interpolation points, so xi(1.0-1e-13) == 0 and xi(1.0+1e-13) == 1.

The interpolation interval in interp1d is defined as closed, [0, 2]. In principle one would expect that there is a single floating point value, x=2.0, which gives the result 2.0

However, as noted in the comments above, the spline implementation here is from FITPACK, which defines the k=0 spline as continuous at the knots from the right, except for the last interval which is different. I do not know the reason --- the Fortran code dates from the 80s. My guess would be there is no specific reason for the way it works like this, potentially apart from the fact that it was slightly more convenient to write the code like this with the B-spline representation.

This behavior in my opinion is a bug/quirk, but since the presence of any rounding error in x-values makes its impact zero, it's been way down in priority to address. (One aspect to consider is also that since there is no specification of what it does, it's defined by the implementation; breaking backward compat may be worse than the issue itself.)

EDIT: as noted in the other answer, the spline is actually constructed by splmake; this routine is not from FITPACK. Can't say without looking whether the final interval behavior is due to fitting or construction.

Answer 2

According to the source the interval is between [0, 2], including 2. The return from creating the spline is the cvals splmake(xx, xx, 0)[1] --> array([0.0, 1.0]). Following the source through we'll get to a call to evaluate the spline at a point which results in spleval(splmake(xx, xx, 0), 2) --> array(1.0). To answer your question, there's not really an answer it was just implemented this way, and if you think if the original array as encompassing the valid range I think it makes sense - but perhaps there should be a note about this particular evaluation, in which case you could always submit a request here. I can't comment otherwise I would but hopefully it answers your question.