JuliaMath / Interpolations.jl

Fast, continuous interpolation of discrete datasets in Julia
http://juliamath.github.io/Interpolations.jl/
Other
531 stars 109 forks source link

Interpolant defined on irregularly spaced grid called on 3D matrix #161

Open shaizaq opened 7 years ago

shaizaq commented 7 years ago

I am trying to define a 9D interpolant on an irregularly spaced grid. In particular, I am defining the interpolant as follows:

itp = interpolate((x1, x2, x3, x4, x5, x6, x7, x8, x9), y, Gridded(Linear())); where x1, x2, x3, x7, x8 and x9 have 3 points each, while x4, x5 and x6 have 2 points each. This creates a 3 3 3 2 2 2 3 3 3 interpolant itp.

The problem is that I need to call itp on 5832 120 7 values for each of the 9 inputs, which is taking a lot of time and memory when I use for loops (over each dimension) to call the interpolant. In Matlab, I am able to pass the 3D array for each of the 9 inputs, and it returns a 5832 by 120 by 7 matrix in 1.4 seconds. What is a quick and easy way to do the same in Julia. Using for loops for each dimension took approximately 21 second in Julia.

floswald commented 7 years ago

Can you provide a minimal working example? Interpolations.jl performs very well in high dimensions if you use FixedSizeArrays.jl or StaticArrays.jl. There was an issue a while back with a response from Tim but I'm on a phone now so won't check. Let's see your example!

shaizaq commented 7 years ago

Currently, what I am doing is the following:

Define the interpolant: x1 = [0, 25000, 60000]; x2 = [-21000, 36000, 875000]; x3 = collect(1:1:3); x4 = collect(0:1:1); x5 = collect(0:1:1); x6 = collect(0:1:1); x7 = collect(1:1:no_types); x8 = [50, 300, 800]; x9 = linspace(2000,28000,3);

emax = rand(3,3,3,2,2,2,3,3,3);

itp = interpolate((x1, x2, x3, x4, x5, x6, x7, x8, x9), emax, Gridded(Linear()));

This generates a (3,3,3,2,2,2,3,3,3) interpolate.

I now need to call itp for 5832 x 120 x 7 = 4891320 values for each of the nine inputs. Currently, what I have done is the following:

Construct function foo which takes in itp and 9 input values as arguments and returns the interpolated value:

function foo(itp, x1, x2, x3, x4, x5, x6, x7, x8, x9) s = Array{eltype(x1)}(length(x1),1); for i=1:length( x1 ) s[i] = itp[x1[i],x2[i],x3[i],x4[i],x5[i],x6[i],x7[i],x8[i],x9[i]] end temp = reshape(s,5832,120,7); return temp end

Generate data to pass to function foo: a1 = rand(5832,120,7); b1 = reshape(a1,4891320); b2 = copy(b1); b3 = copy(b1); b4 = copy(b1); b5 = copy(b1); b6 = copy(b1); b7 = copy(b1); b8 = copy(b1); b9 = copy(b1);

interpolated_value = foo(itp,b1,b2,b3,b4,b5,b6,b7,b8,b9);

This returns a 5832 by 120 by 7 matrix of interpolated values, which is what I want. However, it takes more time to compute that than it does on Matlab (On my laptop, Matlab took 1.4 seconds while Julia took 3.3 seconds). What's a better way of doing the same thing? Is there an example using FixedSizeArrays.jl or StaticArrays.jl that I can refer to?

sglyon commented 7 years ago

Hi @shaizaq, what arguments are you passing to foo? Is it b1 through b9?

shaizaq commented 7 years ago

Yes, I am passing b1, through b9 as arguments to the function foo.