Spline con vincoli sul bordo

Question

Ho misurato i dati su una griglia tridimensionale, ad es. f(x, y, t). Voglio interpolare e uniformare questi dati nella direzione di t con spline. Attualmente, faccio questo con 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) 

Tuttavia, ho bisogno la spline di avere derivati ​​di fuga a t=0. C'è un modo per far rispettare questo vincolo?

Answer 1

La cosa migliore che posso pensare è di fare una minimizzazione con un vincolo con scipy.optimize.minimize. È abbastanza facile prendere la derivata di una spline, quindi il vincolo è semplicemente. Vorrei usare una spline regolare (UnivariateSpline) per ottenere i nodi (t), e tenere i nodi fissi (e il grado k, ovviamente) e variare i coefficienti c. Forse c'è un modo per variare le posizioni dei nodi, ma lo lascerò a te.

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

E poi abbiamo generare alcuni dati falsi che dovrebbe avere pendenza iniziale zero e testarlo:

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

Answer 2

Il tuo esempio non funziona (il pitone 2.7.9), quindi ho solo disegnare la mia idea:

1. Calcolare sp
2. prendere il derivato tramite sp.derivative e valutare nei momenti rilevanti (probabilmente gli stessi tempi in cui hai misurato i tuoi dati)
3. Imposta i punti rilevanti su zero (ad es. il valore su t = 0)
4. Calcola un'altra spline dai valori derivati.
5. Integrare la funzione spline. Immagino che dovrai farlo numericamente, ma non dovrebbe essere un problema. Non dimenticare di aggiungere una costante per ottenere la tua funzione originale.

Answer 3

Ecco un modo per farlo. L'idea di base è ottenere i coefficienti di una spline con splrep e quindi modificarli prima di chiamare splev. I primi nodi nella spline corrispondono al valore più basso nell'intervallo di valori x. Se i coefficienti che corrispondono ad essi sono impostati uguali tra loro, ciò appiattisce completamente la spline a tale fine.

utilizzando gli stessi dati, i tempi, x, y, come nel tuo esempio:

# set up example data 
data = np.arange(1000).reshape((5, 5, 40)) 
times = np.arange(data.shape[-1]) 
y = 3 
x = 3 

# make 1D spline 
import scipy.interpolate 
from pylab import * # for plotting 
knots, coefficients, degree = scipy.interpolate.splrep(times, data[y, x]) 
t = linspace(0,3,100) 
plot(t, scipy.interpolate.splev(t, (knots, coefficients, degree))) 

# flatten out the beginning 
coefficients[:2] = coefficients[0] 
plot(t, scipy.interpolate.splev(t, (knots, coefficients, degree))) 
scatter(times, data[y, x]) 
xlim(0,3) 
ylim(720,723) 

Blu: i punti originali e spline attraverso di loro. Verde: spline modificata con derivata = 0 all'inizio. Entrambi sono ingranditi all'inizio.

plot(t, scipy.interpolate.splev(t, (knots, coefficients, degree), der=1), 'g') 
xlim(0,3) 

chiamata splev(..., der=1) per tracciare la derivata prima. La derivata inizia da zero e supera un po 'il limite in modo che la spline modificata possa recuperare (ciò è inevitabile).

La spline modificato non passa attraverso i primi due punti su cui è basata (esso colpisce ancora tutti gli altri punti esatti). È possibile modificare questo aggiungendo un ulteriore punto di controllo interno accanto all'origine per ottenere sia una derivata zero che passare attraverso i punti originali; sperimenta i nodi e i coefficienti finché non fa ciò che vuoi.