ENH: Add basic interpolation routine

[mmcky]

Add basic jit enabled interpolation routine to quantecon and leave more complex interpolation to the interpolations package.

[mmcky]

@jit
def interp(grid, vals, x):
    """
    Linearly interpolate (grid, vals) to evaluate at x.
    Here grid must be regular (evenly spaced).

    Parameters
    ----------
    grid and vals are numpy arrays, x is a float

    Returns
    -------
    a float, the interpolated value
    """

    a, b, len_g = np.min(grid), np.max(grid), len(grid)
    s = (x - a) / (b - a)
    q_0 = max(min(int(s * (len_g - 1)), (len_g - 2)), 0)
    v_0 = vals[q_0]
    v_1 = vals[q_0 + 1]
    λ = s * (len_g - 1) - q_0
    return (1 - λ) * v_0 + λ * v_1

[jstac]

Thanks @mmcky . It would be really helpful for me if this was added and released soon.

At the same time, it might be best if we change it everwhere so that the order of the positional arguments matches https://numpy.org/doc/stable/reference/generated/numpy.interp.html

That means searching through all lectures to find every place where it's called.

In terms of adding tests, the function should match the output of numpy.interp

[mmcky]

thanks @jstac.

Will prioritise this.

[mmcky]

Thanks @mmcky . It would be really helpful for me if this was added and released soon.

At the same time, it might be best if we change it everwhere so that the order of the positional arguments matches https://numpy.org/doc/stable/reference/generated/numpy.interp.html

That means searching through all lectures to find every place where it's called.

In terms of adding tests, the function should match the output of numpy.interp

@jstac looking at the numpy docs it seems the function splits the x and y coordinates

https://numpy.org/doc/stable/reference/generated/numpy.interp.html

Are you proposing we change from

def interp(grid, vals, x):

to

def interp(x, grid, vals):

and grid is equivalent to numpy xp and vals are fp. Is that correct?

[jstac]

That's right @mmcky .

We might even change grid to xp and vals to fp throughout.

[mmcky]

@jstac

q_0 = max(min(int(s * (len_g - 1)), (len_g - 2)), 0)

would you like this function to broadcast over an array of x values?

numba doesn't know how to cast int over an array.

or is x always just a single float value in this case for use inside another jit routine

I am running it over this example

x = np.linspace(0, 2*np.pi, 10)
y = np.sin(x)
xvals = np.linspace(0, 2*np.pi, 50)
numpy_y = np.interp(xvals, x, y)
qe_y = interp(xvals, x, y)

and it has the following issue

No implementation of function Function(<class 'int'>) found for signature:

 >>> int(array(float64, 1d, C))

[jstac]

Ah, good catch. It might be better if it broadcasts --- and also returns a float if called on a float.

However, if that's too much hassle and the calling code in the lectures only calls on floats I'm fine without broadcasting.

[mmcky]

Closing this as using np.interp