CODEWITHSUNDEEP
pythonscipyinterpolationlinear-interpolation
vice_spirit
vice_spirit

Reputation: 323

Scipys interp1d and infinities

When interpolating linearly with scipy.interpolate.interp1d I encounter the following behaviour:

from scipy.interpolate import interp1d
import numpy as np
x = [0,1,2,3,4]
y = [np.inf, 10, 12, 11, 9]
interpolated_function = interp1d(x,y)

Now of course the interpolation is not welldefined in [0,1), but at 1 i would expect it to be defined. However:

In[2]:  interpolated_function(1.)
        C:\WinPython27\python-2.7.10\lib\site-packages\
        scipy\interpolate\interpolate.py:469:
        RuntimeWarning: invalid value encountered in add
        y_new = slope*(x_new - x_lo)[:, None] + y_lo 
Out[2]: array(nan)

In[3]:  interpolated_function(1.000000000000001)
Out[3]: array(10.000000000000002)

Is this expected behaviour? Shouldn't the interpolated function at 1 be evaluated as 10, as this is a perfectly valid datapoint passed to interp1d?

Upvotes: 1

Views: 827

Answers (1)

ev-br
ev-br

Reputation: 26090

interp1d does not special-case nans or infs. For kind="linear" it simply uses the formula printed in the error message you got, on the interval which contains the input value.

Since exact equality is not reliable in floating point, by evaluating it at 1.0 you are relying on the implementation detail which classifies it to be on either (0, 1) or (1, 2).

For example,

In [26]: np.interp(1, x, y)
Out[26]: 10.0

In [32]: interp1d(x, y, kind='slinear')(1.)
Out[32]: array(10.0)

In the future versions of scipy, the behavior will depend on the fill_value. (In the current dev version you can use fill_value="extrapolate".)

The best thing to do is to filter out non-finite numbers on the user side. 

In [33]: x, y = map(np.asarray, (x, y))

In [34]: mask = np.isfinite(x) & np.isfinite(y)

In [35]: ii = interp1d(x[mask], y[mask])

In [36]: ii(1.)
Out[36]: array(10.0)

Upvotes: 2

Related Questions