OffsetArrays.jl
OffsetArrays.jl copied to clipboard
collapse nested SubArray and OffsetArray
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?
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)
.
Hm, I see your point, but that's also, I'd say, the point of collapsing. OffsetArray
s nest in a collapsing way. SubArray
s nest in a collapsing way (in many cases, though not all).
The parent
property you mentioned is not preserved for nested view
s, 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
This looks like a similar case to https://github.com/JuliaLang/julia/pull/26872 so it might be good to have.
@timholy how do you think of this?
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