interp_new.py

@register_jitable
def np_interp_impl_inner(x, xp, fp, dtype):
    x_arr = np.asarray(x)
    xp_arr = np.asarray(xp)
    fp_arr = np.asarray(fp)

    if len(xp_arr) == 0:
        raise ValueError('array of sample points is empty')

    if len(xp_arr) != len(fp_arr):
        raise ValueError('fp and xp are not of the same size.')

    if xp_arr.size == 1:
        return np.full(x_arr.shape, fill_value=fp_arr[0], dtype=dtype)

    if not np.all(xp_arr[1:] > xp_arr[:-1]):
        msg = 'xp must be monotonically increasing'
        raise ValueError(msg)
        # note: NumPy docs suggest this is required but it is not
        # checked for or enforced; see:
        # https://github.com/numpy/numpy/issues/10448
        # This check is quite expensive.

    out = np.empty(x_arr.shape, dtype=dtype)
    idx = 0
    last_idx = np.nan

    for i in range(x_arr.size):
        # shortcut if possible
        if np.isnan(x_arr.flat[i]):
            out.flat[i] = np.nan
            continue
        if x_arr.flat[i] >= xp_arr[-1]:
            out.flat[i] = fp_arr[-1]
            continue
        if x_arr.flat[i] <= xp_arr[0]:
            out.flat[i] = fp_arr[0]
            continue

        # otherwise, avoid binary search if we're in the same
        # region or an immediately neighbouring region
        if xp_arr[idx - 1] < x_arr.flat[i] <= xp_arr[idx]:
            pass
        elif xp_arr[idx] < x_arr.flat[i] <= xp_arr[idx + 1]:
            idx += 1
        elif xp_arr[idx - 2] < x_arr.flat[i] <= xp_arr[idx - 1]:
            idx -= 1
        else:
            idx = np.searchsorted(xp_arr, x_arr.flat[i])

        if x_arr.flat[i] == xp_arr[idx]:
            # replicate numpy behaviour which is present
            # up to (and including) 1.15.4, but fixed in
            # this PR: https://github.com/numpy/numpy/pull/11440
            numerator = fp_arr[idx + 1] - fp_arr[idx]
            denominator = xp_arr[idx + 1] - xp_arr[idx]
            next_slope = numerator / denominator
            if not np.isfinite(next_slope):
                out.flat[i] = np.nan
            else:
                out.flat[i] = fp_arr[idx]
        else:
            if idx != last_idx:
                numerator = fp_arr[idx] - fp_arr[idx - 1]
                denominator = xp_arr[idx] - xp_arr[idx - 1]
                slope = numerator / denominator

            delta_x = x_arr.flat[i] - xp_arr[idx - 1]
            out.flat[i] = fp_arr[idx - 1] + slope * delta_x

        last_idx = idx
    return out