Smooth linear interpolation using NumPy

Question

I need to recode an equivalent of the NumPy interp function. In order to do that, I thought of building a pairwise distance matrix, and then using this matrix to define "weights" that I use to multiply the target values.

def interp_diff(x, xp, fp):
    dxp = np.diff(xp)
    dxp = np.concatenate(([dxp[0]], dxp))

    pairwise_distance = np.abs(x[:, np.newaxis] - xp)  # Calculate the absolute differences between x and xp

    # Add a small epsilon to pairwise distances to avoid division by zero
    epsilon = 1e-8
    pairwise_distance += epsilon

    weights = 1 / pairwise_distance  # Calculate the weights based on the differences

    # Normalize the weights
    weights /= np.sum(weights, axis=1)[:, np.newaxis]

    # Apply hardness threshold to the weights
    weights = np.where(weights > 0, weights, 0)

    return np.dot(weights, fp)

This produces expected results when xp values are placed on a regular grid, but it does not work if xp values are not evenly spaced. How can I make this work for any spacing of xp?

The constraint I have is that I can't use index related methods (argsort, argwhere, searchsorted, etc...). Which is what makes it a bit challenging.

Example usage on a regular grid:

x_np = np.linspace(0, 5, 4)
y_np = np.sin(x_np)
x_i = x_np
x_i = np.linspace(0, 5, 10)
y_i = interp_diff(x_i, xp=x_np, fp=y_np)
ax = plt.figure().gca()
ax.scatter(x_np, y_np)
ax.scatter(x_i, y_i, marker='+', s=30, alpha=0.7)
ax.plot(x_i, y_i)
plt.show()

Produces the expected result:

However, switching to non evenly spaced grid:

def sinspace(start, stop, num):

    ones = 0 * start + 1
    return start + (stop - start) * (1 - np.cos(np.linspace(
        0 * ones,
        np.pi / 2 * ones,
        num
    )))
x_np = sinspace(0, 5, 4)
y_np = np.sin(x_np)
x_i = x_np
x_i = np.linspace(0, 5, 10)
y_i = interp_diff(x_i, xp=x_np, fp=y_np)
ax = plt.figure().gca()
ax.scatter(x_np, y_np)
ax.scatter(x_i, y_i, marker='+', s=30, alpha=0.7)
ax.plot(x_i, y_i)
plt.show()

This is what I get:

Answer 1

I'm unable to replicate your results with a slightly modified version of the code you supplied:

import numpy
from matplotlib import pyplot

def interp_diff(x, xp, fp):
    dxp = numpy.diff(xp)
    dxp = numpy.concatenate(([dxp[0]], dxp))

    pairwise_distance = numpy.abs(x[:, numpy.newaxis] - xp)  # Calculate the absolute differences between x and xp

    # Add a small epsilon to pairwise distances to avoid division by zero
    epsilon = 1e-8
    pairwise_distance += epsilon

    weights = 1 / pairwise_distance  # Calculate the weights based on the differences

    # Normalize the weights
    weights /= numpy.sum(weights, axis=1)[:, numpy.newaxis]

    # Apply hardness threshold to the weights
    weights = numpy.where(weights > 0, weights, 0)

    return numpy.dot(weights, fp)

# 1 : 1 grid
x_np = numpy.linspace(0,5,4)
y_np = numpy.sin(x_np)
x_i = x_np
y_i = interp_diff(x_i,xp=x_np,fp=y_np)
ax = pyplot.figure().gca()
ax.scatter(x_np,y_np)
ax.scatter(x_i,y_i,marker='+',s=30,alpha=0.7)
ax.plot(x_i,y_i)
pyplot.show()

# 1 : 2.5 grid
x_np = numpy.linspace(0,5,4)
y_np = numpy.sin(x_np)
x_i = numpy.linspace(0,5,10)
y_i = interp_diff(x_i,xp=x_np,fp=y_np)
ax = pyplot.figure().gca()
ax.scatter(x_np,y_np)
ax.scatter(x_i,y_i,marker='+',s=30,alpha=0.7)
ax.plot(x_i,y_i)
pyplot.show()

Which version of the numpy and matplotlib libraries are you using?

Answer 2

The interp_diff function to work with unevenly spaced xp values, i have chnaged with np.subtract.outer function to calculate the pairwise differences between x and xp without relying on index-related methods.

import numpy as np

import matplotlib.pyplot as plt

def interp_diff(x, xp, fp):
    dxp = np.diff(xp)
    dxp = np.concatenate(([dxp[0]], dxp))

pairwise_distance = np.abs(np.subtract.outer(x, xp))  # Calculate the absolute differences between x and xp

# Add a small epsilon to pairwise distances to avoid division by zero
epsilon = 1e-8
pairwise_distance += epsilon

weights = 1 / pairwise_distance  # Calculate the weights based on the differences

# Normalize weights
weights /= np.sum(weights, axis=1)[:, np.newaxis]

# Apply hardness threshold to the weights
weights = np.where(weights > 0, weights, 0)

return np.dot(weights, fp)

x_np = np.linspace(0, 5, 4)
y_np = np.sin(x_np)
x_i = np.linspace(0, 5, 10)
y_i = interp_diff(x_i, xp=x_np, fp=y_np)

ax = plt.figure().gca()
ax.scatter(x_np, y_np)
ax.scatter(x_i, y_i, marker='+', s=30, alpha=0.7)
ax.plot(x_i, y_i)
plt.show()

Answer 3

Not sure this will do exactly what you want, since it uses a few boolean arrays which may cause problems with your gradients. But it produces the correct result (comparing to Numpy) and doesn't use any complicated indexing:

def interp_diff(x, xp, fp):
    diff_filter = np.diff(x[:, np.newaxis] >= xp, axis=-1)
    prev_filter = np.concatenate((diff_filter, ~np.any(
        diff_filter, axis=-1, keepdims=True)), axis=-1)
    next_filter = np.concatenate(
        (np.broadcast_to(False, (diff_filter.shape[0], 1)), diff_filter), axis=-1)

    epsilon = 1e-8
    pairwise_distance = xp - x[:, np.newaxis]
    prev_dist = np.where(prev_filter, -pairwise_distance, 0)
    next_dist = np.where(next_filter, pairwise_distance, 0)
    progress = prev_dist / \
        ((prev_dist + next_dist).sum(axis=-1) + epsilon)[..., np.newaxis]

    return np.sum(np.where(prev_filter, fp, 0) + np.concatenate((np.diff(fp), [0])) * progress, axis=-1)

It works by abusing np.diff on boolean arrays to get the previous and the subsequent point in xp. It then calculates how far between the two points the new point is and uses np.diff again to offset the new point from the previous closest point.