Skip to content

Instantly share code, notes, and snippets.

What would you like to do?
# 2d comparison
import numpy as np
x = np.linspace(-1,1,100)
y = np.linspace(-1,1,100)
f = lambda x,y: np.sinc(x**2+y**2)
vals = np.array( [[f(xx,yy) for yy in y] for xx in x] )
N = 1000000 # number of points to evaluate
eval_points = -1+2*np.random.rand(N*2).reshape((N,2))
# irregular grid (linear)
print(f"# Interpolating a {len(x)}×{len(y)} grid at {N} points.")
print("## Irregular Grid")
import scipy
import scipy.interpolate
bspl = scipy.interpolate.RectBivariateSpline(x,y,vals, kx=1,ky=1)
out1 = bspl.ev(eval_points[:,0], eval_points[:,1])
%time out1 = bspl.ev(eval_points[:,0], eval_points[:,1])
from HARK import interpolation as eai
bli = eai.BilinearInterp(vals, x, y)
out2 = bli(eval_points[:,0], eval_points[:,1])
%time out2 = bli(eval_points[:,0], eval_points[:,1])
# regular spacing (linear)
from interpolation import interp
out3 = interp(x,y,vals,eval_points)
%time out3 = interp(x,y,vals,eval_points)
# new API
from interpolation import mlinterp
grid_uneven = (x,y) # tuple with two vectorw
out4 = mlinterp(grid_even,vals,eval_points)
%time out4 = mlinterp(grid_uneven,vals,eval_points)
print("\n## Regular Grid")
grid_even = (
out5 = mlinterp(grid_even,vals,eval_points)
%time out5 = mlinterp(grid_even,vals,eval_points)
# old API (optimized code generation)
from interpolation.splines import LinearSpline
ls = LinearSpline([-1.0,-1.0],[1.0,1.0], [100,100], vals)
out6 = ls(eval_points)
%time out6 = ls(eval_points)
l = [out1,out2,out3,out4,out5,out6]
print( abs(sum(l)-out5*len(l)).max() )
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment