Cubic spline interpolation

[Tomas Lycken]

Hi,

Is there a (maintained) package somewhere with cubic spline capabilities? I need something that fulfills the following requirements:

* Scalar-valued functions of one variable, f(x), specified on uniform or non-uniform x-grids
* Scalar-valued functions of two variables, f(x,y), at least specified on uniform grids that don't need to have the same spacing in x and y (i.e. rectangular, but not necessarily quadratic, grid cells)
* Evaluation of function value and evaluate up to at least second order derivatives, i.e. both f'x, f'y, f'xx, f'xy and f'yy in the 2D case

The only packages I've found that seem to approach this functionality are

* https://github.com/timholy/Grid.jl - only up to quadratic splines, as far as I can tell from the readme; also unsure on if it can evaluate second order derivatives

* https://github.com/gusl/BSplines.jl - only 1D interpolation, and base splines rather than exact cubic splines

* https://github.com/EconForge/splines - Specifies Julia 0.2- in its require and hasn't been touched in four months => I doubt it works with my Julia installation which is at master. It would also probably take quite a lot of work to learn to use it, since it has no documentation at all.

Are there any others that I've missed? Is there any non-official effort toward creating this functionality in Julia?

Thanks in advance,

// Tomas

[Tim Holy]

Adding support for cubic splines to Grid should be fairly straightforward.
Adding support for 2nd derivatives might be a bit harder, since in dimensions
higher than 1 you really need the Hessian. If you're interested in
implementing it, you'd basically need to create set_hessian_coordinates (like
set_gradient_coordinate) and add an `hcoef1d` field to InterpGridCoefs.

--Tim

[Tim Holy]

I should also add that unless you need a continuous second derivative,
quadratic (while not as popular) is IMO nicer than cubic in several ways---for
example, by requiring fewer points (27 rather than 64 in 3 dimensions), it's
faster to evaluate. I suspect the main reason quadratic isn't popular is that
it's "non-interpolating" (in the terminology of P. Thevenaz, T. Blu, and M.
Unser, 2000), but the infrastructure in Grid (based on interp_invert!)
compensates for that.

--Tim

On Wednesday, March 26, 2014 06:26:12 AM Tomas Lycken wrote:

[Tomas Lycken]

I'd happily use quadratic interpolation if it weren't for that exact limitation - I do need continuous second derivatives in the 2D interpolation...

I've never actually implemented spline interpolation myself, so I'd probably have to spend quite a lot of time delving into the current code to understand how it works and how to add to it. Might be good exercise, but unfortunately not something I have time for at the moment. It's an interesting learning project, though, so I'll definitely consider coming back to it later =)

The reason I started this thread was that I'm currently working on a quite large project (my masters thesis project) in C++, but after realizing today that my current codebase has some deeply nested memory management issues that I'm having difficulties troubleshooting, I wanted to see if it would be possible to port the project to Julia in a relatively short amount of time. For that purpose, I wanted to see if all the stuff I'm doing through external C++ libraries was implemented in Julia packages as well, so that I wouldn't have to write any Julia code for stuff that I haven't had to write C++ code for. Cubic splines was the only missing piece of the puzzle :P However, implementing a cubic spline interpolation routine is, unfortunately, well out of scope for my thesis, so it will have to wait until I have time to spend on it. Until then, I'd better get back to those segfaults...

// T

[Simon Kornblith]

You could write a C interface to the C++ spline library and then ccall that from Julia. That's probably not too much work, but I can't promise that the interface plus the port to Julia would be less work than fixing the segfaults.

Simon

[Stefan Schwarz]

Hello Tim,

this is maybe not the right place to ask this question,

but may you elaborate on that matter, that interp_invert!

compensates "the non-interpolating" issue of 

the quadratic?

Thanks in advance

Stefan

[Tim Holy]

This is spelled out in the paper
P. Thevenaz, T. Blu, and M. Unser (2000). "Interpolation Revisited."
IEEE Transactions on Medical Imaging, 19: 739-758.
Basically, the idea is that interpolation computes values from a set of
coefficients:
y(x) = sum_i coef_i * f_i(dx_i)
where dx_i comes from a displacement of x. For an "interpolating" function,
the coefficients are the actual on-grid values of the function you're trying to
interpolate. However, you can generalize the notion to "non-interpolating"
expressions, where the coefficients are distinct from the on-grid values. You
can pretty quickly convince yourself that this is necessary, for example, if
you want to have 3-point quadratic interpolation with continuous derivatives.
The coefficients turn out to be related to the on-grid values via a simple
tridiagonal matrix equation (which is fast to invert).

Explicitly, for quadratic interpolation you can convince yourself that the
only expression with continuous derivatives at the half-point between integers
is
y(alpha) = (alpha-0.5)^2*ym/2 + (.75-alpha)^2*y0 + (alpha+0.5)^2*yp/2
where |alpha| <= 0.5 and ym,y0,yp are the interpolation coefficients. Plug in
alpha=0, and you can see that you don't recover y0. Hence it's "non-
interpolating."

--Tim

[Tomas Lycken]

@Tim Holy, I've started thinking about how I should go about implementing cubic spline interpolation for Grid.jl over my Easter break - as a start, I'm trying to read through the code that does the quadratic interpolation and understand how all the parts of the library work together, and how I should hook into the current code structure. Of course, questions keep popping up ranging anywhere from "this seems to be dead code, is that intentional or am I missing something?" to "I have no idea what this function does, would you mind explaining it?" What is the best place for me to ask such quesitons? Here on the mailing list, in issues (or one single discussion issue) on Github, or somewhere else?

Thanks in advance!

// Tomas

[Tim Holy]

Probably best is to file an issue in Grid.

Thanks very much for tackling this!

--Tim

On Saturday, April 19, 2014 04:42:06 PM Tomas Lycken wrote:
> *@Tim Holy*, I've started thinking about how I should go about implementing