collapse nested SubArray and OffsetArray

[goretkin]

import OffsetArrays: OffsetArray

a = rand(4,4)
a[:] .= 1:16

sa = view(a, 1:3, 1:3)
osa = OffsetArray(sa, (11:13, 11:13))
sosa = view(osa, 11:12, 11:12)

function Base.view(osa::OffsetArray{T, N, <:SubArray}, I::Vararg{Any,N}) where {T, N}
    view(parent(osa), map((a,b) -> a .- b, I, osa.offsets)...)
end

sosa2 = view(osa, 11:12, 11:12)

julia> sosa
2×2 view(OffsetArray(view(::Array{Float64,2}, 1:3, 1:3), 11:13, 11:13), 11:12, 11:12) with eltype Float64:
 1.0  5.0
 2.0  6.0

julia> sosa2
2×2 view(::Array{Float64,2}, 1:2, 1:2) with eltype Float64:
 1.0  5.0
 2.0  6.0

The map and splat is probably not great for performance. Other than that, does this seem like a reasonable idea?

[johnnychen94]

Doesn't look good to me since parent(sosa) != parent(sosa2) and I believe the behavior of parent(soda) is the correct one; I mean, it doesn't look good to me to strip the offsets when we call parent(::SubArray).

[goretkin]

Hm, I see your point, but that's also, I'd say, the point of collapsing. OffsetArrays nest in a collapsing way. SubArrays nest in a collapsing way (in many cases, though not all).

The parent property you mentioned is not preserved for nested views, e.g.

ssa = SubArray(sa, (1:3, 1:3))
ssa2= view(sa, 1:3, 1:3)

julia> ssa
3×3 view(view(::Array{Float64,2}, 1:3, 1:3), 1:3, 1:3) with eltype Float64:
 1.0  5.0   9.0
 2.0  6.0  10.0
 3.0  7.0  11.0

julia> ssa2
3×3 view(::Array{Float64,2}, 1:3, 1:3) with eltype Float64:
 1.0  5.0   9.0
 2.0  6.0  10.0
 3.0  7.0  11.0

julia> parent(ssa)
3×3 view(::Array{Float64,2}, 1:3, 1:3) with eltype Float64:
 1.0  5.0   9.0
 2.0  6.0  10.0
 3.0  7.0  11.0

julia> parent(ssa2)
4×4 Array{Float64,2}:
 1.0  5.0   9.0  13.0
 2.0  6.0  10.0  14.0
 3.0  7.0  11.0  15.0
 4.0  8.0  12.0  16.0

[johnnychen94]

This looks like a similar case to JuliaLang/julia#26872 so it might be good to have.

@timholy how do you think of this?

[jishnub]

This does have some unexpected consequences

julia> sosa3 = @view osa[Base.IdentityUnitRange(11:12), Base.IdentityUnitRange(11:12)]
2×2 view(::Array{Float64,2}, 1:2, 1:2) with eltype Float64:
 1.0  5.0
 2.0  6.0

julia> axes(sosa3)
(Base.OneTo(2), Base.OneTo(2))

This breaks the rule that if b = a[idx] for some vector indices idx, then b[j] = a[idx[j]]. This also means that slices and views currently have different axes if the indices are not 1-based.

julia> axes(osa[Base.IdentityUnitRange(11:12), Base.IdentityUnitRange(11:12)])
(11:12, 11:12)

julia> axes(@view osa[Base.IdentityUnitRange(11:12), Base.IdentityUnitRange(11:12)])
(Base.OneTo(2), Base.OneTo(2))

This needs some more thinking