Skip to content

Instantly share code, notes, and snippets.

@teoliphant
Last active April 16, 2023 09:30
Show Gist options
  • Save teoliphant/2edb542689af5f4732e83f54cb121d4c to your computer and use it in GitHub Desktop.
Save teoliphant/2edb542689af5f4732e83f54cb121d4c to your computer and use it in GitHub Desktop.
Vectorized linear interpolation
# 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
@dvincentwest
Copy link

dvincentwest commented Aug 30, 2018

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:

image

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment