Create spline from knots and coefficients using scipy

Question

I'm trying to reproduce a function from a paper, which is specified only in terms of spline knots and coefficients. After finding this on stackoverflow, given a scipy interpolation object, from its knots and coefficients, I can recreate the scipy interpolation. However, the approach fails for the function specified in the paper. To reproduce a scipy interpolation I can do this:

using PyCall, PyPlot, Random
Random.seed!(5)
sp = pyimport("scipy.interpolate")
x = LinRange(0,1,50)
y = (0.9 .+ 0.1rand(length(x))).*sin.(2*pi*(x.-0.5))
t = collect(x[2:2:end-1]) # knots
s1 = sp.LSQUnivariateSpline(x, y, t)
x2 = LinRange(0, 1, 200) # new x-grid
y2 = s1(x2) # evaluate spline on that new grid
figure()
plot(x, y, label="original")
plot(x2, y2, label="interp", color="k")
knots = s1.get_knots()
c = s1.get_coeffs()
newknots(knots, k) = vcat(fill(knots[1],k),knots,fill(knots[end],k)) # func for boundary knots of order k
forscipyknots = newknots(knots, 3)
s2 = sp.BSpline(forscipyknots, c, 3)
y3 = s2(x2)
plot(x2,y3,"--r", label="reconstructed \nfrom knots and coeff")
legend()

Which provides the following as expected:

On trying to reproduce a function (image below) with specified knots = [.4,.4,.4,.4,.7] and coefficients c = [2,-5,5,2,-3,-1,2] which is supposed to produce:

With the below code and above knots and coefficients:

knots = [.4,.4,.4,.4,.7]
c = [2,-5,5,2,-3,-1,2]
forscipyknots = newknots(knots, 3)
s2 = sp.BSpline(forscipyknots, c, 3)
figure()
plot(x2, s2(x2))

I get the following (below) instead. I'm sure I'm messing up the boundary knots - how can I fix this?

Answer 1

Short answer:

The inner-knot sequence t=[0.4,0.4,0.4,0.4,0.7] and the parameters c=[2,-5,5,2,-3,-1,2] do not allow a spline to be constructed, the example contains an error (more on this later). The best you can get out of it is to remove one of the 0.4 knots and construct a quadratic (second-degree) spline as follows

tt = [0.0,0.0,0.0,0.4,0.4,0.4,0.7,1.0,1.0,1.0]
c  = [2,-5,5,2,-3,-1,2]
s2 = BSpline(tt,c,2)

This produces the following graph

Long answer:

The knot sequence in Example 3 contains only the inner knots, therefore, you need to add the boundary knots. Since you want to evaluate the spline on the interval [0,1] the full knot sequence needs to cover the points 0 and 1. The simplest is to add 0 to the beginning and 1 to the end of the sequence and replicate them as necessary according to the desired degree of the spline. A cubic (third-degree) spline would require four boundary knots (i.e. four zeros and four ones) and a quadratic spline would require three boundary knots (three zeros and three ones).

There is a problem, however. A cubic spline would require 9 parameters, while Example 3 only gives you 7. Hence, you cannot construct a cubic spline from this. With the seven parameters given, you could construct a quadratic spline, but, there, the problem is that for quadratic splines each point can only appear at most three times in the inner-knot sequence. And 0.4 appears four times (which would suggest a cubic spline). Hence, all you can do is to remove one of the 0.4 knots and construct a second-degree spline as in the short answer above.

Now I will explain what you did wrong. In the first example you obtained the knot sequence from an existing spline using knots = s1.get_knots(), which gave you knots=[0,0.02,0.04,...,0.98,1]. This sequence contains the boundary knots 0 and 1 (although only once). Hence, to construct a cubic spline, you replicated each of these three times to obtain forscipyknots = [0,0,0,0,0.02,0.04,...,0.98,1,1,1,1]. So far so good.

In Example 3, however, the knot sequence does not contain the boundary points. As you did the same as before, you ended up replicating the 0.4 and 0.7 knots three times resulting in forscipyknots = [0.4,0.4,0.4,0.4,0.4,0.4,0.4,0.7,0.7,0.7,0.7]. You cannot construct a spline on this sequence, whatever comes out of this is not a spline. What you needed instead was forscipyknots = [0.0,0.0,0.0,0.0,0.4,0.4,0.4,0.4,0.7,1.0,1.0,1.0,1.0] (which would not have worked because you do not have enough coefficients; but you could try this with your own, for instance, c = [1,2,-5,5,2,-3,-1,2,1]). To do this, you needed to add 0 to the beginning and 1 to the end of the array and only then use your newknots function.

Just as an example, a cubic spline could look like this

tt = [0.0,0.0,0.0,0.0,0.4,0.4,0.4,0.4,0.7,1.0,1.0,1.0,1.0]
c  = [1,2,-5,5,2,-3,-1,2,1]
s2 = BSpline(tt,c,3)