I'm working on a small project here in my uni where I have to create a program that, based on data from a .txt file, perfmorms a certain Method (root finding, interpolaltion, derivation, etc) and pastes the results in another .txt file.
Thing is, I can't get it to work on the Cubic Spline method since the function scipy.interpolate.make_interp_spline returns an object of the BSpline class, and I need it as an array with the coefficients to be able to add it to the .txt file.
I looked around the SciPy library to look for a solution and found the interpolate.PPoly.from_spline() function but it's still not what I need; I don't know how to convert a BSpline or a PPoly class object into a Numpy array
I'm pretty much a beginner with these kinds of things so any help to make it work is appreciated.
The result of PPoly.from_spline() can be a bit confusing, but it is still probably what you need.
I'll go through a simple example to explain.
The following code builds an interpolation cubic spline of the function y(t)=t**2:
t = np.linspace(0, 1, 5)
y = t**2
bsp = make_interp_spline(t, y)
The interpolation of the five points will result in a cubic B-Spline with the knot vector [0, 0, 0, 0, 0.5, 1, 1, 1, 1] as can be seen if you print(bsp.t). This is a specific case of a clamped cubic B-Spline where the first and last four knots are multiplied.
Going into the details of the B-Spline representation is out of the scope of this answer. I recommend reading about it in the many references out there, for example in this online course by Shene.
Running the following command will give you a PPoly object as you suggested.
poly = PPoly.from_spline(bsp)
Now we can analyze the result to understand what's going on.
From the documentation of PPoly it has two relevant attributes: poly.x are the breakpoints between polynomials, and poly.c are the coefficients of the polynomials in local power basis.
The local power basis is important to notice, it means the i'th polynomial is represented in the variable (x-x[i]) (and not in the variable x as maybe one might expect). As the documentation says "the polynomial between x[i] and x[i + 1] is written:"
S = sum(c[m, i] * (xp - x[i])**(k-m) for m in range(k+1))
Now let's see how this is demonstrated in our example.
poly.c is a numpy array of (local) coefficients of shape (4, 8). The 4 comes from number of coefficients of each (cubic) polynomial. The 8 is a bit more confusing. However, if you print poly.x you will see it is [0,0,0,0, 0.5, 1,1,1,1] (taken from bsp's knot vector), so 8 is the number of intervals [x[i], x[i+1]] of domains on which the polynomials are defined.
In this example (and probably every example taken from a clamped B-Spline) the first and last three intervals are degenerate.
If we print the coefficients using print(poly.c.T) (I use the transpose of the coefficient array since it's easier to understand) we'll get the following (on my machine, your's may slightly vary):
[[ 1.62832710e-15 1.00000000e+00 2.10358047e-16 0.00000000e+00]
[ 1.62832710e-15 1.00000000e+00 2.10358047e-16 0.00000000e+00]
[ 1.62832710e-15 1.00000000e+00 2.10358047e-16 0.00000000e+00]
[ 1.62832710e-15 1.00000000e+00 2.10358047e-16 0.00000000e+00]
[-3.55271368e-15 1.00000000e+00 1.00000000e+00 2.50000000e-01]
[-3.55271368e-15 1.00000000e+00 2.00000000e+00 1.00000000e+00]
[-3.55271368e-15 1.00000000e+00 2.00000000e+00 1.00000000e+00]
[-3.55271368e-15 1.00000000e+00 2.00000000e+00 1.00000000e+00]]
As I said above, the first and last three (degenerate) intervals can be ignored in our context, so let's look at the fourth and fifth lines.
Line 4 gives the coefficients of the interval [0, 0.5] and they are basically [0, 1, 0, 0]. These correspond to the polynomial:
0*(x-0)**3 + 1*(x-0)**2 + 0*(x-0) + 0*1 = x**2
So we get that on the interval [0, 0.5] the interpolating polynomial is x**2 as expected (since a cubic spline interpolates x**2 exactly).
Similarly, line 5 gives the coefficients of the interval [0.5, 1] and they are basically [0, 1, 1, 0.25]. These correspond to the polynomial:
0*(x-0.5)**3 + 1*(x-0.5)**2 + 1*(x-0.5) + 0.25*1
And when expanding the terms we get (again as expected) x**2.
So, basically, after getting a PPoly object from a spline using from_spline() you can get a numpy array of the local basis polynomial coefficients over the i'th interval using poly.c[:, i]. If you need the non-local coefficients, you need to expand the local basis from (x-x[i]) to x as demonstrated above.