Support inv on diagonal matrices from unitranges

[DanielVandH]

Fixes #176

[jishnub]

Won't this change the behavior from throwing a SingularException to returning an Inf if one of the values is zero?

[DanielVandH]

Technically yes, although this inv never worked in the first place for nonsingular matrices. Similar methods for inv don't throw, e.g.

julia> A = Diagonal((0:∞) .+ (0:∞)) # force non-infrange type
ℵ₀×ℵ₀ Diagonal{Int64, LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{InfiniteArrays.InfUnitRange{Int64}, InfiniteArrays.InfUnitRange{Int64}}}} with indices OneToInf()×OneToInf():
 0  ⋅  ⋅  ⋅  ⋅   ⋅   ⋅   ⋅   ⋅   ⋅   ⋅   ⋅   ⋅  ⋅  ⋅  ⋅  …
 ⋅  2  ⋅  ⋅  ⋅   ⋅   ⋅   ⋅   ⋅   ⋅   ⋅   ⋅   ⋅  ⋅  ⋅  ⋅
 ⋅  ⋅  4  ⋅  ⋅   ⋅   ⋅   ⋅   ⋅   ⋅   ⋅   ⋅   ⋅  ⋅  ⋅  ⋅
 ⋅  ⋅  ⋅  6  ⋅   ⋅   ⋅   ⋅   ⋅   ⋅   ⋅   ⋅   ⋅  ⋅  ⋅  ⋅
 ⋅  ⋅  ⋅  ⋅  8   ⋅   ⋅   ⋅   ⋅   ⋅   ⋅   ⋅   ⋅  ⋅  ⋅  ⋅
 ⋅  ⋅  ⋅  ⋅  ⋅  10   ⋅   ⋅   ⋅   ⋅   ⋅   ⋅   ⋅  ⋅  ⋅  ⋅  …
 ⋅  ⋅  ⋅  ⋅  ⋅   ⋅  12   ⋅   ⋅   ⋅   ⋅   ⋅   ⋅  ⋅  ⋅  ⋅
 ⋅  ⋅  ⋅  ⋅  ⋅   ⋅   ⋅  14   ⋅   ⋅   ⋅   ⋅   ⋅  ⋅  ⋅  ⋅
 ⋅  ⋅  ⋅  ⋅  ⋅   ⋅   ⋅   ⋅  16   ⋅   ⋅   ⋅   ⋅  ⋅  ⋅  ⋅
 ⋅  ⋅  ⋅  ⋅  ⋅   ⋅   ⋅   ⋅   ⋅  18   ⋅   ⋅   ⋅  ⋅  ⋅  ⋅
 ⋅  ⋅  ⋅  ⋅  ⋅   ⋅   ⋅   ⋅   ⋅   ⋅  20   ⋅   ⋅  ⋅  ⋅  ⋅  …
 ⋅  ⋅  ⋅  ⋅  ⋅   ⋅   ⋅   ⋅   ⋅   ⋅   ⋅  22   ⋅  ⋅  ⋅  ⋅
 ⋅  ⋅  ⋅  ⋅  ⋅   ⋅   ⋅   ⋅   ⋅   ⋅   ⋅   ⋅  24  ⋅  ⋅  ⋅
 ⋮               ⋮                   ⋮                ⋮  ⋱

julia> inv(A)
ℵ₀×ℵ₀ Diagonal{Float64, LazyArrays.BroadcastVector{Float64, typeof(inv), Tuple{LazyArrays.BroadcastVector{Int64, typeof(+), Tuple{InfiniteArrays.InfUnitRange{Int64}, InfiniteArrays.InfUnitRange{Int64}}}}}} with indices OneToInf()×OneToInf():
 Inf    ⋅    ⋅     ⋅         ⋅      ⋅    ⋅         …
   ⋅   0.5   ⋅     ⋅         ⋅      ⋅    ⋅
   ⋅    ⋅   0.25   ⋅         ⋅      ⋅    ⋅
   ⋅    ⋅    ⋅    0.166667   ⋅      ⋅    ⋅
   ⋅    ⋅    ⋅     ⋅        0.125   ⋅    ⋅
   ⋅    ⋅    ⋅     ⋅         ⋅     0.1   ⋅         …
   ⋅    ⋅    ⋅     ⋅         ⋅      ⋅   0.0833333
   ⋅    ⋅    ⋅     ⋅         ⋅      ⋅    ⋅
   ⋅    ⋅    ⋅     ⋅         ⋅      ⋅    ⋅
   ⋅    ⋅    ⋅     ⋅         ⋅      ⋅    ⋅
   ⋅    ⋅    ⋅     ⋅         ⋅      ⋅    ⋅         …
   ⋅    ⋅    ⋅     ⋅         ⋅      ⋅    ⋅
   ⋅    ⋅    ⋅     ⋅         ⋅      ⋅    ⋅
  ⋮                                ⋮               

[jishnub]

Should a change like this go to LinearAlgebra? That way, we won't need to special case this for other infinite array types.

[DanielVandH]

I wasn't too sure, the method in LinearAlgebra does just loop over each element and check iszero for each D.diag[i] which is probably more effective than my inv.(D.diag) for finite arrays. What if I instead did

function inv(D::Diagonal{T, <:InfRanges}) where {T} 
    d = D.diag 
    idx = findfirst(==(zero(T)), d) # iszero is broken for infranges
    isnothing(idx) || throw(SingularException(idx))
    return Diagonal(inv.(d))
end

I do agree that having to special case is a bit annoying but I don't think such a change would get through to LinearAlgebra

[jishnub]

In any case, we can add it here, as we would still need it on older versions of Julia

[codecov]

Codecov Report

All modified and coverable lines are covered by tests ✅

Project coverage is 84.97%. Comparing base (d6eaad4) to head (603baf5).
Report is 3 commits behind head on master.

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #177      +/-   ##
==========================================
- Coverage   86.00%   84.97%   -1.03%     
==========================================
  Files           6        6              
  Lines         743      752       +9     
==========================================
  Hits          639      639              
- Misses        104      113       +9     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[dlfivefifty]

Test failure seems spurious