[SciPy-User] help interpreting univariate spline

[Elliot Hallmark]

On Fri, Apr 23, 2010 at 6:53 PM, Elliot Hallmark <permafacture at gmail.com> wrote:

> I am wanting to use a scipy interpolation routine to generate
> polynomials approximating function I have sampled.
> ...
> import numpy as np
>from scipy.interpolate import splrep, UnivariateSpline

>x = np.linspace(1, 10, 100)
>y = 1/x     #approximate a hyperbola
>g = UnivariateSpline(x, y, w=None, k=3, s=.0005)
>print g.get_knots()
>print g.get_coeffs()

>and the output was:
>[  1.           2.18181818   3.27272727   5.54545455  10.        ]
>
>[ 0.98462432  0.66575572  0.423       0.25461626  0.13929847  0.11763985
 0.09929961]

>That is 5 knots and 7 coefficients for a degree=3 spline.  naively, I
>was expecting 4 coefficents for each interval between knots.

I found this post:

http://mail.scipy.org/pipermail/scipy-user/2009-June/021572.html

This seems very round-a-bout and ignores the piecewise nature of the
spline.  I mean, really? is there no direct or documented way to get
the polynomial coefficents from a spline in python?  If not, then this
code is the best I have found.  For my purposes I have modified it
(see below)


"""
try_spline_coeff.py

explorations into getting polynomial coefficents from a spline.

Modified from http://mail.scipy.org/pipermail/scipy-user/2009-June/021572.html
by josef.pktd at gmail

please email offonoffoffonoff at gmail with improvements to this technique.
"""

import numpy as np
import scipy as sp
from scipy import interpolate

def cubice(x,a):
    "evaluate a cubic"
    return (a[:,0]+a[:,1]*x+a[:,2]*x**2+a[:,3]*x**3)

def deriv2coeff(x2, tck):
    '''construct a piecewise cubic polynomial from derivatives of
spline at given knots (x2)

	output is an array of knot abscissa and coefficent of poly for that section.

	out = (x_knot,a,b,c,d)'''

    if tck[2] != 3:
	print "interpolation must be cubic!"
	return np.inf
    coeffr = np.vstack(interpolate.spalde(x2,tck))
    coeffr[:,3] /= 6.
    coeffr[:,2] = (coeffr[:,2] - 6*coeffr[:,3]*x2)/2.
    coeffr[:,1] = coeffr[:,1] - 2*coeffr[:,2]*x2 - 3*coeffr[:,3]*x2**2
    coeffr[:,0] = coeffr[:,0] - coeffr[:,1]*x2 - coeffr[:,2]*x2**2 -
coeffr[:,3]*x2**3
    return np.concatenate((np.resize(x2,(x2.size,1)),coeffr), axis=1)

def calc_poly(xs, coeff):
	'''determine which piece of the piecewise polynomial the given x
value corisponds to
	and return the value of the approximation based on evaluating the
correct polynomial'''
	result = np.zeros((xs.size,4))	
	for i,x in enumerate(xs):
		a = b = c = d = 0
		for x_knot,a1,b1,c1,d1 in coeff:
			if x>x_knot:		#x_knot is in ascending order (required by splrep also)
				a=a1
				b=b1
				c=c1
				d=d1
		result[i] = a,b,c,d
		#print "single: ",x, result[i]
	return cubice(xs,result)

npoints = 51
x = np.linspace(1, 20, 2*npoints)
y = 1/x
tck = interpolate.splrep(x, y, k=3)

#calculate coefficents at original x points,
# would using tck knot values be more faithful to original spline?
coeffr = deriv2coeff(x, tck)

#print coeffr

#now, compare the polynomial and original spline over a different dataset
x2 = np.linspace(3, 18, 40)
y_spline = interpolate.splev(x2, tck)
y_polynomial = calc_poly(x2,coeffr)
print "max error: ", np.max(np.abs(y_spline-y_polynomial))

#using the actual knots of the spline
coeffk = deriv2coeff(tck[0][3:-3], tck)
y_polynomial = calc_poly(x2,coeffk)
print "max error using spline's knots: ", np.max(np.abs(y_spline-y_polynomial))