Optimized interpolation routines in Python / numba
The library contains:

splines.*
: fast numbacompatible multilinear and cubic interpolation 
multilinear.*
: fast numbacompatible multilinear interpolation (alternative implementation) 
smolyak.*
: smolyak polynomials  complete polynomials
install
Install latest version:
 from conda:
conda install c condaforge 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 numbaaccelerated 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, nonuniform 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))*31
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 vectorvalued 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 errorprone. Eventually, the functions UCGrid
and CGrid
will provide some type check and sensible conversions where it applies. This may change when if a parameterized structurelike 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 # nonuniform 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.