Feature request: Support for interpolation

[moorepants]

Discussed in #473

^{Originally posted by Peter230655 May 25, 2025}
It seems that symbolic interpolation becomes time consuming, if more than a few points are to be interpolated.

[moorepants]

One way to support this would be to create a symbolic unevaluated function that represents interpolation, e.g.:

f = Interpolation(x(t), xdata, ydata)

This function would need to support being differentiated, for example: f.diff(x(t)). For opty's purpose, the function should be continuous.

Both f and f.diff(x(t)) should support being passed to ccode() and pycode().

The basic ideas for custom functions are here: https://docs.sympy.org/latest/guides/custom-functions.html

If the NumPyCodePrinter is used on f, for example, it should generate something like:

>>> NumPyCodePrinter.doprint(f)
'numpy.interp(x, xdata, ydata)'

So lambdify((x, xdata, ydata), f) would generate a function that just calls NumPy's interp() function.

The C code printer would need to output a manually written interpolation routine.

For any function that SymPy does not already have, this is the general process for adding support.

[moorepants]

Also, this can be implemented manually by a user, so I doubt I would add it to opty. If it is added to any software, it should be SymPy.

[moorepants]

An interpolator needs to know that the data are 1D arrays, so this is probably needed:

x_vector = MatrixSymbol('x', n, 1)
y_vector = MatrixSymbol('y', n, 1)
x = Symbol('x')
y = Interpolate(x, x_vector, y_vector)

[Peter230655]

def smooth_hump(x, a, b, steepness):
    """returns approx 1 for a < x < b, 0 otherwise
    except returns 0.5 for x = a or x = b"""
    return 0.5 * (sm.tanh(steepness * (x - a)) - sm.tanh(steepness * (x - b)))

amp = 30.0
omega = 2*np.pi/500.0  # one period every 500 meters
xp = np.linspace(-250.0, 1250.0, num=1499)
yp = amp*np.sin(omega*xp) 
thetap = np.atan(amp*omega*np.cos(omega*xp))
steepness = sm.symbols('steepness', real=True)

def smooth_theta(x, steepness):
    """smooth piecewise linear approximation. It is not good if x = xp[i] for
    some i"""
    result = sum([ smooth_hump(x, xp[i], xp[i + 1], steepness) * (
        ((thetap[i + 1] - thetap[i]) / (xp[i + 1] - xp[i]) * (x - xp[i]) + thetap[i]))
                      for i in range(len(xp) - 1)])
    return result

I tried this 'brute force' approach at linear interpolation, above.
I use a hump function to select which intervall x is in.
If x is far away from a point of support, it is quite accurate, if x = xp[i], some i is is very bad.
As a sum of 1499 functions of tanh(..) its derivative must be huge
I is not practical:

• with amp = 10, it took about 17,500 iterations,
• with amp = 30, 40,000 iterations still were not sufficient,

Update:
After I increased the propulsion power p(t), it only needed around 100 iterations, even for amp = 30.
It solves fast, 0.85 sec, but it takes very long to build Problem.
When I tried backend='numpy' I got this error:
RecursionError: maximum recursion depth exceeded during compilation

Reducing len(xp) from 1499 to 373 reduced the time to build Problem to around 30 sec, solution time around 0.5 sec.
backend='numpy' does not throw an error.

Update
With increased p(t) also @moorepants original version converged rapidly, around 120 iterations with amp = 30 and time for solving only 0.45 sec.

[tvdbogert]

As discussed in #472, we may just want to do this numerically.

sm.interpolate creates a single polynomial that passes exactly through all data points. With more data points, this becomes a very high order polynomial which goes wild near the endpoints of the x range. I would not consider this an acceptable method for interpolation.

The spline interpolation seems more useful. The Piecewise expression is guaranteed to be continuous. And even guaranteed to be smooth when cubic spline is selected. There is just the matter of what it will return when x is out of range.

[Peter230655]

You feel my smooth_theta(...) is not such a brilliant idea for linear interpolation? :-)
Up to about 375 data points it seems to work reasonably fast.
In Jason's plot_hilly_ride.py I found no difference in convergence speed, number of iterations needed and accuracy between

• Jason's original known_trajectory_map idea
• my smooth_theata with 373 data points
• using the symbolic function, this sm.atan(.....)

Of course a very simple example.

Can spline interpolation handle many data points?
I did not try.

[Peter230655]

I wrote this snippet just to see what my differentiable hump linear interpolation would look like.
I think, it behaves as expected. It seems to gets the y coordinates of points which should be on the interpolating graph quite well: the red dots are on the blue lines.
Its derivative looks like expected, I think, since the derivative of the linear interpolation is discontinuous at certain points, the derivative of differentiable hump linear interpolation changes rapidly, but continuously from one constant value to the next.

For x_values < 0.0 there are problems, I hope I will get them ironed out.

I have no idea whether this is a practical approach to linear interpolation. At least for up to about 375 points it works reasonably fast. If practical, it can become part of the equations of motion, as it is symbolic.

import numpy as np
import sympy as sm
import matplotlib.pyplot as plt


def smooth_hump(x, a, b, steepness):
    """returns approx 1 for a < x < b, 0 otherwise
    except returns 0.5 for x = a or x = b"""
    return 0.5 * (sm.tanh(steepness * (x - a)) - sm.tanh(steepness * (x - b)))

# The coordinates of the points, between which the linear interpolation is done.
x_values = [0.0, 2.0, 4.0, 6.0, 8.0, 10.0]
y_values = [0.0, 2.0, -3.0, 5.0, 8.0, 0.0]


def smooth_lin_interp(x, steepness):
    """smooth piecewise linear approximation."""
    result = sum([ smooth_hump(x, x_values[i], x_values[i + 1], steepness) * (
        ((y_values[i + 1] - y_values[i]) / (x_values[i + 1] - x_values[i]) * 
         (x - x_values[i]) + y_values[i]))
                      for i in range(len(x_values) - 1)])
    return result


X, Steepness = sm.symbols('X, Steepness', real=True)
lin_interp_lam = sm.lambdify([X, Steepness],
                        smooth_lin_interp(X, Steepness), cse=True)

lin_interp_dx = smooth_lin_interp(X, Steepness).diff(X)
lin_interp_dx_lam = sm.lambdify([X, Steepness], lin_interp_dx, cse=True)

Steepness = 10.0
fig, ax = plt.subplots(figsize=(8, 4))
ax.plot(x_values, y_values)
for point in (0.0, 2.5, 3.7, 5, 8, 9.5, 10):
    ax.plot(point, lin_interp_lam(point, Steepness), 'ro')
    
DT = []
for point in np.linspace(0.0, 10, 100):
    DT.append(lin_interp_dx_lam(point, Steepness))  
ax.plot(np.linspace(0.0, 10, 100), DT)
_ = ax.set_title((f'The blue line is the linear interpolation between a set of'
              f' given points. \n The red dots are points of the smooth theta'
              f' function. \n The orange line is the derivative of the smooth' 
              f'  theta function'))

[Peter230655]

Seems the problem is not with x_values < 0..0, but with x_values very close to the first and the last value of the interval.

[Peter230655]

I now compared sympy.interpolating_spline(1, x, x_values, y_values) to my hump interpolation.

• Interpolationg_spline has no issues with edge points
• hump interpolation seems a bit faster (I compared essentially the time to lambdify the functions):
  • len(x_values) = 100: 0.5 sec vs. 2.2 sec on my PC
  • len(x_values) = 600: 3.4 sec vs. 57 sec
  • len(x_values) = 1501: 12.4 sec vs. 331 sec. (This was the number used in plot_hilly_ride

So, unless len(x_values) is large, definitely sympy.interpolation_spline is better!

NB:
Being retired, I never deal with real world problems.
@moorepants , @tvdbogert : What would be practical sizes of len(x_values) ? Thanks!

[moorepants]

One of our walking datasets has 50,000 time samples for a single subject.

[Peter230655]

Based on my limited experience neither symbolic interpolation function would finish in our life time.

[moorepants]

For sure, we already know that. When I've said "symbolic interpolator" before, I haven't intended its meaning in the literal sense, but in the sense that we need a symbolic object that represents interpolation (unevaluated), as described above.

[Peter230655]

As discussed in #472, we may just want to do this numerically.

sm.interpolate creates a single polynomial that passes exactly through all data points. With more data points, this becomes a very high order polynomial which goes wild near the endpoints of the x range. I would not consider this an acceptable method for interpolation.

The spline interpolation seems more useful. The Piecewise expression is guaranteed to be continuous. And even guaranteed to be smooth when cubic spline is selected. There is just the matter of what it will return when x is out of range.

Just info:
When I tried sympy.interpolate_spline(...) with opty, I got this error:
ValueError: All Piecewise expressions must contain an (expr, True) statement to be used as a default condition. Without one, the generated expression may not evaluate to anything under some condition.

(This is only the last line, it was a very long error message)