[Numpy-discussion] Interpolation question

Anne Archibald peridot.faceted at gmail.com
Sun Mar 28 04:26:28 EDT 2010

On 27 March 2010 20:24, Andrea Gavana <andrea.gavana at gmail.com> wrote:
> Hi All,
>    I have an interpolation problem and I am having some difficulties
> in tackling it. I hope I can explain myself clearly enough.
> Basically, I have a whole bunch of 3D fluid flow simulations (close to
> 1000), and they are a result of different combinations of parameters.
> I was planning to use the Radial Basis Functions in scipy, but for the
> moment let's assume, to simplify things, that I am dealing only with
> one parameter (x). In 1000 simulations, this parameter x has 1000
> values, obviously. The problem is, the outcome of every single
> simulation is a vector of oil production over time (let's say 40
> values per simulation, one per year), and I would like to be able to
> interpolate my x parameter (1000 values) against all the simulations
> (1000x40) and get an approximating function that, given another x
> parameter (of size 1x1) will give me back an interpolated production
> profile (of size 1x40).

If I understand your problem correctly, you have a function taking one
value as input (or one 3D vector) and returning a vector of length 40.
You want to know whether there are tools in scipy to support this.

I'll say first that it's not strictly necessary for there to be: you
could always just build 40 different interpolators, one for each
component of the output. After all, there's no interaction in the
calculations between the output coordinates. This is of course
awkward, in that you'd like to just call F(x) and get back a vector of
length 40, but that can be remedied by writing a short wrapper
function that simply calls all 40 interpolators.

A problem that may be more serious is that the python loop over the 40
interpolators can be slow, while a C implementation would give
vector-valued results rapidly. To make this work, you're going to have
to find a compiled-code interpolator that returns vector values. This
is not in principle complicated, since they just need to run the same
interpolator code on 40 sets of coefficients. But I don't think many
of scipy's interpolators do this. The only ones I'm sure are able to
do this are the polynomial interpolators I wrote, which work only on
univariate inputs (and provide no kind of smoothing). If you're going
to use these I recommend using scipy's spline functions to construct
smoothing splines, which you'd then convert to a piecewise cubic.

But I'd say, give the 40 interpolators a try. If they're slow, try
interpolating many points at once: rather than giving just one x
value, call the interpolators with a thousand (or however many you
need) at a time. This will speed up scipy's interpolators, and it will
make the overhead of a forty-element loop negligible.


> Something along these lines:
> import numpy as np
> from scipy.interpolate import Rbf
> # x.shape = (1000, 1)
> # y.shape = (1000, 40)
> rbf = Rbf(x, y)
> # New result with xi.shape = (1, 1) --> fi.shape = (1, 40)
> fi = rbf(xi)
> Does anyone have a suggestion on how I could implement this? Sorry if
> it sounds confused... Please feel free to correct any wrong
> assumptions I have made, or to propose other approaches if you think
> RBFs are not suitable for this kind of problems.
> Thank you in advance for your suggestions.
> Andrea.
> "Imagination Is The Only Weapon In The War Against Reality."
> http://xoomer.alice.it/infinity77/
> ==> Never *EVER* use RemovalGroup for your house removal. You'll
> regret it forever.
> http://thedoomedcity.blogspot.com/2010/03/removal-group-nightmare.html <==
> _______________________________________________
> NumPy-Discussion mailing list
> NumPy-Discussion at scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion

More information about the NumPy-Discussion mailing list