supplement to stackoverflow comment on this post: https://stackoverflow.com/questions/75427538/regulargridinterpolator-excruciatingly-slower-than-interp2d

manual_RegularGridInterpolator.py

# supplement to my stackoverflow comment on this post: https://stackoverflow.com/questions/75427538/regulargridinterpolator-excruciatingly-slower-than-interp2d

# since `x_vals` and `y_vals` are decreasing rather than increasing, we do some trickery on top of `np.searchsorted()` to get the desired indices.
i_x = x_vals.shape[0] - np.searchsorted(np.flip(x_vals), x)
i_y = y_vals.shape[0] - np.searchsorted(np.flip(y_vals), y)

'''here's an option with list comprehension
points = [
    (
        x_vals[i_x-1+i],
        y_vals[i_y-1+j],
        z_grid[i_x-1+i, i_y-1+j]
    )
    for i, j in [(i, j) for i in range(2) for j in range(2)]
]
'''

points = (
    (
        x_vals[i_x - 1],
        y_vals[i_y - 1],
        z_grid[i_x - 1, i_y - 1]
    ),
    (
        x_vals[i_x],
        y_vals[i_y - 1],
        z_grid[i_x, i_y - 1]
    ),
    (
        x_vals[i_x - 1],
        y_vals[i_y],
        z_grid[i_x - 1, i_y]
    ),
    (
        x_vals[i_x],
        y_vals[i_y],
        z_grid[i_x, i_y]
    )
)

return bilinear_interpolation(x, y, points)