coefficients d'interpolation spline en scipy
je veux calculer les coefficients d'interpolation spline par scipy. En MATLAB:
x=[0:3];
y=[0,1,4,0];
spl=spline(x,y);
disp(spl.coefs);
et il sera de retour:
ans =
-1.5000 5.5000 -3.0000 0
-1.5000 1.0000 3.5000 1.0000
-1.5000 -3.5000 1.0000 4.0000
4 réponses
Je ne suis pas sûr qu'il y ait un moyen d'obtenir exactement ces coefficients de scipy. scipy.interpolate.splrep donne les coefficients pour les noeuds pour un B-spline. Ce que la spline de Matlab vous donne semble être les coefficients polynomiaux partiels décrivant les équations cubiques reliant les points que vous passez, ce qui me porte à croire que la spline de Matlab est une spline basée sur un point de contrôle tel Qu'une Hermite ou un Catmull-Rom au lieu d'une B-spline.
Cependant, scipy.interpolate.interpolate.spltopp fournit un moyen pour obtenir les coefficients polynomiaux partiels d'une B-spline. Malheureusement, il ne semble pas fonctionner très bien.
>>> import scipy.interpolate
>>> x = [0, 1, 2, 3]
>>> y = [0, 1, 4, 0]
>>> tck = scipy.interpolate.splrep(x, y)
>>> tck
Out:
(array([ 0., 0., 0., 0., 3., 3., 3., 3.]),
array([ 3.19142761e-16, -3.00000000e+00, 1.05000000e+01,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00]),
3)
>>> pp = scipy.interpolate.interpolate.spltopp(tck[0][1:-1], tck[1], tck[2])
>>> pp.coeffs.T
Out:
array([[ -4.54540394e-322, 0.00000000e+000, 0.00000000e+000,
0.00000000e+000],
[ -4.54540394e-322, 0.00000000e+000, 0.00000000e+000,
0.00000000e+000],
[ -4.54540394e-322, 0.00000000e+000, 0.00000000e+000,
0.00000000e+000],
[ 0.00000000e+000, 0.00000000e+000, 0.00000000e+000,
0.00000000e+000],
[ 0.00000000e+000, 0.00000000e+000, 0.00000000e+000,
0.00000000e+000]])
Notez qu'il existe un ensemble de coefficients par noeud, pas un seul pour chacun des points originaux passé. De plus, la multiplication des coefficients par la matrice de base de B-spline ne semble pas très utile.
>>> bsbm = array([[-1, 3, -3, 1], [ 3, -6, 3, 0], [-3, 0, 3, 0],
[ 1, 4, 1, 0]]) * 1.0/6
Out:
array([[-0.16666667, 0.5 , -0.5 , 0.16666667],
[ 0.5 , -1. , 0.5 , 0. ],
[-0.5 , 0. , 0.5 , 0. ],
[ 0.16666667, 0.66666667, 0.16666667, 0. ]])
>>> dot(pp.coeffs.T, bsbm)
Out:
array([[ 7.41098469e-323, -2.27270197e-322, 2.27270197e-322,
-7.41098469e-323],
[ 7.41098469e-323, -2.27270197e-322, 2.27270197e-322,
-7.41098469e-323],
[ 7.41098469e-323, -2.27270197e-322, 2.27270197e-322,
-7.41098469e-323],
[ 0.00000000e+000, 0.00000000e+000, 0.00000000e+000,
0.00000000e+000],
[ 0.00000000e+000, 0.00000000e+000, 0.00000000e+000,
0.00000000e+000]])
Le FORTRAN Polynomiale par morceaux Paquet, PPPack, a une commande bsplpp qui passe de B-spline à piecewise forme polynomiale, qui peut répondre à vos besoins. Malheureusement, il n'y a pas D'emballage Python pour PPPack en ce moment.
Voici comment je pourrais obtenir des résultats similaires à MATLAB:
>>> from scipy.interpolate import PPoly, splrep
>>> x = [0, 1, 2, 3]
>>> y = [0, 1, 4, 0]
>>> tck = splrep(x, y)
>>> tck
Out: (array([ 0., 0., 0., 0., 3., 3., 3., 3.]),
array([ 3.19142761e-16, -3.00000000e+00, 1.05000000e+01,
0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
0.00000000e+00, 0.00000000e+00]),
3)
>>> pp = PPoly.from_spline(tck)
>>> pp.c.T
Out: array([[ -1.50000000e+00, 5.50000000e+00, -3.00000000e+00,
3.19142761e-16],
[ -1.50000000e+00, 5.50000000e+00, -3.00000000e+00,
3.19142761e-16],
[ -1.50000000e+00, 5.50000000e+00, -3.00000000e+00,
3.19142761e-16],
[ -1.50000000e+00, 5.50000000e+00, -3.00000000e+00,
3.19142761e-16],
[ -1.50000000e+00, -8.00000000e+00, -1.05000000e+01,
0.00000000e+00],
[ -1.50000000e+00, -8.00000000e+00, -1.05000000e+01,
0.00000000e+00],
[ -1.50000000e+00, -8.00000000e+00, -1.05000000e+01,
0.00000000e+00]])
si vous avez la version de scipy >= 0.18.0 installée, vous pouvez utiliser la fonction CubicSpline de scipy.interpoler pour l'interpolation cubique spline.
vous pouvez vérifier la version de scipy en exécutant les commandes suivantes en python:
#!/usr/bin/env python3
import scipy
scipy.version.version
si votre version de scipy est >= 0.18.0 vous pouvez exécuter l'exemple de code suivant pour l'interpolation cubique de spline:
#!/usr/bin/env python3
import numpy as np
from scipy.interpolate import CubicSpline
# calculate 5 natural cubic spline polynomials for 6 points
# (x,y) = (0,12) (1,14) (2,22) (3,39) (4,58) (5,77)
x = np.array([0, 1, 2, 3, 4, 5])
y = np.array([12,14,22,39,58,77])
# calculate natural cubic spline polynomials
cs = CubicSpline(x,y,bc_type='natural')
# show values of interpolation function at x=1.25
print('S(1.25) = ', cs(1.25))
## Aditional - find polynomial coefficients for different x regions
# if you want to print polynomial coefficients in form
# S0(0<=x<=1) = a0 + b0(x-x0) + c0(x-x0)^2 + d0(x-x0)^3
# S1(1< x<=2) = a1 + b1(x-x1) + c1(x-x1)^2 + d1(x-x1)^3
# ...
# S4(4< x<=5) = a4 + b4(x-x4) + c5(x-x4)^2 + d5(x-x4)^3
# x0 = 0; x1 = 1; x4 = 4; (start of x region interval)
# show values of a0, b0, c0, d0, a1, b1, c1, d1 ...
cs.c
# Polynomial coefficients for 0 <= x <= 1
a0 = cs.c.item(3,0)
b0 = cs.c.item(2,0)
c0 = cs.c.item(1,0)
d0 = cs.c.item(0,0)
# Polynomial coefficients for 1 < x <= 2
a1 = cs.c.item(3,1)
b1 = cs.c.item(2,1)
c1 = cs.c.item(1,1)
d1 = cs.c.item(0,1)
# ...
# Polynomial coefficients for 4 < x <= 5
a4 = cs.c.item(3,4)
b4 = cs.c.item(2,4)
c4 = cs.c.item(1,4)
d4 = cs.c.item(0,4)
# Print polynomial equations for different x regions
print('S0(0<=x<=1) = ', a0, ' + ', b0, '(x-0) + ', c0, '(x-0)^2 + ', d0, '(x-0)^3')
print('S1(1< x<=2) = ', a1, ' + ', b1, '(x-1) + ', c1, '(x-1)^2 + ', d1, '(x-1)^3')
print('...')
print('S5(4< x<=5) = ', a4, ' + ', b4, '(x-4) + ', c4, '(x-4)^2 + ', d4, '(x-4)^3')
# So we can calculate S(1.25) by using equation S1(1< x<=2)
print('S(1.25) = ', a1 + b1*0.25 + c1*(0.25**2) + d1*(0.25**3))
# Cubic spline interpolation calculus example
# https://www.youtube.com/watch?v=gT7F3TWihvk