julia - 在域外正确使用 Interpolations.jl

Question

我正在将 Matlab 代码移植到 julia 中，到目前为止我得到了惊人的结果:一个在 Matlab 中运行 5 个多小时的代码，julia 只用了 8 分钟多一点!但是我有一个问题......在 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

只要 interpolar_a 中的点超出 NivelAltura_int 的范围，就会产生 NaN。

在 Julia 中，我正在尝试对以下内容做同样的事情:

for xx in 1:xlong
    for yy in 1:ylong
        AltInterp = interpolate((Znw_r,),A_s_c_r,Gridded(Linear()));
        NivelAltura_int[xx,yy,1:end] = AltInterp[Znu[1:end]]
        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()));
        U_alturas[xx,yy,1:end] = Uinterp[Alturas[1:end]];
        V_alturas[xx,yy,1:end] = Vinterp[Alturas[1:end]];
    end
end

使用包 Interpolations.jl。每当该点在域外时，此程序包都会进行推断，这对我的目的来说是不正确的。我可以添加几行代码来检查域外的值并将其替换为 NaN，但我相信这会增加一些计算时间并且不是很优雅。

在包的文档中，提到了这样一种对象:

        Uextrap = extrapolate(Uinterp,NaN)

为了控制域外的行为，但我还没有找到如何使用它，我试过在 Uinterp 下添加它，我试过评估它，但它自然不会那样工作。

你能帮我解决这个问题吗？

谢谢!

Answer 1

看来您可能遇到了两个问题。首先，最近有一些关于网格外推 (#101) 的工作可能还没有在标记版本中。如果你愿意生活在边缘，你可以Pkg.checkout("Interpolations")使用开发版(Pkg.free("Interpolations")会让你再次回到稳定版本)。

其次，看起来仍然缺少向量值网格外推法 (issue #24):

julia> using Interpolations
       itp = interpolate((collect(1:10),), collect(.1:.1:1.), Gridded(Linear()))
       etp = extrapolate(itp, NaN);

julia> etp[.5:1:10.5]
ERROR: BoundsError: # ...
 in throw_boundserror at abstractarray.jl:156
 in getindex at abstractarray.jl:488

如您所见，它试图对所有抽象数组使用通用定义，这当然会引发边界错误。插值只需要添加自己的定义。

与此同时，您可以使用带有标量索引的理解:

julia> [etp[x] for x=.5:1:10.5]
11-element Array{Any,1}:
 NaN
   0.15
   0.25
   0.35
   0.45
   0.55
   0.65
   0.75
   0.85
   0.95
 NaN