Why does Matlab interp1 produce different results than numpy interp?

Question

EDIT: Code edited to produce results consistent with Matlab. See below.

I am converting Matlab scripts to Python and the linear interpolation results are different in certain cases. I wonder why and if there is any way to fix this?

Here is the code example in both Matlab and Python and the resulting output (Note that t just so happens to be equal to tin in this case):

MATLAB:

t= [ 736696., 736696.00208333, 736696.00416667, 736696.00625, 736696.00833333, 736696.01041667, 736696.0125];
tin =[ 736696., 736696.00208333, 736696.00416667, 736696.00625, 736696.00833333, 736696.01041667,  736696.0125];
xin = [   nan , 1392.,  1406. , 1418.  ,  nan , 1442. ,   nan];

interp1(tin,xin,t)

ans =

 NaN        1392        1406        1418         NaN        1442         NaN

Python (numpy):

(scipy interpolate.interp1d produces the same result as numpy)

t= [ 736696., 736696.00208333, 736696.00416667, 736696.00625, 736696.00833333, 736696.01041667, 736696.0125];
tin =[ 736696., 736696.00208333, 736696.00416667, 736696.00625, 736696.00833333, 736696.01041667,  736696.0125];
xin = [   nan , 1392.,  1406. , 1418.  ,  nan , 1442. ,   nan];

x = np.interp(t,tin,xin)

array([   nan,  1392.,  1406.,    nan,    nan,    nan,    nan]) 

# Edit
# Find indices where t == tin and if the np.interp output 
# does not match the xin array, overwrite the np.interp output at those 
# indices 
same = np.where(t == tin)[0]
not_same = np.where(xin[same] != x[same])[0]
x[not_same] = xin[not_same]

Answer 1

It appears as if Matlab includes an additional equality check in it's interpolation.

Linear 1-D interpolation is generally done by finding two x values which span the input value x and then calculating the result as:

y = y1 + (y2-y1)*(x-x1)/(x2-x1)

If you pass in an x value which is exactly equal to one of the input x coordinates, the routine will generally calculate the correct value since x-x1 will be zero. However, if your input array has a nan as y1 or y2 these will propagate to the result.

Based on the code you posted, my best guess would be that Matlab's interpolation function has an additional check that is something like:

if x == x1:
    return y1

and that the numpy function does not have this check.

To achieve the same effect in numpy you could do:

np.where(t == tin,xin,np.interp(t,tin,xin))