[SciPy-Dev] Multilinear interpolation
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 <email@example.com> wrote:
>> 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
>> - 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.
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
More information about the SciPy-Dev