Interpolation of 3D Arrays does not work since v0.9.0

[schustermartin]

Since the release of v0.9.0, which introduced Bounds Errors, the interpolation of 3d grids does not work anymore. The following code is a minimum example extended from the 2d example code. It runs just fine in v0.8.0 but not in later releases:

A = rand(8,8,20)
knots = ([x^2 for x = 1:8], [x^2 for x = 1:8], [0.2y for y = 1:20])
itp = interpolate(knots, A, Gridded(Linear()))
itp(4,4,1.2)
A[2,2,6]

At line 3 where I initialize the interpolations object I get the following error message

Error showing value of type Interpolations.GriddedInterpolation{Float64,3,Float64,Gridded{Linear},Tuple{Array{Int64,1},Array{Int64,1},Array{Float64,1}}}: ERROR: BoundsError: attempt to access 8×8×20 interpolate((::Array{Int64,1},::Array{Int64,1},::Array{Float64,1}), ::Array{Float64,3}, Gridded(Linear())) with element type Float64 at index [Base.OneTo(8), Base.OneTo(8), 18]

[timholy]

Put a semicolon after itp = interpolate(knots, A, Gridded(Linear())) to suppress printing---it's just the printing of itp causing this. This will be fixed by #227 and/or #242.

[timholy]

I'm beginning to think we may need to create a custom show method for AbstractInterpolation to circumvent this bug, just to hold us over during the deprecation period. Thoughts?

[natschil]

I think it's a good idea. The current behaviour is confusing and annoying. Creating a show() method does a lot to reduce that confusion, and has no real disadvantages as far as I can tell.

[fratrik]

The test in test/io.jl for this issue passes for me.

However, I still run into the same problem in the REPL.

julia> include("io.jl")
Test Summary: | Pass  Total
IO            |   17     17
Test.DefaultTestSet("IO", Any[DefaultTestSet("BSpline", Any[], 6, false), DefaultTestSet("Gridded", Any[], 4, false), DefaultTestSet("scaled", Any[], 3, false), DefaultTestSet("Monotonic", Any[], 1, false), DefaultTestSet("Extrapolation", Any[], 2, false), DefaultTestSet("Combinations", Any[], 1, false), DefaultTestSet("Issue #260", Any[], 0, false)], 0, false)

julia> show(scale(interpolate(zeros(3,3,3),BSpline(Cubic(Free(OnGrid())))),2:4,2:4,2:4))
ERROR: BoundsError: attempt to access 3×3×3 scale(interpolate(OffsetArray(::Array{Float64,3}, 0:4, 0:4, 0:4), BSpline(Cubic(Free(OnGrid())))), (2:4, 2:4, 2:4)) with element type Float64 at index [Base.OneTo(3), Base.OneTo(3), 1]
Stacktrace:
 [1] throw_boundserror(::ScaledInterpolation{Float64,3,Interpolations.BSplineInterpolation{Float64,3,OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}},BSpline{Cubic{Free{OnGrid}}},Tuple{Base.OneTo{Int64},Base.OneTo{Int64},Base.OneTo{Int64}}},BSpline{Cubic{Free{OnGrid}}},Tuple{UnitRange{Int64},UnitRange{Int64},UnitRange{Int64}}}, ::Tuple{Base.OneTo{Int64},Base.OneTo{Int64},Int64}) at ./abstractarray.jl:484
 [2] checkbounds at ./abstractarray.jl:449 [inlined]
 [3] view(::ScaledInterpolation{Float64,3,Interpolations.BSplineInterpolation{Float64,3,OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}},BSpline{Cubic{Free{OnGrid}}},Tuple{Base.OneTo{Int64},Base.OneTo{Int64},Base.OneTo{Int64}}},BSpline{Cubic{Free{OnGrid}}},Tuple{UnitRange{Int64},UnitRange{Int64},UnitRange{Int64}}}, ::Base.OneTo{Int64}, ::Base.OneTo{Int64}, ::Int64) at ./subarray.jl:147
 [4] show_nd(::IOContext{Base.TTY}, ::ScaledInterpolation{Float64,3,Interpolations.BSplineInterpolation{Float64,3,OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}},BSpline{Cubic{Free{OnGrid}}},Tuple{Base.OneTo{Int64},Base.OneTo{Int64},Base.OneTo{Int64}}},BSpline{Cubic{Free{OnGrid}}},Tuple{UnitRange{Int64},UnitRange{Int64},UnitRange{Int64}}}, ::getfield(Base, Symbol("##364#365")){String}, ::Bool) at ./arrayshow.jl:293
 [5] _show_nonempty at ./arrayshow.jl:403 [inlined]
 [6] show(::Base.TTY, ::ScaledInterpolation{Float64,3,Interpolations.BSplineInterpolation{Float64,3,OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}},BSpline{Cubic{Free{OnGrid}}},Tuple{Base.OneTo{Int64},Base.OneTo{Int64},Base.OneTo{Int64}}},BSpline{Cubic{Free{OnGrid}}},Tuple{UnitRange{Int64},UnitRange{Int64},UnitRange{Int64}}}) at ./arrayshow.jl:421
 [7] show(::ScaledInterpolation{Float64,3,Interpolations.BSplineInterpolation{Float64,3,OffsetArrays.OffsetArray{Float64,3,Array{Float64,3}},BSpline{Cubic{Free{OnGrid}}},Tuple{Base.OneTo{Int64},Base.OneTo{Int64},Base.OneTo{Int64}}},BSpline{Cubic{Free{OnGrid}}},Tuple{UnitRange{Int64},UnitRange{Int64},UnitRange{Int64}}}) at ./show.jl:313
 [8] top-level scope at none:0

[dkarrasch]

It's due to the fact that there is no show method for ScaledInterpolation objects. Without the scaling, everything works fine:
There is a show method for ScaledInterpolation objects, but these

Interpolations.jl/src/io.jl

Lines 115 to 117 in 56a4569

	Base.show(io::IO, ::MIME{Symbol("text/plain")}, X::ScaledInterpolation) = show_ranged(io, X, getknots(X))
	Base.show(io::IO, ::MIME{Symbol("text/plain")}, X::GriddedInterpolation) = show_ranged(io, X, getknots(X))
	Base.show(io::IO, ::MIME{Symbol("text/plain")}, X::AbstractExtrapolation) = show_ranged(io, X, getknots(X))

are not reached when calling show(sipt). The pipeline is show(sitp) -> show(IOBuffer(), sitp) and then directly to arrayprint or whatever. It seems that adding

Base.show(io::IO, X::ScaledInterpolation)   = show_ranged(io, X, getknots(X))

and similar for the other two would fix the problem. Here is my local output:

show(scale(interpolate(zeros(3,3,3),BSpline(Cubic(Free(OnGrid())))),2:4,2:4,2:4))
3×3×3 scale(interpolate(OffsetArray(::Array{Float64,3}, 0:4, 0:4, 0:4), BSpline(Cubic(Free(OnGrid())))), (2:4, 2:4, 2:4)) with element type Float64:
[:, :, 1] =
 0.0  0.0  0.0
 0.0  0.0  0.0
 0.0  0.0  0.0

[:, :, 2] =
 0.0  0.0  0.0
 0.0  0.0  0.0
 0.0  0.0  0.0

[:, :, 3] =
 0.0  0.0  0.0
 0.0  0.0  0.0
 0.0  0.0  0.0

Needless to say that I have no clue when the 3-args show gets called and how.