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?
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]