[SciPy-Dev] Multilinear interpolation

Robert Kern robert.kern@gmail....
Thu Feb 21 08:14:00 CST 2013

On Thu, Feb 21, 2013 at 12:08 PM, Pablo Winant <pablo.winant@gmail.com> wrote:
> Hi,
> Continuing a series of thread about interpolation, I have an almost
> ready implementation of multilinear implementation, that is linear
> interpolation in each dimension.  I have been toying with several
> implementations
> - in pure numpy (quite memory hungry)
> - in cython
> - using pycuda
> Obviously, if there is some interest for it in scipy, it would be a good
> motivation to finish it.

Yes, please! I have had to rewrite such an implementation a couple of
times already. Just to be clear, by "multilinear" interpolation, do
you mean that you are given a set of values on an orthogonal, but
otherwise irregular N-dimensional grid? Then on the interior of each
grid cell, you do linear interpolation for each axis in turn using the
grid points of the cell, generalizing bilinear interpolation?

Here is the implementation that I keep using, which (ab)uses

import numpy as np
from scipy import ndimage

def nlinear(x, yp, axes):
    coords = []
    for i in range(x.shape[-1]):
        coords.append(np.interp(x[..., i], axes[i], np.arange(len(axes[i]))))
    coords = np.array(coords)
    return ndimage.map_coordinates(yp, coords, order=1)

> There are two main differences with respect to
> the existing interp routine :
> - it is independent of fitpack
> - it extrapolates linearly while the existing routine returns a constant
> (or NaN)
> Now I have some questions about it:
> - in this case I several implementations, in the same library that can
> potentially be enabled on different systems (with a compiler, with a
> gpu). I made a python file that tries to import one version after all
> shadowing slower ones by faster ones. It is meant to replicate the
> .m/.mex file behavior of Matlab. Is there a better approach ?

I think for scipy, we would just use the plain Cython implementation.
Everyone who installs scipy from source needs a compiler anyways. I
personally don't like GPU versions being used implicitly just because
a GPU is available. I suspect we just wouldn't include the GPU version
in scipy.

> - the cython code is meant to operate on double* inputs, and is limited
> to 4 dimensions. It is however very easy to generate single or higher
> order code. What would you use for that purpose ?

You can take a look at how ndimage does it. It's not pretty.


Robert Kern

More information about the SciPy-Dev mailing list