# Optimized interpolation routines in Python / numba

The library contains:

`splines.*`

: fast numba-compatible multilinear and cubic interpolation`multilinear.*`

: fast numba-compatible multilinear interpolation (alternative implementation)`smolyak.*`

: smolyak polynomials- complete polynomials

## install

Install latest version:

- from conda:
`conda install -c conda-forge interpolation`

- from PyPI:
`pip install interpolation`

Latest development version from git:

```
pip install poetry # if not already installed
pip install git+https://github.com/econforge/interpolation.py.git/
```

The project uses a `pyproject.toml`

file instead of `setup.py`

and other legacy configuration files. For those used to development installation, this is feasible using `dephell`

:

```
pip install dephell # if not already installed
dephell --from=pyproject.toml --to=setup.py # only once
pip install -e . # like old times
```

## multilinear and cubic interpolation

Fast numba-accelerated interpolation routines for multilinear and cubic interpolation, with any number of dimensions. Several interfaces are provided.

### eval_linear

Preferred interface for multilinear interpolation. It can interpolate on uniform and nonuniform cartesian grids. Several extrapolation options are available.

```
import numpy as np
from interpolation.splines import UCGrid, CGrid, nodes
# we interpolate function
f = lambda x,y: np.sin(np.sqrt(x**2+y**2+0.00001))/np.sqrt(x**2+y**2+0.00001)
# uniform cartesian grid
grid = UCGrid((-1.0, 1.0, 10), (-1.0, 1.0, 10))
# get grid points
gp = nodes(grid) # 100x2 matrix
# compute values on grid points
values = f(gp[:,0], gp[:,1]).reshape((10,10))
from interpolation.splines import eval_linear
# interpolate at one point
point = np.array([0.1,0.45]) # 1d array
val = eval_linear(grid, values, point) # float
# interpolate at many points:
points = np.random.random((10000,2))
eval_linear(grid, values, points) # 10000 vector
# output can be preallocated
out = np.zeros(10000)
eval_linear(grid, values, points, out) # 10000 vector
## jitted, non-uniform multilinear interpolation
# default calls extrapolate data by using the nearest value inside the grid
# other extrapolation options can be chosen among NEAREST, LINEAR, CONSTANT
from interpolation.splines import extrap_options as xto
points = np.random.random((100,2))*3-1
eval_linear(grid, values, points, xto.NEAREST) # 100
eval_linear(grid, values, points, xto.LINEAR) # 10000 vector
eval_linear(grid, values, points, xto.CONSTANT) # 10000 vector
# one can also approximate on nonuniform cartesian grids
grid = CGrid(np.linspace(-1,1,100)**3, (-1.0, 1.0, 10))
points = np.random.random((10000,2))
eval_linear(grid, values, points) # 10000 vector
# it is also possible to interpolate vector-valued functions in the following way
f = lambda x,y: np.sin(x**3+y**2+0.00001)/np.sqrt(x**2+y**2+0.00001)
g = lambda x,y: np.sin(x**3+y**2+0.00001)/np.sqrt(x**2+y**2+0.00001)
grid = UCGrid((-1.0, 1.0, 10), (-1.0, 1.0, 10))
gp = nodes(grid) # 100x2 matrix
mvalues = np.concatenate([
f(gp[:,0], gp[:,1]).reshape((10,10))[:,:,None],
g(gp[:,0], gp[:,1]).reshape((10,10))[:,:,None]
],axis=2) # 10x10x2 array
points = np.random.random((1000,2))
eval_linear(grid, mvalues, points[:,1]) # 2 elements vector
eval_linear(grid, mvalues, points) # 1000x2 matrix
out = np.zeros((1000,2))
eval_linear(grid, mvalues, points, out) # 1000x2 matrix
# finally, the same syntax can be used to interpolate using cubic splines
# one just needs to prefilter the coefficients first
# the same set of options apply but nonuniform grids are not supported (yet)
f = lambda x,y: np.sin(x**3+y**2+0.00001)/np.sqrt(x**2+y**2+0.00001)
grid = UCGrid((-1.0, 1.0, 10), (-1.0, 1.0, 10))
gp = nodes(grid) # 100x2 matrix
values = f(gp[:,0], gp[:,1]).reshape((10,10))
# filter values
from interpolation.splines import filter_cubic
coeffs = filter_cubic(grid, values) # a 12x12 array
from interpolation.splines import eval_cubic
points = np.random.random((1000,2))
eval_cubic(grid, coeffs, points[:,1]) # 2 elements vector
eval_cubic(grid, coeffs, points) # 1000x2 matrix
out = np.zeros((1000,2))
eval_cubic(grid, coeffs, points, out) # 1000x2 matrix
```

*Remark*: the arguably strange syntax for the extapolation option comes from the fact the actualy method called must be determined by type inference. So `eval_linear(..., extrap_method='linear')`

would not work because the type of the last argument is always a string. Instead, we use opts.CONSTANT and opts.LINEAR for instance which have different numba types.

Despite what it looks `UCGrid`

and `CGRid`

are not objects but functions which return very simple python structures that is a tuple of its arguments. For instance, `((0.0,1.0, 10), (0.0,1.0,20))`

represents a 2d square discretized with 10 points along the first dimension and 20 along the second dimension. Similarly `(np.array([0.0, 0.1, 0.3, 1.0]), (0.0, 1.0, 20))`

represents a square nonuniformly discretized along the first dimension (with 3 points) but uniformly along the second one. Now type dispatch is very sensitive to the exact types (floats vs ints), (tuple vs lists) which is potentially error-prone. Eventually, the functions `UCGrid`

and `CGrid`

will provide some type check and sensible conversions where it applies. This may change when if a parameterized structure-like object appear in numba.

### interp

Simpler interface. Mimmicks default `scipy.interp`

: mutlilinear interpolation with constant extrapolation.

```
### 1d grid
from interpolation import interp
x = np.linspace(0,1,100)**2 # non-uniform points
y = np.linspace(0,1,100) # values
# interpolate at one point:
interp(x,y,0.5)
# or at many points:
u = np.linspace(0,1,1000) # points
interp(x,y,u)
```

### object interface

This is for compatibility purpose only, until a new jittable model object is found.

```
from interpolation.splines import LinearSpline, CubicSpline
a = np.array([0.0,0.0,0.0]) # lower boundaries
b = np.array([1.0,1.0,1.0]) # upper boundaries
orders = np.array([50,50,50]) # 50 points along each dimension
values = np.random.random(orders) # values at each node of the grid
S = np.random.random((10**6,3)) # coordinates at which to evaluate the splines
# multilinear
lin = LinearSpline(a,b,orders,values)
V = lin(S)
# cubic
spline = CubicSpline(a,b,orders,values) # filter the coefficients
V = spline(S) # interpolates -> (100000,) array
```

### development notes

Old, unfair timings: (from `misc/speed_comparison.py`

)

```
# interpolate 10^6 points on a 50x50x50 grid.
Cubic: 0.11488723754882812
Linear: 0.03426337242126465
Scipy (linear): 0.6502540111541748
```

More details are available as an example notebook (outdated)

Missing but available soon: - splines at any order - derivative

Feasible (some experiments) - evaluation on the GPU (with numba.cuda) - parallel evaluation (with guvectorize)

## smolyak interpolation

See testfile for examples.