Coverage remained the same when pulling 8bfe266 on ev-br:bsplines into fde1826 on scipy:master.

The last commit fixes the out-of-bounds behavior to be almost consistent with the b-spline definition in the literature: the interpolation interval is closed, so that the spline is left-continuous at the right-hand endpoint; outside of the interval it is zero, unless extrapolation is requested.

A corner case is nan handling: at the moment, nan is not treated specially, so that b(nan)=0.

Coverage remained the same when pulling 28ec512 on ev-br:bsplines into fde1826 on scipy:master.

The from_fitpack_tck and get_fitpack_tck are potentially confusing. Does it make sense to have two incompatible tck representations? Can we have another naming convention for b-splines?

Have you considered how to replace the spline implementations in the signal and ndimage modules with this?

Have you seen this thread: http://thread.gmane.org/gmane.comp.python.scientific.devel/17349 ? Both @pv and @charris were also working on / thinking about b-splines.

@rgommers in fact this code heavily uses things done by Pauli, borderline shameless stealing.

Thanks for the link, I'll have w closer look --- I only saw gh-1408 and scattered older ML discussions. The point of this PR is mostly to provide a stab at an implementation and to have something concrete to discuss. I'll happily scrap all this if there's a better implementation and/or if there are objections to details of the format (tck padding etc).

I think this can be used for both replacing much of the fitpack in the interpolate module (interp1d, shudder), either as a replacement or making user-facing interp1d, UnivariateSpline use these under the hood. Same for cardinal splines in signal and ndimage.
Doing this all in a single PR is a bit too much, so for the moment I'm trying to concentrate on the low level functionality, then move on the interpolating splines and fitting (have bits and pieces, not ready for a PR), and then start replacing and consolidating the user-facing code. It all relies on a consensus on the basic format --- hence this PR.

I'll address your comments, but for the moment I've run out of time for today, time for a New Year :-).

Agreed that signal/ndimage is not for this PR. My question on tck is not an objection, more like something I was wondering about.

Happy New Year!

tck: it's ugly, I agree. The problem, as far as I understand it, is with how exactly to define a spline, esp at the edges of the knot vector. Dierckx's convention is probably documented is his book, which costs more than a hundred quid and is not widely available from what I can see. Extracting it from the fitpack code is fun :-). My convention here is, I think, compatible with the reference in the docstring. Will attempt at documenting it better. Whether this is acceptable/typical, I'd appreciate an input on.

Happy New Year!

Coverage remained the same when pulling e3a96ce on ev-br:bsplines into fde1826 on scipy:master.

Have added a couple of commits with improvements to documentation [shameless stealing from Pauli, again], small tweaks in the basis_element signature, and tests.

Also timed the evaluations (in a rather ad hoc way) against fitpack:
http://nbviewer.ipython.org/gist/ev-br/8260758
Have to admit, I don't really understand these timings. For scalar input, splev wins hands down (python overhead of BSpline.__call__?). For 1D input it seems to depend on whether the input is sorted: if it is, fitpack shines by a factor of 3--5; otherwise, it loses by a factor of about two. At the same time, cython evaluations are barely sensitive to the ordering of input.
Clearly, fitpack does something smart somewhere (and I do something dumb). Will try investigating.

Coverage remained the same when pulling 05bf090 on ev-br:bsplines into fde1826 on scipy:master.

Sorry for the delay. I have one bigger comment about the spline representation:

I really think we should go with the representation used in FITPACK. The argument for this is:

• It's what was available previously, and we retain tck compatibility with splrep et al.
• It's the simplest and most general representation (if you assume zero coefficients outside the array). The finite interval spline can be obtained by adding coinciding knots at the endpoints (which is what __init__ currently does).
• The only small quirk of the FITPACK repr is that the c array is padded at the end to the same size as t rather than to the number len(t)-k-1 of the actual basis functions. (We could avoid this by allowing any size c array an assuming zero coefficients out of bounds.)
• The assertion p = BSpline(t, c, k); assert_allclose(p.t, t); assert_allclose(p.c, c); assert_allclose(p.k, k) is retained, as for PPoly and BPoly. I think having this assertion be true is somewhat important in the general case.
• We can drop the get_fitpack_* functions completely.

A separate constructor could be added for the finite interval case.

I'll just add a few comments about things that should be kept in mind for the future.

• Tensor products. This will in general require several copies of tck.
• Vector valued splines, needed for curves.
• Periodic splines, needed for closed curves.

These can probably be built on top of 1D, scalar valued splines, but some thought up front may make that easier to do.

@charris: Tensor product support should be done similarly as for piecewise polynomials in gh-3104 (not merged yet, if you have ideas how to improve that, it would be appreciated).

This PR already supports vector-valued splines, if I understand what you mean correctly.

Periodic splines: for evaluation, does this mean that the knot set is treated as periodic, identifying t[0] and t[-1]?

Vector valued splines is essentially a stack of splines with the same knot points and end conditions, i.e., a map of a line segment into a multidimensional space. If they are already there, great. It also makes solving for multiple rhs more efficient.

Yes, the coefficients are periodic. This mostly affects the spline construction, it is a boundary condition. I've looked at the Sherman-Morrison formula and its generalizations for doing this in the uniform spline case, but I don't know what the current state of the art is.

I'll add that when I was looking at splines, my impression was that cubic splines were nearly optimal for most problems, apart from the results of differentiation and integration, and some boundary conditions -- not-a-knot, etc. -- made more sense with that restriction. So it might be worth while offering some options that are only available for cubic splines.

OK, time to stop this from bitrotting completely. Rebased, squashed, addressed the review comments.
This version has tck format directly compatible with fitpack.

Have also updated the notebook with timings (as crude as they are). Which indicate this is still very much WIP: being two to three times slower than fitpack is no good at all. The prime suspect is my half-baked evaluation routine, will investigate.

Coverage remained the same when pulling 980ef5f on ev-br:bsplines into 87df1db on scipy:master.

The last commit adds the constructor for interpolating splines. Several boundary conditions are supported (natural, not-a-knot, explicit derivatives, periodic). All linear algebra uses banded matrices. Still WIP: periodic boundary conditions need a bit more work.

Coverage increased (+0.06%) when pulling 934a316 on ev-br:bsplines into 93656a8 on scipy:master.

Coverage increased (+0.06%) when pulling 934a316 on ev-br:bsplines into 93656a8 on scipy:master.

Thanks. Brief look at the last commits seems to indicate that this is close to the minimal set of feature completeness. I'll hope to look at this closer soon(TM) (probably not within the next week however...)

Coverage increased (+0.09%) when pulling 97e111b on ev-br:bsplines into 93656a8 on scipy:master.

OK, pushed an update where I believe I addressed @pv's review comments. Apparently Github does not hide outdated comments anymore. Relevant commits start from e6c5b72.

The main change is that splrep and splprep again return tck tuples and not BSpline objects (0131d72). I've to say I actually like it more. The interface with spl*-functions is still not perfect for multidimentional splines, but at least the logic seems clear: splint(BSpline, ...) is consistent with BSpline.integrate(), and splint(tck, ...) is backwards compatible. So far I did not implement a keyword switch, as suggested in #3174 (comment), but it's easy to bolt on.

Also, tuple-like unpacking of BSplines is gone. This was needed for backwards compat if splrep returns a BSpline instance, to support legacy code which does t, c, k = splrep(x, y). But as long as we're not trying to sneak BSplines in place of tuples, I'd rather avoid that sort of magic. Knots, coefficients and the spline order are anyway available as b.tck or as individual attributes b.t, b.c b.k.

Ok, looks good to me currently, and OK to merge from my side.
If we come up with better names than make_* or other tweaks before 0.19.0, we can make those.

OK, I'm going to hit the green button in a couple of days, unless further comments.

Great - thanks very much for your persistence on this one.

Thanks @ev-br, let's get it in!

Thanks all for the code, for the reviews, for encouragement and for evolving this from where it started to a useful state!

Yay! Great job & persistence @ev-br.

Maybe a PR to say "this is the way of the future" and clarify some of the existing text on splines in the roadmap would be a good idea?

Here's a PR: #6715

Note that this introduced some performance regressions:
https://pv.github.io/scipy-bench/regressions.xml
e.g. https://pv.github.io/scipy-bench/index.html#interpolate.Interpolate1d.time_interpolate?commits=a247ff7f61e6d88a02db3f451c52b378a3238c29&idx=3

Or more precisely, the slinear and zero in interp1d got slower --- the cubic and quadratic got faster (the previous scaling with full matrices was of course ridiculous :).
Maybe there's some trivial special casing possible.

Hmm, evaluations for slinear and zero got about an order of magnitude faster, (time_interpolate_eval benchmark), but construction (time_interpolate benchmark) got an order of magnitude slower, or am I misreading it?