How to Optimize Grid Lines and Support Points for Piecewise Linear Interpolation?
I am working on a piecewise linear interpolation problem in Python where I optimize the placement of support points (in x- and y-direction) to fit some data. The 1D version of my code works as expected, distributing the support points nicely to fit the data, especially placing more support points in areas with smaller bending radii.
import numpy as np
from scipy.interpolate import interpn
from scipy.optimize import minimize
import matplotlib.pyplot as plt
# Erzeuge Testdaten: parabel um den nullpunkt
x_data = np.linspace(-3.15/2, 3.15/2, 800)
y_data = np.sin(x_data)
N = 8 # Anzahl der Stützstellen
def piecewise_linear_interp_error(support_points, x_data, y_data):
x_support = support_points[:N]
y_support = support_points[N:]
idx = np.argsort(x_support)
x_support = x_support[idx]
y_support = y_support[idx]
grid = (x_support,)
points = x_data[:, None]
# Interpolation der Stützstellen
y_interp = interpn(grid, y_support, points, method='linear', bounds_error=False, fill_value=None).flatten()
# Berechne den quadratischen Fehler
error = np.sqrt(np.mean((y_data - y_interp) ** 2))
return error
# Startwerte für die Optimierung (gleichmäßig verteilte Punkte)
initial_x_support = np.linspace(np.min(x_data), np.max(x_data), N)
initial_y_support = np.interp(initial_x_support, x_data, y_data)
initial_support_points = np.concatenate([initial_x_support, initial_y_support])
result = minimize(piecewise_linear_interp_error, initial_support_points, args=(x_data, y_data), method='SLSQP')
optimized_x_support = result.x[:N]
optimized_y_support = result.x[N:]
plt.plot(x_data, y_data, label='Original Data')
plt.plot(optimized_x_support, optimized_y_support, 'ro-', label='Optimized Support Points')
plt.legend()
plt.show()
print("Optimized Support Points:")
for x, y in zip(optimized_x_support, optimized_y_support):
print(f"({x:.2f}, {y:.2f})")
I modeled the 1-D and n-D versions to be similar. Most of the extra code is for piercing together, separating the optimization vector, and handling the grid for the multidimensional case.
import numpy as np
from scipy.interpolate import interpn
from scipy.optimize import minimize
import matplotlib.pyplot as plt
def display_support_points(src_grid, target_vec):
print("Optimized Support Points:")
for i, points in enumerate(src_grid):
formatted_points = ', '.join(f'{p:.2f}' for p in points)
print(f"({formatted_points}, {target_vec[i]})")
class PiecewiseLinearInterpolation:
def __init__(self, source_values, target_values, source_resolution):
self.n_dimensions = source_values.shape[1]
self.src_vals = source_values
self.target_values = target_values
self.src_resolution = source_resolution # List of support points for each dimension
self.src_vec_shape = None
self.initial_support_vector = self.generate_inital_support_vector()
def generate_inital_support_vector(self):
initial_src_support_vec = []
for i, x in enumerate(self.src_vals.T):
initial_src_support_vec.append(np.linspace(np.min(x), np.max(x), self.src_resolution[i]))
initial_src_support_vec = np.concatenate(initial_src_support_vec)
# Create a grid for each dimension based on the resolutions
src_grids = [np.linspace(np.min(x), np.max(x), res) for x, res in
zip(self.src_vals.T, self.src_resolution)]
# Create a meshgrid for interpolation
src_grid = np.array(np.meshgrid(*src_grids, indexing='ij')).T.reshape(-1, self.n_dimensions)
orig_source_idx = []
for dim in range(self.n_dimensions):
orig_source_idx.append(np.unique(self.src_vals[:, dim], ))
# reshape original target_values to the shape of src_vals
self.src_vec_shape = [len(np.unique(x)) for x in orig_source_idx]
target_vec = self.target_values.reshape(self.src_vec_shape)
initial_target_support = interpn(
orig_source_idx, # The grid for each dimension
target_vec, # The data to interpolate
src_grid, # The points where interpolation is needed
method='linear', # Interpolation method
bounds_error=False, # Do not raise an error if points are out of bounds
fill_value=None # Use NaN for out-of-bounds points
).flatten() # Flatten the result to get the initial support points for y
return np.concatenate([initial_src_support_vec, initial_target_support])
def calc_interpolation_error(self, support_vector):
src_vec, target_vec = self.split_support_vector(support_vector)
src_vec_sorted, target_vec_sorted = self.reorder_support_vectors(src_vec, target_vec)
points = np.array(self.src_vals).reshape(-1, self.n_dimensions)
target_interp = interpn(src_vec_sorted, target_vec_sorted, points, method='linear', bounds_error=False,
fill_value=None)
interp_err = np.sqrt(np.sum((self.target_values - target_interp) ** 2))
return interp_err
def reorder_support_vectors(self, src_vec, target_vec):
for i in range(self.n_dimensions):
idx = np.argsort(src_vec[i])
src_vec[i] = src_vec[i][idx]
target_vec = np.take(target_vec, idx, axis=i)
return src_vec, target_vec
def split_support_vector(self, support_vec):
src_vec = []
start = 0
for res in self.src_resolution:
end = start + res
src_vec.append(support_vec[start:end])
start = end
target_vec = support_vec[start:]
target_vec = target_vec.reshape(tuple(self.src_resolution))
return src_vec, target_vec
def optimize_support_points(self):
result = minimize(self.calc_interpolation_error, self.initial_support_vector, method='SLSQP')
result_src_vec, result_target_vec = self.split_support_vector(result.x)
res_src_vec_sorted, res_target_vec_sorted = self.reorder_support_vectors(result_src_vec, result_target_vec)
result_grid = np.array(np.meshgrid(*res_src_vec_sorted, indexing='ij')).T.reshape(-1, self.n_dimensions)
return result_grid, res_target_vec_sorted
def display_results(self, src_grid, target_vec):
plt.figure(figsize=(10, 8))
if src_grid.shape[1] == 1:
plt.plot(src_grid[:, 0], target_vec, 'ro-', label='Optimized Support Points')
elif src_grid.shape[1] == 2:
ax = plt.axes(projection='3d')
ax.scatter(self.src_vals[:, 0], self.src_vals[:, 1], self.target_values, label='Original Data')
ax.scatter(src_grid[:, 0], src_grid[:, 1], target_vec,
color='red', label='Optimized Support Points')
else:
raise ValueError("Plotting for dimensions higher than 2 is not supported.")
plt.legend()
plt.show()
# This is how its supposed to be used:
def main():
# Create a dataset with two dimensions for src_vals
x = np.linspace(-3.15 / 2, 3.15 / 2, 20)
y = np.linspace(-3.15 / 2, 3.15 / 2, 20)
X, Y = np.meshgrid(x, y)
Z = np.sin(X) + np.sin(Y)
raw_src_vec = np.array([X.flatten(), Y.flatten()]).T
raw_target_vec = Z.flatten()
x_resolutions = [5, 10] # Individual resolutions for each dimension
interpolator = PiecewiseLinearInterpolation(raw_src_vec, raw_target_vec, x_resolutions)
src_grid, target_vec = interpolator.optimize_support_points()
interpolator.display_results(src_grid, target_vec)
#display_support_points(src_grid, target_vec)
if __name__ == "__main__":
main()
To me, it seems as if the optimizer messed it up; however, it is more likely that I messed up with the multi-dimensional stuff. I wrote tests and verified pieces that seemed suspicious, and those were correct. How can I streamline the code to be more straightforward? Where is my data corruption?
There are two problems.
The first is that you're misunderstanding the meaning of the data (or using the wrong data to plot).
Take another look at the error function which you are trying to minimize.
def calc_interpolation_error(self, support_vector):
src_vec, target_vec = self.split_support_vector(support_vector)
src_vec_sorted, target_vec_sorted = self.reorder_support_vectors(src_vec, target_vec)
points = np.array(self.src_vals).reshape(-1, self.n_dimensions)
target_interp = interpn(src_vec_sorted, target_vec_sorted, points, method='linear', bounds_error=False,
fill_value=None)
interp_err = np.sqrt(np.sum((self.target_values - target_interp) ** 2))
return interp_err
This function calculates the error between the target values and the output of the interp function.
In other words, support_vector is just a parameter that allows the interp function to predict correctly, and there is no guarantee that it will be aligned with the Original Data.
The following plots the output of the interp function for randomly generated points.
def display_results(self, src_grid, target_vec):
plt.figure(figsize=(10, 8))
if src_grid.shape[1] == 1:
plt.plot(src_grid[:, 0], target_vec, "ro-", label="Optimized Support Points")
elif src_grid.shape[1] == 2:
ax = plt.axes(projection="3d")
ax.scatter(self.src_vals[:, 0], self.src_vals[:, 1], self.target_values, label="Original Data")
ax.scatter(src_grid[:, 0], src_grid[:, 1], target_vec.ravel(),
color="red", label="Optimized Support Points")
# Generate random points to predict.
x_min, x_max = self.src_vals[:, 0].min(), self.src_vals[:, 0].max()
y_min, y_max = self.src_vals[:, 1].min(), self.src_vals[:, 1].max()
n = 1000
points = np.random.default_rng(0).random((n, 2)) * ([x_max - x_min, y_max - y_min]) + [x_min, y_min]
# Predict (interpolate) the data for the random points.
src_vec = src_grid[:self.src_resolution[0], 0], src_grid[::self.src_resolution[0], 1]
predicted_data = interpn(src_vec, target_vec, points, method="linear", bounds_error=False,
fill_value=None)
# Plot the predicted data.
ax.scatter(points[:, 0], points[:, 1], predicted_data.ravel(),
color="orange", label="Predicted Data")
else:
raise ValueError("Plotting for dimensions higher than 2 is not supported.")
plt.legend()
plt.show()
Here is the result. The orange dots are predicted data using support points. As I will explain later, please ignore the support points for now.
As you can see, the results are fairly good.
The second problem is much simple. There is a bug in the plot.
target_vec should be transposed.
ax.scatter(src_grid[:, 0], src_grid[:, 1], target_vec.T.ravel(),
color='red', label='Optimized Support Points')
I believe this is (roughly) what you were expecting.
Note that if you analyze the result carefully, you will find that half of the support points are not being used because they are outside the area. This may be solved by limiting the range of the support points or by allowing the error function to evaluate points outside the source grid. But this is something you will have to experiment with, so I'll leave it at that.
As a side note, using the nearest interpolation instead of linear made it easier to understand what is happening, and easier to debug.
This may also help your further development :)