python - 边界处有约束的样条

Question

我在三维网格上测量了数据，例如f(x, y, t)。我想用样条在 t 的方向上插入和平滑这些数据。目前，我使用 scipy.interpolate.UnivariateSpline 执行此操作:

import numpy as np
from scipy.interpolate import UnivariateSpline

# data is my measured data
# data.shape is (len(y), len(x), len(t))
data = np.arange(1000).reshape((5, 5, 40))  # just for demonstration
times = np.arange(data.shape[-1])
y = 3
x = 3
sp = UnivariateSpline(times, data[y, x], k=3, s=6)

但是，我需要样条曲线在 t=0 处具有消失的导数。有没有办法强制执行此约束？

Answer 1

我能想到的最好的办法是使用 scipy.optimize.minimize 进行约束最小化。样条的导数很容易，所以约束很简单。我会使用常规样条拟合 (UnivariateSpline) 来获得结点 (t)，并保持结点固定(度数 k，当然)，并改变系数 c。也许还有一种方法可以改变结的位置，但我会把它留给你。

import numpy as np
from scipy.interpolate import UnivariateSpline, splev, splrep
from scipy.optimize import minimize

def guess(x, y, k, s, w=None):
    """Do an ordinary spline fit to provide knots"""
    return splrep(x, y, w, k=k, s=s)

def err(c, x, y, t, k, w=None):
    """The error function to minimize"""
    diff = y - splev(x, (t, c, k))
    if w is None:
        diff = np.einsum('...i,...i', diff, diff)
    else:
        diff = np.dot(diff*diff, w)
    return np.abs(diff)

def spline_neumann(x, y, k=3, s=0, w=None):
    t, c0, k = guess(x, y, k, s, w=w)
    x0 = x[0] # point at which zero slope is required
    con = {'type': 'eq',
           'fun': lambda c: splev(x0, (t, c, k), der=1),
           #'jac': lambda c: splev(x0, (t, c, k), der=2) # doesn't help, dunno why
           }
    opt = minimize(err, c0, (x, y, t, k, w), constraints=con)
    copt = opt.x
    return UnivariateSpline._from_tck((t, copt, k))

然后我们生成一些应该具有零初始斜率的假数据并对其进行测试:

import matplotlib.pyplot as plt

n = 10
x = np.linspace(0, 2*np.pi, n)
y0 = np.cos(x) # zero initial slope
std = 0.5
noise = np.random.normal(0, std, len(x))
y = y0 + noise
k = 3

sp0 = UnivariateSpline(x, y, k=k, s=n*std)
sp = spline_neumann(x, y, k, s=n*std)

plt.figure()
X = np.linspace(x.min(), x.max(), len(x)*10)
plt.plot(X, sp0(X), '-r', lw=1, label='guess')
plt.plot(X, sp(X), '-r', lw=2, label='spline')
plt.plot(X, sp.derivative()(X), '-g', label='slope')
plt.plot(x, y, 'ok', label='data')
plt.legend(loc='best')
plt.show()