[SciPy-Dev] Multilinear interpolation

Pablo Winant pablo.winant@gmail....
Thu Feb 21 08:37:27 CST 2013

On 21/02/2013 15:14, Robert Kern wrote:
> 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?

Yes exactly.

> Here is the implementation that I keep using, which (ab)uses
> ndimage.map_coordinates():
> 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.
> https://github.com/scipy/scipy/blob/master/scipy/ndimage/src/ni_interpolation.c#L349
I hesitate between two options that seem be to be used somewhere in 
scipy : Mako templates and code generation via a simple python module. 
Do you have a preference ?
> --
> Robert Kern
> _______________________________________________
> SciPy-Dev mailing list
> SciPy-Dev@scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-dev

More information about the SciPy-Dev mailing list