How to translate IDL's SPLINE function to Python [particularly for the case we have 3 data points]

Question

The SPLINE function in IDL allows for cubic interpolation on data (with at at least 3 data points). While the Scipy library in Python can perfomr similar computations with the UnivariateSpline and the splrep functions, these break if the interpolations are set to cubic ones and we have just 3 data points (something that doesn't happen with SPLINE).

This is a simple example of what SPLINE does in IDL:

> x = [2., 3., 4.]
> y = (x-3)^2
> t = FINDGEN(20)/10.+2
> z = SPLINE(x, y, t)
> print, z

1.00000     0.810662     0.642087     0.493590     0.364684     0.255081     0.164684    0.0935898    0.0420876    0.0106618      0.00000    0.0106618    0.0420876    0.0935898     0.164684     0.255081     0.364684     0.493590      0.642087     0.810662

But if I try to do this in Python with from scipy.interpolate import splrep, splev

x = np.array([2., 3., 4.])
y = (x-3)**2
t = np.arange(20)/10.+2
z = splev(t, splrep(x, y, k = 3))

or with from scipy.interpolate import UnivariateSpline

x = np.array([2., 3., 4.])
y = (x-3)**2
t = np.arange(20)/10.+2
z = UnivariateSpline(x, y, k = 3)(t)

I always get this error message:

TypeError: m > k must hold

Which I understand since there can't be a unique k-degree polynomial solution when we have to fit m data points if m ≤ k. But then it begs the question... How does SPLINE in IDL performs this calculation? And how can I reporduce it in Python?

I can try by lowering the polynomical to k = 2 (a quadratic interpolation), like so

z = splev(t, splrep(x, y, k = 2))

or

z = UnivariateSpline(x, y, k = 2)(t)

And I will get in both cases:

> print(z)
[1. 0.81 0.64 0.49 0.36 0.25 0.16 0.09 0.04 0.01 0. 0.01 0.04 0.09 0.16 0.25 0.36 0.49 0.64 0.81]

Which is certainly similar to the output in IDL unless we ignore everything below the second decimal place.

How can I perform the same calculations as SPLINE in Python, even in the case that k = m like SPLINE does?

Answer 1

Since this remained unanswered, I've translated the 170 lines of code of SPLINE.pro to Python3, as so:

import numpy as np

def spline(X, Y, T, sigma = 1.0):
    n = min(len(X), len(Y))
    if n <= 2:
        print('X and Y must be arrays of 3 or more elements.')
    if sigma != 1.0:
        sigma = min(sigma, 0.001)
    yp = np.zeros(2*n)
    delx1 = X[1]-X[0]
    dx1 = (Y[1]-Y[0])/delx1
    nm1 = n-1
    nmp = n+1
    delx2 = X[2]-X[1]
    delx12 = X[2]-X[0]
    c1 = -(delx12+delx1)/(delx12*delx1)
    c2 = delx12/(delx1*delx2)
    c3 = -delx1/(delx12*delx2)
    slpp1 = c1*Y[0]+c2*Y[1]+c3*Y[2]
    deln = X[nm1]-X[nm1-1]
    delnm1 = X[nm1-1]-x[nm1-2]
    delnn = X[nm1]-X[nm1-2]
    c1 = (delnn+deln)/(delnn*deln)
    c2 = -delnn/(deln*delnm1)
    c3 = deln/(delnn*delnm1)
    slppn = c3*Y[nm1-2]+c2*Y[nm1-1]+c1*Y[nm1]
    sigmap = sigma*nm1/(X[nm1]-X[0])
    dels = sigmap*delx1
    exps = np.exp(dels)
    sinhs = 0.5*(exps-1/exps)
    sinhin = 1/(delx1*sinhs)
    diag1 = sinhin*(dels*0.5*(exps+1/exps)-sinhs)
    diagin = 1/diag1
    yp[0] = diagin*(dx1-slpp1)
    spdiag = sinhin*(sinhs-dels)
    yp[n] = diagin*spdiag
    delx2 = X[1:]-X[:-1]
    dx2 = (Y[1:]-Y[:-1])/delx2
    dels = sigmap*delx2
    exps = np.exp(dels)
    sinhs = 0.5*(exps-1/exps)
    sinhin = 1/(delx2*sinhs)
    diag2 = sinhin*(dels*(0.5*(exps+1/exps))-sinhs)
    diag2 = np.concatenate([np.array([0]), diag2[:-1]+diag2[1:]])
    dx2nm1 = dx2[nm1-1]
    dx2 = np.concatenate([np.array([0]), dx2[1:]-dx2[:-1]])
    spdiag = sinhin*(sinhs-dels)
    for i in range(1, nm1):
        diagin = 1/(diag2[i]-spdiag[i-1]*yp[i+n-1])
        yp[i] = diagin*(dx2[i]-spdiag[i-1]*yp[i-1])
        yp[i+n] = diagin*spdiag[i]
    diagin = 1/(diag1-spdiag[nm1-1]*yp[n+nm1-1])
    yp[nm1] = diagin*(slppn-dx2nm1-spdiag[nm1-1]*yp[nm1-1])
    for i in range(n-2, -1, -1):
        yp[i] = yp[i]-yp[i+n]*yp[i+1]
    m = len(T)
    subs = np.repeat(nm1, m)
    s = X[nm1]-X[0]
    sigmap = sigma*nm1/s
    j = 0
    for i in range(1, nm1+1):
        while T[j] < X[i]:
            subs[j] = i
            j += 1
            if j == m:
                break
        if j == m:
            break
    subs1 = subs-1
    del1 = T-X[subs1]
    del2 = X[subs]-T
    dels = X[subs]-X[subs1]
    exps1 = np.exp(sigmap*del1)
    sinhd1 = 0.5*(exps1-1/exps1)
    exps = np.exp(sigmap*del2)
    sinhd2 = 0.5*(exps-1/exps)
    exps = exps1*exps
    sinhs = 0.5*(exps-1/exps)
    spl = (yp[subs]*sinhd1+yp[subs1]*sinhd2)/sinhs+((Y[subs]-yp[subs])*del1+(Y[subs1]-yp[subs1])*del2)/dels
    if m == 1:
        return spl[0]
    else:
        return spl

Now everything works as expected and the result are the exact same as in IDL.

x = np.array([2., 3., 4.])
y = (x-3)**2
t = np.arange(20)/10.+2
z = spline(x, y, t)

> print(z)
[1.         0.81066204 0.64208752 0.49359014 0.36468451 0.25508134
 0.16468451 0.09359014 0.04208752 0.01066204 0.         0.01066204
 0.04208752 0.09359014 0.16468451 0.25508134 0.36468451 0.49359014
 0.64208752 0.81066204]