如何将 IDL 的 SPLINE 函数转换为 Python [特别是对于我们有 3 个数据点的情况]

Question

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).

这是 SPLINE 在 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

但是如果我尝试在 Python 中使用来自 scipy.interpolate 导入 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))

或与从 scipy.interpolate 导入 UnivariateSpline

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

我总是收到此错误消息：

TypeError: m > k must hold

我明白这一点，因为如果 m ≤ k，当我们必须拟合 m 个数据点时，不可能有唯一的 k 次多项式解。但随之而来的问题是……IDL 中的 SPLINE 如何执行此计算？我怎样才能在 Python 中重新制作它？

我可以尝试将多项式降低到 k = 2（二次插值），就像这样

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

或

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

在这两种情况下我都会得到：

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

这肯定与 IDL 中的输出类似，除非我们忽略小数点后第二位以下的所有内容。

如何在 Python 中执行与 SPLINE 相同的计算，即使 k = m 就像 SPLINE 一样？

Answer 1

由于此问题仍未得到解答，我已将 SPLINE.pro 的 170 行代码翻译为 Python3，如下所示：

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

现在一切都按预期工作，结果与 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]

如何将 IDL 的 SPLINE 函数转换为 Python [特别是对于我们有 3 个数据点的情况]

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

python

scipy

interpolation

idl-programming-language