Search code examples
pythonscipybspline

Cannot reconstruct scipy's spline

I am constructing a spine that interpolates well some data I have and I want to store this in a database. I thought that this would be straightforward since I can just take the knots and the coefficients and just create a new spline with them. However, this goes horribly wrong. Here is a minimal (non) working example.

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

xs = np.linspace(0,2,101)
ys = xs*np.cos(2*np.pi*xs)

knotNumber = 9
knots = np.linspace(0,2,knotNumber)[1:-1]
spline = interpolate.LSQUnivariateSpline(xs,ys,knots,k=3)
re_spline = interpolate.BSpline(spline.get_knots(),spline.get_coeffs(),3)

sp_ys = spline(xs)
re_ys = re_spline(xs)

plt.plot(xs, ys, '-', color='black');
plt.plot(xs, sp_ys, '-', color='blue');
plt.plot(xs, re_ys, '-', color='red');

This results to the following figure:

splines

As you can see the LSQUnivariateSpline blue curve covers completely the black. The red curve that comes from the spline I created seems completely unrelated.

Is there something obvious I am missing here? How can I reconstruct the spline?

Solution

  • BSpline needs k extra knots at the start and end (the fixed points); different than (LSQ)`UnivariateSpline.

    I found this by looking at the result of

    interpolate.splrep(xs, ys, k=3, task=-1, t=knots)

    which shows 3 extra zeros at the start of the outputted knots, and 3 extra twos at the end:

    [0. 0. 0. 0. 0.25 0.5 0.75 1. 1.25 1.5 1.75 2. 2. 2. 2. ].

    (The amount of extra endpoints obviously depends on the degree.)

    With that given, the following should work:

    xs = np.linspace(0,2,101)
ys = xs*np.cos(2*np.pi*xs)

k = 3
knotNumber = 9
knots = np.linspace(0,2,knotNumber)
spline = interpolate.LSQUnivariateSpline(xs, ys, knots[1:-1], k=k)
bknots = np.array(k * [knots[0]] + spline.get_knots().tolist() + k * [knots[-1]])
re_spline = interpolate.BSpline(bknots, spline.get_coeffs(), k)

sp_ys = spline(xs)
re_ys = re_spline(xs)

    (Or with any other, prettier, way to prepend and append the extra knots to spline.get_knots().)