Renarks on the spline implementation

General:

The implementation follows the same approach as https://github.com/EconForge/python_interpolation , but uses Numba instead of Cython. It is also similar to the Julia version https://github.com/EconForge/python_interpolation (see the notebook http://nbviewer.ipython.org/github/EconForge/splines/blob/master/cubic_splines_julia.ipynb)

It is currently located in dolo.numeric.interpolation but the goal is to move the implementation to an independent package.

Structure

The filtering steps are implemented in filter_cubic_splines.py

The evaluation is made at three different levels:

• compiled code: eval_cubic_splines_numba.py, where actual computations are made

• low-level API: eval_cubic_splines.py, simple dimension independent API to evaluate splines. Optional inplace computations. Example:

    eval_cubic_spline(a, b, orders, coefs, points, values)
  

• high-level API: splines, interpolation object. Example:

    sp = MultivariateSpline(a, b, orders)
    sp.set_values(values)
    sp(points)
  

Filtering steps

The filtering step are done in file filter_cubic_splines.py and could probably be improved by:

• using a band-solver for the sake of clarity
• avoid recomputing the system to be solved many times (the banded matrix)
• find a way to write the filtering steps in a dimension-independent way
• filter multi-splines at once (maybe)
• parallelize the steps (how ?)

Code generation

The compiled code is generated using tempita, a simple template library for python. The splines functions are generated by: python gen_splines.py --multivariate=True --max_order=4 eval_cubic_splines_numba.py.in eval_cubic_splines_numba.py or python gen_splines.py --multivariate=True --max_order=4 eval_cubic_splines_cython.py.in eval_cubic_splines_cython.py.

Not all calculations can be made in the template because of this bug in tempita (https://bitbucket.org/ianb/tempita/issue/4/defining-and-using-functions-inside-py) that prevents the definition of recursive functions.

Alternatives to tempita:

• define a tensor reduction function and call it at the end of the function
• manipulate the AST, before compiling the function, maybe using https://github.com/lihaoyi/macropy

Numba/Cython

In the following remarks, "should" denotes an assertion that remains to be tested.

• There is no post-installation process for numba code, while an extension must be built for Cython (fails on some corporate machines)
• The performance cost of using Numba instead of Cython should be small, especially since LLVM now has a compiler.
• The LLVM compiler should optimize loops across functions. (i.e. for i in range(k): f(i)) We need to check that this is indeed the case and that we really get autovectorization in the current implementation (SIMD)
• Cython supports a parallel context manager that makes it easy to parallelize loops with openmp. With several processors, this leads to big performance gains.
• open question: what is the good way to jit-compile dimension independent functions ?

A good option would be to have both numba and cython implementation of the compiled routines and import only one of them for the low-level API.

Compilation of the Cython version

on pablo's machine gcc -pthread -fopenmp -fno-strict-aliasing -g -DNDEBUG -g -fwrapv -O3 -Wall -Wstrict-prototypes -fPIC -I/home/pablo/.local/opt/anaconda/include/python2.7 -c eval_cubic_splines_cython.c -o eval_cubic_splines_cython.o

then

gcc -pthread -shared eval_cubic_splines_cython.o -L/home/pablo/.local/opt/anaconda/lib -lpython2.7 -o eval_cubic_splines_cython.so

Open question: how to to incorporate the splines package as a submodule ? Add a setup.py file ?