[SciPy-Dev] scipy.interpolate - other use cases

Tue May 28 08:01:03 EDT 2019

Hi Sci-Py Devs,

I actually had this use case during my Ph.D. thesis, and have had variants
on it many times since.  Back then I needed interpolation between two sets
of 2D, regularly structured but skewed, grids.  The data changed but the
grid layout did not.  I also hit this one when computing with
semi-Lagrangian methods where I have a fixed advective term and time step;
not a particularly difficult case but one that chews up a lot of CPU in my
work.

I find saving the interpolant as an [ia,ja,a], (row index, column number,
value) sparse matrix works well, however, that means writing the whole
interpolant myself, since digging out coefficients is often as much work as
coding from scratch.  The matrix representation also means that a parabolic
solution step over the interpolation ends up looking like $L \cdot M^{-1}
\cdot K$ where L and K are interpolation matrices, so I get to use
carefully optimized matrix-vector operations once the interpolant is
computed.  I've done this for both regular and barycentric co-ordinates on
an underlying triangular/tetrahedral grid.  The resulting speed-up over
dynamically computing the interpolation weights is striking (can't tell you
exactly how much, off the top of my head).

This has been on my wish list for a while; I'll see if I can dig out some
code for it.

Best regards
-- Simon

On Tue, May 28, 2019 at 1:12 AM Dieter Werthmüller <Dieter at werthmuller.org>
wrote:

>  > Did you profile this? What's the relative time saved here for your
> use case?
>
> I did profile it. The speed-up factor in my tests is between 1.8 and 3,
> depending on problem size. Interestingly, for bigger problem sizes the
> factor is decreasing.
>
> - Factor 3 was for interpolating y-z slices from a 16x16x16 to a
> 32x32x32 grid by looping over x. In absolute values: 1.62 ms instead of
> 4.91 ms.
>
> - Factor 1.8 was for interpolating y-z slices from a 128x128x128 to a
> 256x256x256 grid by looping over x. In absolute values: 1.03 s instead
> of 1.67 s.
>
> I put the notebook here:
>
> https://github.com/prisae/tmp-share/blob/master/Time_RegGridProlongator_vs_scipyRegGridInterpolator.ipynb
>
>  > This hasn't come up before AFAIK, so doing this in scipy.interpolate
> may not be justified.
>
> Fair enough. It thought I'd ask. If it comes up again in the future we
> can still think about it. (I am working with a multigrid solver, so
> jumping a lot between different grid sizes, and these times add up, so
> it was worth exploring it for me.)
>
> Dieter
>
>
> On 27/05/2019 22:26, Ralf Gommers wrote:
> >
> >
> > On Fri, May 24, 2019 at 10:51 PM Dieter Werthmüller
> > <Dieter at werthmuller.org <mailto:Dieter at werthmuller.org>> wrote:
> >
> >     Dear Devs,
> >
> >     I have a particular issue with scipy.interpolate for which I would
> like
> >     to know:
> >
> >     - Is there already a better way to do it which I miss?
> >     - Is there interest to add this functionality to SciPy? or
> >     - Is there interest to extend the functionality of the existing fct?
> >     - None of the above or other ideas/comments?
> >
> >     My example deals with scipy.interpolate.RegularGridInterpolator;
> >     however, it probably is the same for many of the interpolators.
> >
> >     The usage of it in, e.g., 2D is
> >
> >           func = RegularGridInterpolator((x, y), data)
> >           func(pts)
> >
> >     where `x` and `y` define the regular grid on which `data` is defined,
> >     and we interpolate it at `pts`-coordinates. This is ideal if you
> >     want to
> >     interpolate the same data for different points again and again.
> >
> >     However, there are other use cases. Imagine you have a coarse grid,
> >     defined by `x`, `y`, and a fine grid, defined by the
> `pts`-coordinates.
> >     Both grids are regular grids, but not equidistant. Now we want to
> >     interpolate a changing `data` from the coarse grid to the fine grid.
> >     The
> >     current implementation is not good for that, as you have to create a
> >     new
> >     RegularGridInterpolator-instance for each new data-set, which
> involves
> >     the calculation of the indices and weights every time, even though
> they
> >     remain the same. So in this case, the implementation should instead
> >     look
> >     like:
> >
> >           func = RegularGridInterpolator((x, y), pts)
> >           func(data)
> >
> >     hence, `pts` and `data` interchanged. In this way, the required
> indices
> >     and weights can all be calculated during the initiation, which will
> >     significantly speed-up the interpolation. And `func` can
> >     subsequently be
> >     called relatively cheaply in, for instance, a loop which changes the
> >     data every time, but not the coordinates.
> >
> >
> > Did you profile this? What's the relative time saved here for your use
> case?
> >
> >
> >     This behaviour can be obtained by simply moving things between the
> >     `__init__` and `__call__` functions in `RegularGridInterpolator`,
> which
> >     could probably be done as a switch triggered by a new parameter.
> >
> >     Any thoughts?
> >
> >
> > This hasn't come up before AFAIK, so doing this in scipy.interpolate may
> > not be justified.
> >
> > Cheers,
> > Ralf
> >
> >
> >
> >
> > _______________________________________________
> > SciPy-Dev mailing list
> > SciPy-Dev at python.org
> > https://mail.python.org/mailman/listinfo/scipy-dev
> >
> _______________________________________________
> SciPy-Dev mailing list
> SciPy-Dev at python.org
> https://mail.python.org/mailman/listinfo/scipy-dev
>

-- 
1129 Ibbetson Lane
Mississauga, Ontario
L5C 1K9       Canada
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/scipy-dev/attachments/20190528/f8173d45/attachment-0001.html>