How to set manual knots in scipy.interpolate.splprep in case of periodic boundary conditions?

Question

or

Why splprep does not work with its own knots?

I am trying to figure out how to set the knots in scipy.interpolate.splprep. In a non-periodic case I succeed, i.e. I can reproduce this SE example.

In case of periodic boundary conditions (PBC) I have a problem, though. Here scipy.interpolate.splprep doesn't even work with its own knots as this example shows:

import numpy as np
from scipy.interpolate import splev, splrep
import scipy

print "my scipy version: ", scipy.__version__

srate = 192000.
freq = 1200.

timeList = [ i / srate for i in range( int( srate / freq + 1 ) ) ]
signal = np.fromiter( ( np.sin( 2 * np.pi * freq * t ) for t in timeList ), np.float )

spl = splrep( timeList, signal, per=1 )
knots = spl[0]
spl = splrep( timeList, signal, t=knots, per=1 )

gives:

my scipy version:  1.0.0
Traceback (most recent call last):
  File "splrevtest.py", line 15, in <module>
      spl = splrep( timeList, signal, t=knots, per=1 )
  File "/somepath/python2.7/site-packages/scipy-1.0.0-py2.7-linux-/python2.7/site-packages/scipy-1.0.0-py2.7-linux-x86_64.egg/scipy/interpolate/fitpack.py", line 289, in splrep
      res = _impl.splrep(x, y, w, xb, xe, k, task, s, t, full_output, per, quiet)
  File "/somepath/python2.7/site-packages/scipy-1.0.0-py2.7-linux-x86_64.egg/scipy/interpolate/_fitpack_impl.py", line 514, in splrep
      raise _iermess[ier][1](_iermess[ier][0])
  ValueError: Error on input data

Moreover, if you plot the first and the last knots like:

print timeList[-1]
print knots[:4]
print knots[-4:]

you get

>> 0.000833333333333
>> [-1.56250000e-05 -1.04166667e-05 -5.20833333e-06  0.00000000e+00]
>> [0.00083333 0.00083854 0.00084375 0.00084896]

meaning that there are six dots out of the actual data range, three before and three after. This is probably OK as we have PBC and the knots are identical to knots inside the data range modulus the period, but strange anyway. (BTW, the example does not even work when setting per=0.) Moreover, setting points outside the data range fails if I set the dots manually. It works, even with per=1, if the dots are well inside the data range. Like e.g.:

import matplotlib.pyplot as plt
import numpy as np
from scipy.interpolate import splev, splrep
import scipy

x = np.linspace( 0, 1, 120 )
timeList = np.linspace( 0, 1, 15 )
signal = np.fromiter( ( np.sin( 2 * np.pi * t +.3 ) for t in timeList ), np.float )
signal[ -1 ] -= .15
myKnots=np.linspace( .05, .95, 8 )
spl = splrep( timeList, signal, t=myKnots, per=1 )
fit = splev(x,spl)

fig = plt.figure()
ax = fig.add_subplot( 1, 1, 1 )
ax.plot( timeList, signal, marker='o' )
ax.plot( x, fit , 'r' )
for i in myKnots:
    ax.axvline( i )
plt.show()

providing:

wehre one can also see that the last point is ignored in PBC.

So what is scipy.interpolate.splprep actually doing here, and why is it not accepting a similar knot construction for manual set knots if per=1. Is it a bug?

I hope the answer to the last question is 'no', as otherwise I'd be wrong asking this here.

Answer 1

In order to make splrep work with your own knots, you should change:

spl = splrep(timeList, signal, t=knots, per=1)

to use the interior knots:

spl = splrep(timeList, signal, t=knots[4:-4], per=1)

This interface isn't so intuitive but it appears in the documentation (albeit under the task parameter and not under t) (my emphasis):

If task=-1 find the weighted least square spline for a given set of knots, t. These should be interior knots as knots on the ends will be added automatically.

The additional knots are also documented there. This is consistent with the general connection between the number of knots and the number of coefficients: number_of_knots = number_of_coefficients + degree + 1.

The last point in the PBC construction is ignored (as shown in your figure) because the periodic definition requires that the value of the last point in the period be equal to the first. This is documented under the per parameter where it says "Values of y[m-1] and w[m-1] are not used." (again, not a very intuitive interface I guess, but that's probably how the Fortran code was implemented).