Splines

[Pablo Winant]

Hi,

Maybe some mildly ambitious project that we could start now, would be to create some good spline routines.

There are some existing libraries in Fortran (notably fitpack) and in C (einspline).
Python/Scipy has some wrappers around fitpack, and there has been some discussion on the scipy mailing list this spring to create a new implementation in Cython (a compiled dialect of python so to say), but it hasn't produced any output yet, that I am aware of. There is also a Julia package called BSplines, that is currently restricted to dimension 1.

Dolo already includes cubic splines with derivatives and extrapolation, on regular grids. These are generated from a template file that generate Cython code, which is very similar in essence to the einspline code, except that it can handle dimensions >3. Templating is used only for the evaluation, while traditional coding is used for the filtering step. It would be very easy to adapt the templating, to produce Julia code instead. It could be generated once for all, or replaced later by Julia macros. That would give us (natural) cubic splines in Julia with little work to do.

We could then work on improving the generality of the approach, by allowing for higher order splines, non regular grids, or arbitrary tensor product of linear bases (for instance splines X linear), but the cubic splines (regular and natural) should already cover a lot of cases very fast. We should also monitor the other ongoing efforts so that we don't duplicate efforts to much.

What do you think ?

Pablo

[Spencer Lyon]

Would the primary use of the splines be function approximation and interpolation?

[Pablo Winant]

That would be my use for it. I guess regression splines would also be useful.

[Derek Thomas]

Hi,

Spencer asked me to comment on this discussion.  I have been actively involved in the development of a library for what we call hierarchical spline forests (hsf).  The primary motivation for this is isogeometric analysis (IGA) where we use splines as the basis functions for analysis.  While some of the features in the library aren't necessary for your projects such as adaptive local refinement and the arbitrary topology forest constructions, the core of the library is very flexible and general and I think it would provide some useful features for you.  It is capable of representing splines of arbitrary order and dimension.  It also provides some significant unique features.

The primary concept/feature that you should consider is Bezier extraction (see Borden 2011 DOI: 10.1002/nme.2968).  This is a method for constructing splines without recursion using the Bernstein polynomials and the knot vector for each dimension.  It generalizes nicely to tensor product constructions and can be optimized.  One significant feature of our hsf library is a cache system for the Bezier extraction operators.  This is especially important for high dimension, high order, potentially nonuniform splines.

As it currently stands, the code would need to be expanded slightly for the cases with more than 3 dimensions.  The core structures make no assumptions about dimension, only the interface needs to be developed.  It would also be necessary to build up some interpolation functions.  None of this represents a significant difficulty and we would be happy to work with you to incorporate this functionality.  An additional feature is that there is a python interface to the majority of the code that would be needed for this project (thanks to Spencer).  I hope this helps,

Derek Thomas

[Pablo Winant]

Hi,

Thank you Derek for these informations, and for participating to that discussion.

I am not familiar with isogeometric analysis, but I'm very curious about it. To date, the level of sophistication finite elements approach in economics is lagging well behind what is done in physics. It is partly due to inherent complications arising from economic modellng, for instance the high dimensionality nature of certain problems or the impossibility to define an efficient mesh without knowing the "geometry" of the solution to be determined. On the other hand many hand there are still many approaches (including mesh refinement) that have not been exhausted yet. Thinking about that, I have the (yet unproven) intuition that your library could actually be very useful. And certainly worth investigating.

Please tell us when your library is opensourced. At that point in time, we'll be happy to test it, use it, or to help extending it if needed.  In the meantime, good luck with your submission !

Best,

Pablo

[Pablo Winant]

I have made a tentative implementation for cubic splines in Julia here: https://github.com/EconForge/splines_julia . This is very preliminary, incomplete and inelegant, but since I won't be making much progress on it in the following weeks, I thought it was better to share it now.

Here is what it contains:

- a function mapping function values with the spline coefficients (splines_filter.jl)
- functions to evaluate the splines at the given points (csplines.jl).  This file is programmatically generated using the script gen_splines.py and the template csplines.jl.in  Since that part of the code was ported from python I have used a Julia macro to locally use 0-based indexing. The macro itself is in zero_indexing.jl
- an example file (test.jl)

I have kept the interface to the bare minimum, there is a lot of room for improvement. With very minimal effort it can be extended to support the computation of derivatives since the code is already there on the python side. There is only the template to update a bit.
It will be interesting to see how it compares in terms of performance with it's low-level equivalent or with its python/numba counterpart. We can probably also use it as a benchmark to test higher-level implementation, when we try a more general approach.

Pablo

Le mercredi 25 septembre 2013 20:19:54 UTC+2, Pablo Winant a écrit :

Hi,

Maybe some mildly ambitious project that we could start now, would be to create some good spline routines.

There are some existing libraries in Fortran (notably fitpack) and in C (einspline).
Python/Scipy has some wrappers around fitpack, and there has been some discussion on the scipy mailing list this spring to create a new implementation in Cython (a compiled dialect of python so to say), but it hasn't produced any output yet, that I am aware of. There is also a Julia package called BSplines, that is currently restricted to dimension 1.

Dolo already includes cubic splines with derivatives and extrapolation, on regular grids. These are generated from a template file that generate Cython code, which is very similar in essence to the einspline code, except that it can handle dimensions >3. Templating is used only for the evaluation, while traditional coding is used for the filtering step. It would be very easy to adapt the templating, to produce Julia code inead. It could be generated once for all, or replaced later by Julia macros. That would give us (natural) cubic splines in Julia with little work to do.

[Sébastien Villemot]

Le lundi 30 septembre 2013 à 14:46 -0700, Pablo Winant a écrit :
> I have made a tentative implementation for cubic splines in Julia
> here: https://github.com/EconForge/splines_julia . This is very
> preliminary, incomplete and inelegant, but since I won't be making
> much progress on it in the following weeks, I thought it was better to
> share it now.

Thanks for contributing this!

How would you want to proceed for the next steps? I see the following
possible improvements:

- improve the interface so that instead of having two functions in the
API (one returning an array of coefficients, the second for evaluating
the spline), we would have only one function in the API which would
return a function object (that we can call directly later to evaluate
the spline). This is the approach taken by the BSplines package and I
think it is much more julia'ish. It may also involve some performance
improvements since it some constants may be optimized out in the process
by the LLVM JIT (to be verified)

- add derivatives

- rewrite the macro system in Julia (and drop the python code)

- think about how to integrate this code into the BSplines package

- benchmark it

I am willing to give some help but I would need a paper/reference to
understand better what the code is doing.

Also I ran the test file and got an error for the 2D code:

2D interpolation

ERROR: invperm: input is not a permutation
in invperm at combinatorics.jl:110
in \ at linalg/factorization.jl:372
in \ at linalg/dense.jl:519
in filter_coeffs at /home/sebastien/splines_julia/splines_filter.jl:3
in include at boot.jl:238
in include_from_node1 at loading.jl:96
at /home/sebastien/splines_julia/test.jl:52

Best,

--
Sébastien Villemot
Researcher in Economics
Dynare developer
http://www.dynare.org/sebastien

[Pablo Winant]

Le vendredi 4 octobre 2013 17:58:46 UTC+2, Sébastien Villemot a écrit :

Le lundi 30 septembre 2013 à 14:46 -0700, Pablo Winant a écrit :
> I have made a tentative implementation for cubic splines in Julia
> here: https://github.com/EconForge/splines_julia . This is very
> preliminary, incomplete and inelegant, but since I won't be making
> much progress on it in the following weeks, I thought it was better to
> share it now.

Thanks for contributing this!

How would you want to proceed for the next steps? I see the following
possible improvements:

- improve the interface so that instead of having two functions in the
API (one returning an array of coefficients, the second for evaluating
the spline), we would have only one function in the API which would
return a function object (that we can call directly later to evaluate
the spline). This is the approach taken by the BSplines package and I
think it is much more julia'ish. It may also involve some performance
improvements since it some constants may be optimized out in the process
by the LLVM JIT (to be verified)

It's a good idea, but I don't know how you can implement that. Is it possible to attach some data to functions ?

 

- add derivatives

This I can do easily following the same approach as before. Ther is also the multivariate intepolation that can be added very easily. Then I would also like to have cartesian products of irregular grids, but there are some more formulas to change.
 

- rewrite the macro system in Julia (and drop the python code)

I agree. There is still one drawback here: the macro code will probably more convoluted that the generated one. Maybe it is still possible to keep it somewhere when calling the macro, in order to see the output.
Another thing: it is possible to avoid macros completely if we allow for recursive calls for each evaluation.

- think about how to integrate this code into the BSplines package

This require careful thinking. Currently, the BSpline package doesn't do interpolation at all, it just provides a way to evaluate the basis functions and the rest is left as a reader's exercise.
It would make sense to add to the bspline package some operations to filter the coeffs, and a function to evaluate splines. This would be a different approach, that would also happen to be more extensible.

There is also the Grid.jl package, that actually has interpolation routines and is pobably worth investigating. In terms of functionalities, it looks like it is currently more developed than the bsplines package, although it lacks extrapolation and natural BC. It should be possible to add it here. According to https://github.com/JuliaLang/julia/issues/2557, this package aims at performing simple gridded interpolation.

I think there is no complete consensus on how the interpolation/splines packages should be separated. At any rate, it is a good idea to start with a good api that suits our needs, to see where this api would fit more naturally in the Julia ecosystem (from a user perspective) and see were the computational code should be contributed.

- benchmark it

I'll do that. So far, I have seen that the Julia version as similar timings to the python/numba version and that both lags very slightly behind an equivalent C code. What will ultimately more interesting will be to compare the generated code to a more generic approach. Another interesting question would be to compare the vectorized evaluation (interpolate at each point) versus a repeated evaluation (one point after all),
 

I am willing to give some help but I would need a paper/reference to
understand better what the code is doing.

I'll try to document it soon and find some references. Basically, it uses blending formulas. In 1d, it gives an explicit cubic polynomial expression whose coefficients are determined by 4 values on the grid of the function to interpolate.
 

Also I ran the test file and got an error for the 2D code:

2D interpolation

ERROR: invperm: input is not a permutation
 in invperm at combinatorics.jl:110
 in \ at linalg/factorization.jl:372
 in \ at linalg/dense.jl:519
 in filter_coeffs at /home/sebastien/splines_julia/splines_filter.jl:3
 in include at boot.jl:238
 in include_from_node1 at loading.jl:96
at /home/sebastien/splines_julia/test.jl:52

Which system/juliia version do you use ? In principle the solution operator\ is supposed to be operated on sparse arrays, which is not consistent with the error message. What I don't understand is why you get the message only for 2d interpolation.
Are you still using julia 0.1 ? (https://groups.google.com/forum/#!searchin/julia-dev/deprecate$20julia$200.1/julia-dev/l4lZKoJcn94/W6QBcdxRZ7kJ)

Best,

Pablo
 

[Sébastien Villemot]

Yes, this is called "closures" in computer science, and Julia supports
them, as shown by this example:

julia> function f(x)
x = 1
g = y -> x + y
return(g)
end
# method added to generic function f

julia> f(1)(2)
3

> - rewrite the macro system in Julia (and drop the python
> code)
>
> I agree. There is still one drawback here: the macro code will
> probably more convoluted that the generated one. Maybe it is still
> possible to keep it somewhere when calling the macro, in order to see
> the output.
> Another thing: it is possible to avoid macros completely if we allow
> for recursive calls for each evaluation.

Which approach is preferable in your opinion?

I am willing to do the Julia macro stuff if we choose that path.

> Also I ran the test file and got an error for the 2D code:
>
> 2D interpolation
>
> ERROR: invperm: input is not a permutation
> in invperm at combinatorics.jl:110
> in \ at linalg/factorization.jl:372
> in \ at linalg/dense.jl:519
> in filter_coeffs
> at /home/sebastien/splines_julia/splines_filter.jl:3
> in include at boot.jl:238
> in include_from_node1 at loading.jl:96
> at /home/sebastien/splines_julia/test.jl:52
>
>
> Which system/juliia version do you use ? In principle the solution
> operator\ is supposed to be operated on sparse arrays, which is not
> consistent with the error message. What I don't understand is why you
> get the message only for 2d interpolation.
> Are you still using julia 0.1 ? (https://groups.google.com/forum/#!
> searchin/julia-dev/deprecate$20julia
> $200.1/julia-dev/l4lZKoJcn94/W6QBcdxRZ7kJ)

The test was done using a recent git snapshot (commit b3be803*
2013-10-04 13:15:40 UTC)

[Pablo Winant]

Yes of course. I don't know why I didn't think of it. That's very natural.

 

>         - rewrite the macro system in Julia (and drop the python
>         code)
>
> I agree. There is still one drawback here: the macro code will
> probably more convoluted that the generated one. Maybe it is still
> possible to keep it somewhere when calling the macro, in order to see
> the output.
> Another thing: it is possible to avoid macros completely if we allow
> for recursive calls for each evaluation.

Which approach is preferable in your opinion?

I am willing to do the Julia macro stuff if we choose that path.

I don't have a strong preference here. There is basically one single problematic line of code. It is a tensor reduction
of the form:

M_ijk x^i y^j z^k

where each index varies over 4 positions. Using a macro, one can expand this operation to perform it in one line. I believe it will always be more efficient than performing several reduction (the recursive way). 

So I think the macro is a good and simple way to achieve that and it would be great if you can do it. Maybe there is a way to avoid defining several functions for various dimensions, and have the specialization occur at compile-time when the dimension of arguments are determined).

Depending on how good julia can optimize across functions, we could isolate the tensor reduction as a function and compare the recursive definition of that function with a macro generated one.

>         Also I ran the test file and got an error for the 2D code:
>        
>         2D interpolation
>        
>         ERROR: invperm: input is not a permutation
>          in invperm at combinatorics.jl:110
>          in \ at linalg/factorization.jl:372
>          in \ at linalg/dense.jl:519
>          in filter_coeffs
>         at /home/sebastien/splines_julia/splines_filter.jl:3
>          in include at boot.jl:238
>          in include_from_node1 at loading.jl:96
>         at /home/sebastien/splines_julia/test.jl:52
>        
>
> Which system/juliia version do you use ? In principle the solution
> operator\ is supposed to be operated on sparse arrays, which is not
> consistent with the error message. What I don't understand is why you
> get the message only for 2d interpolation.
> Are you still using julia 0.1 ? (https://groups.google.com/forum/#!
> searchin/julia-dev/deprecate$20julia
> $200.1/julia-dev/l4lZKoJcn94/W6QBcdxRZ7kJ)

The test was done using a recent git snapshot (commit b3be803*
2013-10-04 13:15:40 UTC)

I just saw the bug. The error is triggered by the operation "(smax-smin)/(orders-1)" where smax,smin are supposed to be float vectors and orders is an int vectors.
The division here is supposed to be vectorized, hence it should be ./  instead of / .
This is a mistake that may come up quite often when code is ported from python to julia, since in python operations are vectorized by default. I just fixed it.

[Pablo Winant]

Some update on that front: I have added the computation of the gradient, and updated the api a bit. There are now many options to interpolate data that are demonstrated in this notebook: http://nbviewer.ipython.org/urls/raw.github.com/EconForge/splines_julia/master/cubic_splines_julia.ipynb . These options are:

• the low-level routines: there is a different one for each dimension, and all operations are inplace. These are generated with the template in csplines.jl
• an intermediate routine eval_UC_spline and eval_UC_spline_G (with the gradient). These call the relevant low-level function according to the dimension, an allocate the output.
• a high-level routine, that returns an interpolating function

There are still some pitfalls I don’t like about the implementation, but it would be good, if the high-level api could be tested a bit to see whether we like it or not.

Pablo