reshape may accept AbstractUnitRanges that begin at 1 as arguments

[jishnub]

Currently this works:

julia> reshape(1:2, 1, Base.OneTo(2))
1×2 reshape(::UnitRange{Int64}, 1, 2) with eltype Int64:
 1  2

but this doesn't

julia> reshape(1:2, 1, 1:2)
ERROR: MethodError: no method matching reshape(::UnitRange{Int64}, ::Tuple{Int64, UnitRange{Int64}})
Closest candidates are:
  reshape(::AbstractVector, ::Colon) at reshapedarray.jl:115
  reshape(::AbstractArray, ::Int64...) at reshapedarray.jl:116
  reshape(::AbstractArray, ::Union{Int64, AbstractUnitRange}...) at reshapedarray.jl:110
  ...
Stacktrace:
 [1] reshape(::UnitRange{Int64}, ::Int64, ::UnitRange{Int64})
   @ Base ./reshapedarray.jl:110
 [2] top-level scope
   @ REPL[2]:1

After this PR the latter will work:

julia> reshape(1:2, 1, 1:2)
1×2 reshape(::UnitRange{Int64}, 1, 2) with eltype Int64:
 1  2

Performance seems identical:

julia> @btime reshape($(Ref(1:2))[], 1, :, Base.OneTo(2));
  36.971 ns (0 allocations: 0 bytes)

julia> @btime reshape($(Ref(1:2))[], 1, :, 1:2);
  36.971 ns (0 allocations: 0 bytes)

Aside from this, these now behave correctly:

julia> reshape(1:2, 1, :, Base.OneTo(2))
1×1×2 reshape(::UnitRange{Int64}, 1, 1, 2) with eltype Int64:
[:, :, 1] =
 1

[:, :, 2] =
 2

julia> reshape(1:2, 1, Int32(2))
1×2 reshape(::UnitRange{Int64}, 1, 2) with eltype Int64:
 1  2

These added methods to reshape also ensure that permutedims works correctly for OffsetVectors:

julia> permutedims(OffsetArray(1:2, 3:4))
1×2 OffsetArray(reshape(::UnitRange{Int64}, 1, 2), 1:1, 3:4) with eltype Int64 with indices 1:1×3:4:
 1  2

[jishnub]

This uses a general Tuple instead of one with a union of valid types to allow OffsetArrays to define methods for AbstractUnitRange arguments that return OffsetArrays

[mcabbott]

This means that loading OffsetArrays will change some results. That's a stronger form of piracy than changing e.g. zeros(1:5) from an error into valid input. Are there other places it already does this?

julia> reshape('a':'h', 1:2, 1:4)  # with this PR
2×4 reshape(::StepRange{Char, Int64}, 2, 4) with eltype Char:
 'a'  'c'  'e'  'g'
 'b'  'd'  'f'  'h'

julia> using OffsetArrays

julia> reshape('a':'h', 1:2, 1:4)
2×4 OffsetArray(reshape(::StepRange{Char, Int64}, 2, 4), 1:2, 1:4) with eltype Char with indices 1:2×1:4:
 'a'  'c'  'e'  'g'
 'b'  'd'  'f'  'h'

[jishnub]

Yes loading OffsetArrays will indeed change the type of the array. This is unfortunate but I don't know if there is a type-stable way to resolve this? I think reshape is where this conflict is the worst, while such piracy similar and other functions have mostly been resolved. Tim's attitude in OffsetArrays has often been that types are an implementation detail, and it's better to think of arrays as ordered (key=>value) collections. The piracy in OffsetArrays is in the same spirit as it does not change the axes or values of the arrays, only the type. I don't fully subscribe to this view as the types often dictate performance, and would be happy to see this resolved. Perhaps a trait might be used to distinguish ranges that are known to start at 1 at compile time from those that aren't, and Base may define methods for the former while OffsetArrays may define methods for the latter? This was the goal of Base.OneTo, except the addition of types such as StaticArrays' SOneTo complicates the scene.

The goal of this PR was mainly to allow custom 1-based range types such as the axes of StaticArrays to be used in a reshape. On nightly this doesn't work:

julia> s = SVector{2}(1:2);

julia> reshape(s, axes(s,1), 1)
ERROR: MethodError: no method matching reshape(::SVector{2, Int64}, ::Tuple{SOneTo{2}, Int64})

With this PR this will work:

julia> s = SVector{2}(1:2);

julia> reshape(s, axes(s,1), 1)
2×1 reshape(::SVector{2, Int64}, 2, 1) with eltype Int64:
 1
 2

The axes of StaticArrays should work identically to Base.OneTo aside from the static size information.The PR also changes the behavior for UnitRanges that start at 1 as a side-effect.

[mcabbott]

unfortunate but I don't know if there is a type-stable way to resolve this? ... similar and other functions have mostly been resolved

I thought this was resolved precisely by not defining such methods: similar(ones(2), 1:3) is an error because similar(ones(2), 2:4) must be an OffsetArray. But I don't know if that's an official stance, or just something I assumed from trying things!

For StaticArrays, isn't deciding what to do with its own static range types something it ought to handle? In many cases it does so:

julia> using StaticArrays

julia> s = SA[1,2];

julia> m = similar(1:0, axes(s,1), axes(s,1))  # correctly overloaded
2×2 MMatrix{2, 2, Int64, 4} with indices SOneTo(2)×SOneTo(2):
  25  4607885968
 100          26

julia> vec(m)  # correctly overloaded -- this is preferable to reshape(::MMatrix{...})
4-element MVector{4, Int64} with indices SOneTo(4):
         25
        100
 4607885968
         26

julia> axes(reshape(m, :)) # not overloaded, hence loses static information
(Base.OneTo(4),)

julia> reshape(s, axes(s,1), axes(s,2))  # also not overloaded, that's a bug perhaps?
ERROR: MethodError: no method matching reshape(::SVector{2, Int64}, ::Tuple{SOneTo{2}, Base.OneTo{Int64}})

julia> using OffsetArrays

julia> reshape(s, axes(s,1), axes(s,2))  # piracy! But at least it was an error before
2×1 OffsetArray(reshape(::SVector{2, Int64}, 2, 1), 1:2, 1:1) with eltype Int64 with indices 1:2×1:1:
 1
 2

julia> axes(ans)  # static structure is not preserved
(OffsetArrays.IdOffsetRange(values=1:2, indices=1:2), OffsetArrays.IdOffsetRange(values=1:1, indices=1:1))

julia> similar(1:0, axes(s,1))  # this is not affected by piracy, good
2-element MVector{2, Int64} with indices SOneTo(2):
 4631963440
 5544374880

Maybe there are stranger examples than mine. You can still use dispatch to catch (Base.OneTo(1), Base.OneTo(2), SOneTo(3), Base.OneTo(4)) with a union, if you know that Base will catch the all-OneTo case. (In my example above, making axes(s,2) === SOneTo(1) would avoid such mixtures.)

Perhaps a trait might be used to distinguish ranges that are known to start at 1 at compile time from those that aren't, and Base may define methods for the former

This sounds like an argument for an AbstractOneTo, which both Base.OneTo and SOneTo could subtype. But the list of functions any package must overload when defining a new axis type (like reshape for SOneTo above) does not seem so very long.

[jishnub]

So let me try to summarize the status:

1. permutedims does not preserve the axes of OffsetVectors, as permutedims for an AbstractVector is a reshape, and reshape in Base is not axis-aware
2. Changing permutedims to use the axis in reshape will break current code using permutedims(::StaticArray) as StaticArrays don't define methods for reshape that accepts a mix of Integers and SOneTo. This was not required till now (and the reshape operation led to an error in any case)
3. Adding methods to reshape in Base that accept ranges will get around this. While this can not be done in general, this may be done for ranges that start at 1. This fixes the issue for StaticArrays
4. However adding methods to reshape would exacerbate the type-piracy in OffsetArrays, as it relies on these methods not being available in Base.
5. Ideally StaticArrays should define methods for reshape that accept SOneTo and we don't have to add methods to Base. However since this doesn't exist currently, there is no way to change permutedims to accept axes without breaking StaticArrays. I wonder if such breakages due to missing methods are acceptable?

[jishnub]

Incidentally reshape(s, axes(s,1), axes(s,2)) should be fixed by JuliaArrays/StaticArrays.jl#916

[mcabbott]

I would say that permutedims(v::AbstractVector) = reshape(v, 1, length(v)) not using axes is just an oversight which should be fixed, as I see you do here. Sorry I missed this at first. Whether the new dimension should be 1 or Base.OneTo(1) or axes(v,2) perhaps needs discussion. (I didn't know you were allowed to mix sizes & axes in reshape.)

Fixing that will mean that it just works on a 1D OffsetArray, since OffsetArrays already overloads all relevant reshape methods. Although arguably OffsetArrays should pass this command though to the array it wraps, in case that can do something better than reshape -- for example permutedims(OffsetArray(Diagonal(1:3),4,5)) now gives an OffsetArray(::SparseMatrixCSC) which is a bit ugly.

Fixing that will also turn permutedims(SA[1,2,3]) into an error, or an OffsetArray(reshape(::SVector)) if this is loaded. But the problem is just that this isn't implemented yet, but would be easy to fix. There was discussion about how to express the permutation statically for higher dimensions, but for vectors & matrices this doesn't matter.

julia> permutedims(SA[1 2; 3 4])  # not implemented, goes via similar, fails to preserve SMatrix
2×2 MMatrix{2, 2, Int64, 4} with indices SOneTo(2)×SOneTo(2):
 1  3
 2  4

julia> adjoint(SA[1 2; 3 4])  # is implemented, preserves SMatrix
2×2 SMatrix{2, 2, Int64, 4} with indices SOneTo(2)×SOneTo(2):
 1  3
 2  4

StaticArrays chooses to leave all reshapes which aren't fully static as errors, e.g. reshape(SA[1 2; 3 4], SOneTo(2), :). I suppose it could catch all those with only one non-static factor. Perhaps nobody has needed that.

[jishnub]

permutedims isn't overloaded for Diagonal anyway, as most array types resort to similar to generate the container.

julia> permutedims(Diagonal(1:4))
4×4 SparseArrays.SparseMatrixCSC{Int64, Int64} with 4 stored entries:
 1  ⋅  ⋅  ⋅
 ⋅  2  ⋅  ⋅
 ⋅  ⋅  3  ⋅
 ⋅  ⋅  ⋅  4

The OffsetArray wrapper doesn't change much in this particular example.

OffsetVectors may pass on the permutedims operation to the parent (and they probably should), but this creates an issue downstream, as any other array type wrapping an OffsetVector will also need to overload permutedims to avoid the generic fallback and lose the offset. Ideally this won't be necessary and the fallback method would handle offset axes correctly.

This PR only changes the permutedims method for an AbstractVector as this is the only one that relies on a reshape and loses the axis information. Arrays of other dimensionalities call similar which preserves axes. The docstring for permutedims states that the 1D case is a reshape (which seems to indicate that the underlying data is shared? Maybe it uses reshape as a generic term).

[mcabbott]

My unwrapping example requires 1.7, it seems.

Fixing permutedims seems fine. But I don't see the argument that this requires defining reshape(rand(4), 1:2, 1:2), which I tried to argue changes the terms of OffsetArrays's licensed piracy.

[jishnub]

I agree wrt the reshape and am happy to drop that if we may fix permutedims without adding methods to reshape.
I see this behavior currently:

julia> A = ones(2:3)
2-element OffsetArray(::Vector{Float64}, 2:3) with eltype Float64 with indices 2:3:
 1.0
 1.0

julia> S = StructArray{ComplexF64}((A, A))
2-element StructArray(OffsetArray(::Vector{Float64}, 2:3), OffsetArray(::Vector{Float64}, 2:3)) with eltype ComplexF64 with indices 2:3:
 1.0 + 1.0im
 1.0 + 1.0im

julia> T = ComplexF64[1+im for _ in axes(A,1)]
2-element OffsetArray(::Vector{ComplexF64}, 2:3) with eltype ComplexF64 with indices 2:3:
 1.0 + 1.0im
 1.0 + 1.0im

julia> S == T
true

julia> permutedims(S) == permutedims(T)
false

julia> permutedims(S)
1×2 StructArray(::Matrix{Float64}, ::Matrix{Float64}) with eltype ComplexF64:
 1.0+1.0im  1.0+1.0im

julia> permutedims(T)
1×2 OffsetArray(::Matrix{ComplexF64}, 1:1, 2:3) with eltype ComplexF64 with indices 1:1×2:3:
 1.0+1.0im  1.0+1.0im

Unless the permutedims in Base is made axis-aware, every wrapper type will have to add methods to deal with offset axes.