Vectorized linear interpolation

simple_interpolation.py

# The kernel function is 
#   roughtly equivalent to new[:] = np.interp(xnew, xvals, yvals)
#   But this is broadcast so that it can be run many, many times 
#   quickly.

# Call using ynew = interp1d(xnew, xdata, ydata)
#   ynew.shape will be xnew.shape
#   Also, ydata.shape[-1] must be xdata.shape[-1]
#   and if ydata or xdata have ndim greater than 1, the initial dimensions
#   must be xnew.shape:
#     xnew.shape = (...)
#     xvals.shape = (n,) or (...,n)
#     yvals.shape = (n,) or (...,n)
#     ynew.shape = (...)
#
#   assumes xvals are sorted along last dimension
#   assumes there are no NaNs in xvals or yvals

# Simple 1-d linear interpolation 
@guvectorize(['float64[:], float64[:], float64[:], float64[:]'],
             "(),(n),(n) -> ()")
def interp1d(xnew, xvals, yvals, ynew):
    i = 0
    N = len(xvals)
    if xnew < xvals[0]:
        x_a = 0.0
        y_a = 0.0
        x_b = xvals[0]
        y_b = yvals[0]
    else:
        while i > 0 and i < N and xnew >= xvals[i]:
            i += 1
        if xnew == xvals[i]:
            ynew[0] = yvals[i]
            return
        if i == N:
            i = N-1
        x_a = xvals[i-1]
        y_a = yvals[i-1]
        x_b = xvals[i]
        y_b = yvals[i]
    slope = (xnew - x_a)/(x_b - x_a)
    ynew[0] = slope * (y_b-y_a) + y_a
    return

I believe you have a small mistake on line 31:

while i > 0 ...

should be

while i >= 0 ...

otherwise, when I do an interpolation of a sine wave I get this: