Using Interpolations.jl How to define the behavior when out of bounds?

[argel....@ciencias.unam.mx]

Hi!
I'm porting a Matlab code into Julia, so far my results are amazing:

The program in Matlab takes 5-6 hours to run while the program in Julia takes a little more than 8 minutes! However, there is a difference in the results...

In julia i'm using the Interpolations package to do the same thing i did with interp1 in matlab:

  for xx=1:xlong
            for yy = 1:ylong
                U_alturas(xx,yy,:) = interp1(squeeze(NivelAltura_int(xx,yy,:)),squeeze(U(xx,yy,:)), interpolar_a);
                V_alturas(xx,yy,:) = interp1(squeeze(NivelAltura_int(xx,yy,:)),squeeze(V(xx,yy,:)), interpolar_a);
            end
        end.

And in julia:

            Uinterp = interpolate((squeeze(NivelAltura_int[xx,yy,1:end],(1,2)),),squeeze(U[xx,yy,1:end],(1,2)),Gridded(Linear()));
            Vinterp = interpolate((squeeze(NivelAltura_int[xx,yy,1:end],(1,2)),),squeeze(V[xx,yy,1:end],(1,2)),Gridded(Linear()));
            Uextrap = extrapolate(Uinterp,NaN);
            Vextrap = extrapolate(Vinterp,NaN);
            U_alturas[xx,yy,1:end] = Uinterp[Alturas[1:end]];
            V_alturas[xx,yy,1:end] = Vinterp[Alturas[1:end]];

Where the tu lines in the middle are intended to produce NaN if the point is out of the domain.

The matter is that Matlab, for some values of xx,yy, produces NaN, while Julia produces numbers, which is not correct. What could i  be missing?

If the two middle lines are the correct way to establish the behavior for points outside the domain, maybe the difference comes from another reason.

Can you help me?

Thanks!

[Tomas Lycken]

Yes, that's the correct way to specify the extrapolation behavior, but you don't seem to use it afterwards; your last two lines refer to Uinterp and Vinterp, rather than Uextrap and Vextrap. To get the desired behavior for OOB indices, you must index into the result of extrapolate(...).

If you want to avoid this, you can of course do both interpolation and extrapolation in one step:

itp = extrapolate(interpolate(...), NaN)
itp[-3.5] # NaN

// T

[argel....@ciencias.unam.mx]

Thanks!

[Carlo Galli]

Thanks! Very helpful.

Looks like however evaluating the extrapolation object below on a Array ( itp = extrapolate(interpolate(...), NaN)[-3.5,-3.4] for example ) throws an error, while this works well with interpolation objects.

A solution would be to use a comprehension
[extrapolate(interpolate(...), NaN)[x] for x in Array]

but for some reason this becomes a type Any array, so it requires
Array{Float64}([extrapolate(interpolate(...), NaN)[x] for x in Array])

Have I missed something, or this is just the way extrapolate works? Apologies if this has already been posted somewhere, am quite new to Julia and sometimes I find that documentation seems not too easy to access.

Thanks again!

[Tomas Lycken]

Sorry for not getting back to you sooner on this.

Indexing into an extrapolation object like that should work - if you haven't already solved that problem, could you please file an issue with a complete (runnable) code sample?

Regarding the comprehension, did you do that in the repl or in a script? Does it still yield Any[] if you wrap it in a function? If so, it might indicate a type stability issue in extrapolate - please file an issue for that too then.

Thanks!

// T