Cubic spline interpolation get coefficients

Question

Say I have two arrays in python and I wish to get (and actually use) the cubic spline interpolation between those points. (IE: I wish to integrate the function). I would strongly prefer a way using numpy scipy.

I know about scipy.interpolate.interp1d. However that only allows me to evalute the points, for say the very simply function of:

Now I could do something simply like:

import numpy as np
import scipy.interpolate
import matplotlib.pyplot as plt

y = np.array([0,2,3,4,8,10,12,12,12,10,9,8,7,1,0,0,1,2])
x = np.array(range(len(y)))
xvals = np.linspace(0, len(y)-1, len(y)*100, endpoint = False)
func = scipy.interpolate.interp1d(x, y, kind = "cubic")
yvals = func(xvals)
plt.plot(xvals,yvals)
plt.plot(x,y, "o")

However I wish to further process this cubic spline (ie I need to get the integration).. For manual doing things I need to get the factors, so:

a_i * x^3 + b_i * x^2 + c_i * x + d_i where i goes from 0 to n/3 

(n = number of elemetns - this is just the definition of the ith cubic)

I hence expect a list of tuples (or 2d array) describing all splines. - Or a way to get the ith cubic, and really, really would love to get have a convenience "x-to-i" to find in which spline I am currently.

(Though of course this latter problem is a simple search for first value larger than reference in a sorted list - I could do that by hand easily if need be).

Answer 1

For interpolation, you can use scipy.interpolate.UnivariateSpline(..., s=0).

It has, among other things, the integrate method.

EDIT: s=0 parameter to UnivariateSpline constructor forces the spline to pass through all the data points. The result is in the B-spline basis, you can get the knots and coefficients with get_coefs() and get_knots() methods. The format is the same as used in FITPACK, dierckx, on netlib. Note though that the tck format used internally by interp1d (which relies on splmake ATM) and UnivariateSpline (alternatively, splrep/splev) are inconsistent.

EDIT2: you can get a piecewise polynomial representation of a spline with PPoly.from_spline --- but this is not consistent with iterp1d. Use splrep(..., s=0) to get an interpolating spline, then convert the result.