JuliaArrays / OffsetArrays.jl

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

Why is Diagonal(::OffsetArray) not supported? #112

Open dlfivefifty opened 4 years ago

dlfivefifty commented 4 years ago

At the moment Diagonal(::OffsetArray) is explicitly forbidden:

julia> using OffsetArrays

julia> x = OffsetArray(randn(5),-1)
5-element OffsetArray(::Array{Float64,1}, 0:4) with eltype Float64 with indices 0:4:
 -1.1571843743724437
  1.182586195865172
 -1.7939602910206534
 -0.42805091122795136
 -1.2399140299181175

julia> Diagonal(x)
ERROR: ArgumentError: offset arrays are not supported but got an array with index other than 1
Stacktrace:
 [1] require_one_based_indexing at ./abstractarray.jl:89 [inlined]
 [2] Diagonal at /Users/solver/Projects/julia-1.4/usr/share/julia/stdlib/v1.4/LinearAlgebra/src/diagonal.jl:9 [inlined]
 [3] Diagonal(::OffsetArray{Float64,1,Array{Float64,1}}) at /Users/solver/Projects/julia-1.4/usr/share/julia/stdlib/v1.4/LinearAlgebra/src/diagonal.jl:13
 [4] top-level scope at REPL[3]:1

However, supporting it seems to require just 2 changes in Base:

julia> Base.require_one_based_indexing(::OffsetArray) = nothing
julia> Base.axes(D::Diagonal) = (axes(D.diag,1),axes(D.diag,1))

julia> Diagonal(x)
5×5 Diagonal{Float64,OffsetArray{Float64,1,Array{Float64,1}}} with indices 0:4×0:4:
 -1.15718   ⋅         ⋅         ⋅          ⋅ 
   ⋅       1.18259    ⋅         ⋅          ⋅ 
   ⋅        ⋅       -1.79396    ⋅          ⋅ 
   ⋅        ⋅         ⋅       -0.428051    ⋅ 
   ⋅        ⋅         ⋅         ⋅        -1.23991

julia> Diagonal(x)[0,0]
-1.1571843743724437

Is there any reason this was not done?

One reason I ask is that the definition Base.axes(D::Diagonal) = (axes(D.diag,1),axes(D.diag,1)) is helpful for BlockArrays.jl, so that Diagonal(::BlockArray) automatically inherits the block structure of the diagonal. I was thinking of making a PR into Base to propose this change but wanted to know context why it wasn't done already.

timholy commented 4 years ago

Is there any reason this was not done?

Probably just because doing it requires writing good tests, etc., and until someone needs this enough to bite the bullet it's much better to have it error. require_one_based_indexing was added in the closing days of the Julia 1.0 merge window as a way to protect against more serious bugs.

dlfivefifty commented 4 years ago

I see, I'll consider doing a PR when I get the time.

Is there a way to bypass an internal constructor, that is, to make a Diagonal(::OffsetArray) without wrongly overloading Base.require_one_based_indexing?

timholy commented 4 years ago

It's in the inner constructor precisely to make that "not possible." (In quotes because of course, it actually is possible, but I ain't tellin' :smile:. Ordinary Julia code should not be calling builtins directly.)