Remarks on the spline implementation

General:

This package was initially a copy of dolo.numeric.interpolation. It is also similar to the Julia version https://github.com/EconForge/splines (see the notebook http://nbviewer.ipython.org/github/EconForge/splines/blob/master/cubic_splines_julia.ipynb)

It currently contains:

• cubic splines
• multilinear interpolation (splines of order 1).

Structure

The filtering steps are implemented in filter_cubic_splines.py

The evaluation is made at three different levels:

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

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

• low-level API: eval_cubic_splines.py, simple dimension independent API to evaluate splines. Optional inplace computations. The following, evaluates a spline on several points.

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

• compiled code: eval_cubic_splines_numba.py or eval_cubic_splines_numba.py where actual computations are made. There are Numba and Cython versions. We currently keep both of them, because the numba version seems a bit slower. We keep numba's version because cython does not work on all systems.

Filtering steps

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

• implement other termination conditions
• 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 ?)
• have several versions of the filter, probably a numba and a cython version (understand why numba's version seems to be much slower

Note: there maybe some inspiration to take from https://github.com/timholy/Grid.jl/blob/master/src/interp.jl , for instance the use of Woodbury matrices.

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

The Cython files can be compiled by performing:

python setup.py build_ext --inplace

This command outputs, the whole Gcc command, which can also be issued manually. For instance, on pablo's machine, the shared libraries can be compiled with:

• cython eval_cubic_splines_cython.pyx
• 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
• 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 of another package (quantecon or olo) ? Do we keep the setup file ? Is there any problem with nested submodules ? (for instance if splines is a submodule of dolo, which is also included as a submodule in quantecon)