Splines interpolation
Spline interpolation can be performed using two functions prefilter
and eval_splines
.
from interpolation.splines import prefilter, eval_splines
These two functions are jit compiled using numba, and should be fast, without big performance loss when used in a loop provided the loop is compiled.
The behaviour of these two functions is determined by the type of their argument. In the text below the types provided are numba types obtained using numba.typeof
function.
grids
Splines are defined by their values on a cartesian product of 1-dimensional grids.
Regular one dimensional grids, are represented by a Tuple(float64, float64, int64)
where first element is the lower bound, the second the upper bound and the latest the number of points.
An irregular one dimensional grid, is represented by a numpy array array(float64, 1d, C)
. It is assumed to be in increasing order. These are currently supported only for multilinear splines (k=1
)
A multidimensional grid, is represented as a tuple of one-dimensional grids. For instance, ((0.0, 1.0, 10), (0.0, 1.0, 10))
. Pay attention to one element tuples: (np.array([0.0, 0.1, 0.5]))
is not a correct grid (it is one-dimensional). It should be (np.array([0.0, 0.1, 0.5]),)
values
Given a d-dimensional grid of dimensions n_1 \times ... \times n_d. Values to be interpolated can be specified as:
- a n_1 \times ... \times n_d numpy array: scalar values
- a n_1 \times ... \times n_d \times n_x numpy array: vector values where n_x variables are defined on each node of the grid
prefiltering variables
For splines of order greater than 1, the coefficients used for interpolation must be prefiltered, for the resulting interpolant to inteporpolate exactly at grid points:
C = prefilter(G, V, k=3)
where G
is a grid and V
the values as described above.Currently, prefilter
is implemented only for k=3
and k=1
. In the k=1
it does nothing as prefiltering is not required.
The coefficient array is of size (n_1+k-1)\times ... \times (n_d+k-1) for scalar values and (n_1+k-1)\times ... \times (n_d+k-1) \times n_x for vector values.
Inplace calculations can be performed as:
prefilter(G, V, k=3, out=C)
Currently, prefiltering cubic splines, always used natural boundaries conditions (f''=0).
interpolating the function
To interpolate values:
eval_spline(G, C, P, out=None, k=3, diff="None", extrap_mode="linear")
where:
G
is a multi-dimensional grid as specified aboveC
an array of coefficientsP
denotes the locations at which to interpolate:- tuple or array of size
d
: point at which to interpolate - 2d array with column size
d
: list of points at which to interpolate
- tuple or array of size
out
:- if
None
, interpolated values are returned - if array: where to store the result inplace. The dimensions must be exactly equal to the array that would be returned by the function (see below)
- if
k: int
: spline order (currently (1 or 3))extrap_mode: str
: how to extrapolate outside of the grid. Either:- '"constant"': 0 outside of the grid
- '"nearest"': takes value from nearest point on the grid
- '"linear"' (default): projects to nearest points and use derivative at this point to extrapolate linearly
diff: str
: specifies which derivatives to compute- '"None"': no derivative
- string representing a tuple of tuple (see below)
just in time compilation and literal values
In the current eval_spline
specification, keyword arguments k
, extrap_mode
, and diff
, are used to control the generation of just in time code and must therefore be known at compile-time. They are ultimately treated as literal values (1 and 3 are of different types for instance.)
This currently implies a few limitations:
diff
must be passed as strings even though it represents a tuple of tuples- keyword arguments cannot be ommited. Error message is especially confusing, when they are.
- there might be a penalty cost in running this function outside of a numba jitted context. It is a known limitation in numba ~0.50, which will ultimately go away
specifying partial derivatives
Provided the interpolated function f(x) is defined on a d
dimensional space, with arguments $x_1, ... x_d$
. A partial derivative of any order (\partial^{k_1} ... \partial^{k_d}) is denoted by d-element tuple (k_1, ..., k_d)
.
The partial derivatives to be computed by eval_spline
can be specified as a tuple of tuples.
For instance, in a two dimensional space, to compute the value and the jacobian, one can pass:
( (0,0), (1,0), (0,1) )
Note that we don't pass the tuple directly to eval_spline
but a string. For the same example we would then pass:
"( (0,0), (1,0), (0,1) )"
dimensions of the output
Depending on the the type of arguments, the output (or the supplied array to perform inplace operations will have different dimensions).
Here is a summary. When relevant n_x
denotes the number of approximated variable (for vector valued interpolation), N
the number of points where the interpolant is approximated (for vectorized operations), n_j
the number of partial derivatives to evaluate.
Vectorized | Vector Valued | Derivatives | Output |
---|---|---|---|
no | no | no | float (no inplace) |
no | n_x |
no | n_x array |
N |
no | no | N array |
N |
n_x |
no | N.n_x array |
no | no | n_j |
n_j tuple (no inplace) |
no | n_x |
n_j |
n_x.n_j array |
N |
no | n_j |
N . n_j array |
N |
n_x |
n_j |
N.n_x.n_j array |