JuliaArrays / OffsetArrays.jl

Fortran-like arrays with arbitrary, zero or negative starting indices.
Other
195 stars 40 forks source link

Implement a default offset option for all higher dimensions, to support a more complete `centered` #340

Closed BioTurboNick closed 8 months ago

BioTurboNick commented 9 months ago

Fixes #339 , but requires support from Julia Base to call the right axes methods. The API for this is probably not right, and obviously needs tests and docs. But would like feedback first.

jishnub commented 9 months ago

One concern that I see is that a lot of code uses Base.require_one_based_indexing to shield against offset axes, but this does not check the axes beyond the ndims of an array. Perhaps we need better support from Base for offset trailing axes. IIUC, this PR would mean that the transpose of a one-based OffsetVector would end up being offset if the default offset is non-zero.

BioTurboNick commented 9 months ago

IUC, this PR would mean that the transpose of a one-based OffsetVector would end up being offset if the default offset is non-zero.

Yes, exactly.

One concern that I see is that a lot of code uses Base.require_one_based_indexing to shield against offset axes, but this does not check the axes beyond the ndims of an array.

Though couldn't OffsetArrays provide an implementation of that?

BioTurboNick commented 9 months ago

require_one_based_indexing relies on has_offset_axes, which I provided an implementation of here.

BioTurboNick commented 9 months ago

Seems I have a StackOverflow issue in one of the constructors. Would it be better as a kwarg?

jishnub commented 9 months ago

I'm a bit unsure of this PR's implications, e.g.

julia> v = OffsetVector(1:2, (0,), 3)
1:2 with indices 1:2

julia> v'
1×2 adjoint(OffsetArray(::UnitRange{Int64}, 1:2)) with eltype Int64 with indices 4:4×1:2:
 1  2

julia> view(v,:)'
1×2 adjoint(view(OffsetArray(::UnitRange{Int64}, 1:2), :)) with eltype Int64 with indices Base.OneTo(1)×OffsetArrays.IdOffsetRange(values=1:2, indices=1:2):
 1  2

julia> v == view(v,:)
true

julia> v' == view(v,:)'
false

so wrappers will behave differently as they are not inheriting the trailing axes of the parent

BioTurboNick commented 9 months ago

Hmm, so it's using the axes(A::AbstractArray{T,N}, d) method, which returns hardcoded OneTo(1) for trailing dimensions.

So there could be an axes(AA::SubArray, d) that passes through to the parent.

Although views of OffsetArrays seem to behave inconsisently?

julia> v = OffsetVector(1:2, (-1,))
1:2 with indices 0:1

julia> axes(view(v, :))
(Base.IdentityUnitRange(OffsetArrays.IdOffsetRange(values=0:1, indices=0:1)),)

julia> axes(view(v, 0:1))
(Base.OneTo(2),)

Why does : preserve parent indexing while specifying the range resets them to 1:n?

jishnub commented 9 months ago

This is because the axes of an indexing operation are the axes of the indices, and a : is transformed to a Base.Slice internally, which has the same indices as its values. On the other hand, the axes of a UnitRange are 1-based, and are unrelated to the starting point of the range. This is one of the fundamental (and less intuitive) aspects of offset indexing, and I had this same question answered by Tim Holy three years back

BioTurboNick commented 9 months ago

Oy. Ok.

Given how tangled this is... maybe the easiest solution is to take defensive measures in the ImageFiltering functions.

jishnub commented 8 months ago

Closing, as this doesn't seem like something that we'll do here.