interpolation_performance_comparison

# 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()
print("## Irregular Grid")
print("\nscipy.interpolate.RectBivariateSpline")
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])

print("\nHARK.interpolate.BilinearInterp")
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)
print("\ninterpolation.interp")
from interpolation import interp
out3 = interp(x,y,vals,eval_points)
%time out3 = interp(x,y,vals,eval_points)



# new API

print("\ninterpolation.mlinterp")

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")

print("\ninterpolation.mlinterp")
grid_even = (
    (-1.0,1.0,100),
    (-1.0,1.0,100)
)
out5 = mlinterp(grid_even,vals,eval_points)
%time out5 = mlinterp(grid_even,vals,eval_points)


# old API (optimized code generation)
print("\ninterpolation.LinearSpline")
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() )