Add one and copy for Diagonal

[putianyi889]

#50
JuliaArrays/FillArrays.jl#126

[dlfivefifty]

This counts as "type piracy": it changes behaviour of FillArrays.jl which is bad.

Instead restrict to the infinite case:

one(D::Diagonal{T,<:AbstractFill{T,Tuple{OneToInf{Int}}}) where T = Eye{T}(size(D,1))
copy(D::Diagonal{T,<:AbstractFill{T,Tuple{OneToInf{Int}}}) where T = D

Though this could be refined to allow for infinite offset indexing.

[dlfivefifty]

Add tests for one and copy directly as well.

[dlfivefifty]

Add tests that Eye(∞) == Eye(∞)^2

[dlfivefifty]

Add tests for Diagonal(Fill(2,∞))

[putianyi889]

^2 shouldn't be a problem because LinearAlgebra deals with matrix powers like this:

# base\intfuncs.jl:243
function power_by_squaring(x_, p::Integer)
    x = to_power_type(x_)
    if p == 1
        return copy(x)
    elseif p == 0
        return one(x)
    elseif p == 2
        return x*x
    elseif p < 0
        isone(x) && return copy(x)
        isone(-x) && return iseven(p) ? one(x) : copy(x)
        throw_domerr_powbysq(x, p)
    end
    t = trailing_zeros(p) + 1
    p >>= t
    while (t -= 1) > 0
        x *= x
    end
    y = x
    while p > 0
        t = trailing_zeros(p) + 1
        p >>= t
        while (t -= 1) >= 0
            x *= x
        end
        y *= x
    end
    return y
end

I added the test anyway.

[putianyi889]

Actually I encountered this problem when doing M^0 where M is infinite and banded. Shall we fix this problem in a more general way?

Tridiagonal has the same issues.

[dlfivefifty]

The right place to do this is in ArrayLayouts.jl. It's not possible to do this for every array but we can do it for LayoutArray, which is the super type for BandedMatrix, BlockArray, etc. and enables adding features to arrays that aren't supported for base. So we could add to ArrayLayouts.jl something like

_one(_, _, d::AbstractMatrix{T}) where T = Base.invoke(one, Tuple{AbstractMatrix{T}}, d) # standard def
function _one(::DiagonalLayout, _, d::AbstractMatrix{T}) where T
    m,n = size(d)
     m==n || throw(DimensionMismatch("multiplicative identity defined only for square matrices"))
    Eye{T}(m)
end
# ... Tridiagonal, etc., special cases
one(d::LayoutMatrix) = _one(MemoryLayout(d), axes(d), d)
one(d::SubArray{<:Any,2,<:LayoutMatrix}) = _one(MemoryLayout(d), axes(d), d)
one(d::Diagonal{<:Any,<:LayoutMatrix}) = _one(MemoryLayout(d), axes(d), d)
one(d::Tridiagonal{<:Any,<:LayoutVector}) = _one(MemoryLayout(d), axes(d), d)
one(d::SymTridiagonal{<:Any,<:LayoutVector}) = _one(MemoryLayout(d), axes(d), d)

Then in BandedMatrices.jl:

function _one(::AbstractBandedLayout, _, d::AbstractMatrix{T}) where T
    m,n = size(d)
     m==n || throw(DimensionMismatch("multiplicative identity defined only for square matrices"))
    convert(BandedMatrix, Eye{T}(m))
end

[codecov]

Codecov Report

Merging #52 (2dddc9d) into master (a9437da) will increase coverage by 0.07%.
The diff coverage is 100.00%.

@@            Coverage Diff             @@
##           master      #52      +/-   ##
==========================================
+ Coverage   77.44%   77.51%   +0.07%     
==========================================
  Files           4        4              
  Lines         625      627       +2     
==========================================
+ Hits          484      486       +2     
  Misses        141      141              

Impacted Files	Coverage Δ
src/infarrays.jl	76.28% <100.00%> (+0.30%)	⬆️

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update a9437da...2dddc9d. Read the comment docs.