[SciPy-Dev] Splines

[Pauli Virtanen]

Hi,

Charles R Harris <charlesr.harris <at> gmail.com> writes:
> There have been several threads on the list about splines
> and consolidation of splines. For instance, there are several
> uniform spline implementations for images and signal processing,
> various low level functions in fitlib that are unexposed, and
> perhaps useful altenatives to b-splines for some applications
> like straight interpolation. For myself, I've started implementing
> several functions in pure Python with an eye to converting them
> to Cython once the interface and documentation is in place,
> mostly for doing things that fitpack doesn't do because it is
> very integrated, as opposed to supplying a basic toolset.
> As part of this project, I'd like to get some feedback
> on which functions people use most and what features they
> would like to see. I'm not interested in the high level classes
> at this point, either the current classes or combo functions
> like fpcurf, but rather a collection of good lower level function
> that could be used to implement the higher level functions
> in a more flexible way. Thoughts?

Great! I was going to start on this for 0.13.0, but this should
speed things up considerably :)

Overall, I think what we need is are (i) a well-specified spline
format, and (ii) solid basic functions for evaluating and
manipulating them, (iii) ditto for tensor product splines.

How to construct the splines (interpolation, smoothing, etc.) should
be considered as a separate problem. We can turn to FITPACK for
smoothing splines, but it should not be used for interpolating 
splines.

    ***

Some misc thoughts on this:

* The spline data format should be documented and set in
stone as a first step. Users (and future developers) will
want to toy around with it.

Also, the data format for tensor product N-dim splines needs
to be set. They are what we are missing, and what people are
constantly asking for. We don't want them turn to
`ndimage.map_coordinates`, which is clunky to use.

The Fitpack tck format looks like this:

https://github.com/pv/scipy-work/commits/spline-unify

Currently, there's also a second B-spline data format used in
scipy.interpolate with different padding, but we may want to
stick with the FITPACK one, it's probably as good as any.

* Functions for splines with varying x-coordinates are needed.
Uniform grid splines would be a nice bonus as a speed-gain
optimization.

* The 1-D spline routines should be able to work over an
arbitrary axis of multidimensional data. Even better if this
can be done without reshaping and copying the input data
(e.g. with Numpy iterators).

This sounds like providing strided 1-D loops for heavy lifting,
and bolting array iterators on top.

* Functions for integration + differentiation of splines
as as abstract objects would be useful. For efficiency,
evaluation of derivatives & integrals probably might need
to be provided separately.

* For tensor product splines, evaluation on a scattered point
set + on a grid would be useful.

* It will probably be easiest to start from a clean slate,
rather than trying to reuse scipy.interpolate.

* FITPACK should not be used as a basis for this work, no
need to use ancient FORTRAN 77 code for simple stuff. We can
use its routines for generating smoothing splines, though.

* Routines for constructing interpolating splines --- I think most
of the time people will use these for simple gridded data interpolation
rather than smoothing. FITPACK's knot selection is nice when it works,
but often it doesn't (or requires careful fiddling), so we should have
something simple and robust as a default algorithm.

* I'm not sure what to do with the various boundary conditions
and out-of-bounds value handling. It's probably best to leave room
for various ways to do this...

* Scipy also has two implementations of piecewise polynomials
--- these should be consolidated into one, too.

Cheers,
Pauli