Type instability in hessian when interpolating a 3D field

[touste]

Hi, see the following MWE:

using Interpolations
nx = 5
A = rand(Float64, nx, nx, nx) * 100
itp = interpolate(A, BSpline(Quadratic(Reflect(OnCell()))))
@code_warntype Interpolations.hessian(itp, 1., 1., 1.)

This gives a type instability on the wis variable:
wis::Tuple{Tuple{Interpolations.WeightedAdjIndex{3, Float64}, Interpolations.WeightedAdjIndex{3, Float64}, Interpolations.WeightedAdjIndex{3, Float64}}, Tuple{Interpolations.WeightedAdjIndex{3, Float64}, Interpolations.WeightedAdjIndex{3, Float64}, Interpolations.WeightedAdjIndex{3, Float64}}, Tuple{Interpolations.WeightedAdjIndex{3, Float64}, Interpolations.WeightedAdjIndex{3, Float64}, Interpolations.WeightedAdjIndex{3, Float64}}, Tuple{Interpolations.WeightedAdjIndex{3, Float64}, Interpolations.WeightedAdjIndex{3, Float64}, Interpolations.WeightedAdjIndex{3, Float64}}, Tuple{Interpolations.WeightedAdjIndex{3, Float64}, Interpolations.WeightedAdjIndex{3, Float64}, Interpolations.WeightedAdjIndex{3, Float64}}, Vararg{Any}}

This in turn leads to a return type Any for the function.

[mkitti]

See JuliaLang/julia#43368

[touste]

Are you sure this is related? I see the same issue in 1.6.4 and 1.5.4

[touste]

That was fast! Thank you very much for fixing this.

[N5N3]

Well, the above fix seems only work for 3d case. Some more effort is needed.
Has been extended to 7d, should be enough for most cases.