Feature request: Interpolate one array, "re-use" weights on other arrays

[eirikbrandsaas]

In my code I need to interpolate many different arrays at the same points. A signficant part of my computation time is spent on interpolating, but I can cut the time down significantly if it would be possible to "reuse" interpolation weights and positions. The following example illustrates.

something not to different from this:

## Create arrays, interpolation objects etc
xgrid = 1:1:10
ygrid = 1:1:10
array1 = [xgrid[ix] + 1*ygrid[iy] for ix=1:length(xgrid), iy = 1:length(ygrid) ]
array2 = [xgrid[ix] + 2*ygrid[iy] for ix=1:length(xgrid), iy = 1:length(ygrid) ]
array3 = [xgrid[ix] + 3*ygrid[iy] for ix=1:length(xgrid), iy = 1:length(ygrid) ] # Three arrays defined over the same grids
array1_intrp = LinearInterpolation((xgrid,ygrid),array1,extrapolation_bc=Line())
array2_intrp = LinearInterpolation((xgrid,ygrid),array2,extrapolation_bc=Line())
array3_intrp = LinearInterpolation((xgrid,ygrid),array3,extrapolation_bc=Line())  # Three interpolations objects

## Do lots of calculations that require calling the interpolation objects many times

for i = 1:1000 # big loop
  i = 1
  x = 10.0/i # Some values which depends on which iteration I am on
  y = 10.0/i
  val1 = array1_intrp(x,y)
  val2 = array2_intrp(x,y)
  val3 = array3_intrp(x,y)
  val = function(val1,val2,val3) # some function that combines the three interpolated values
end

Now, the point is that since all of the interpolations are evaluated at the same values (x,y) and are "defined" on the same grids, we don't really need to do all the calculations to find val2, val3, since we can reuse the weights that were used in finding val1.

Then the question becomes:

1. Is this currently possibly? If yes, then
   1.1 Is it possible to make interpolations.jl export the weights and closest grid points so that I can manually perfrom the interpolations from val2,val3?
   1.2 Is there a way to tell interpolations to just interpolate three objects at the same time?
2. Is this something you would consider adding in the future?

[timholy]

Yes, it's possible, see the WeightedIndex framework. There's an example in RegisterHindsight, specifically https://github.com/HolyLab/RegisterHindsight.jl/blob/803b42cbaec4c271fb7d743bcb8fe502c1582d56/src/RegisterHindsight.jl#L71-L113. It would be great to have someone contribute this to the docs in this package.

[eirikbrandsaas]

Wow, that's amazing. I'll take a look at it in the coming days, and if I'm able to get it working I'll try to create a small example that could be used in the documentation.

[eirikbrandsaas]

So I tried playing around with the WeightedIndex stuff which worked great*, but I still can't figure out how to make the interpolations package return the i and f (as they are called in the documentation in WeightedIndex) I need if I give it an array and a point I want the array valued at.** I tried looking at your code, but frankly, there was too much going on for me to decipher.

*: At least until you give it an index i that is out of bounds because then it gives you bogus values without errors.
**: Of course I could find those myself, but I assume there is a way to make interpolations return those values.

[timholy]

Try Juno.@enter array1_intrp(1.3, 1.8) and see what Interpolations itself does. You should pretty quickly find your way to

Interpolations.jl/src/b-splines/indexing.jl

Lines 5 to 9 in 44ad9a6

	@inline function (itp::BSplineInterpolation{T,N})(x::Vararg{Number,N}) where {T,N}
	@boundscheck (checkbounds(Bool, itp, x...) || Base.throw_boundserror(itp, x))
	wis = weightedindexes((value_weights,), itpinfo(itp)..., x)
	itp.coefs[wis...]
	end

[timholy]

At least until you give it an index i that is out of bounds because then it gives you bogus values without errors.

Bounds checking happens before WeightedIndex computation (yes, would be good to document this too).

[timholy]

See #365. The preview of the new devdocs should help a lot. Good luck!

[pvillacorta]

Hi, did you get anywhere with this? I need to solve (almost) exactly the same problem as the one illustrated in the first comment.
Thank you :)

[eirikbrandsaas]

No, I never made any progress (see #484 )

[Balinus]

This feature would be awesome for weather and climate data. For example, I often have thousands of arrays that needs to be interpoolated on a new grid and the source and destination grid is static.

[mkitti]

I'll have some time in about a month to think about this a bit more. There's probably a way to hack thisnow, but a stable interface to do this would be better.