Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Deprecate getindex(::Factorization, ::Symbol) in favor of dot overloading #25184

Merged
merged 3 commits into from
Dec 22, 2017

Conversation

andreasnoack
Copy link
Member

@andreasnoack andreasnoack commented Dec 19, 2017

Now that we have dot overloading, it is possible to use that instead of the special getindex methods for Factorizations. Furthermore, with a little care in the getproperty methods, it is possible to make the constant propagation work such that qr and lu are now inferred even though they extract the factors from the Factorization.

julia> A = randn(3,3);

julia> F = lufact(A);

julia> F[:L]
3×3 Array{Float64,2}:
  1.0        0.0       0.0
 -0.368965   1.0       0.0
  0.414762  -0.531993  1.0

becomes

julia> F.L
3×3 Array{Float64,2}:
  1.0        0.0       0.0
 -0.368965   1.0       0.0
  0.414762  -0.531993  1.0

Fixes #25159

@andreasnoack andreasnoack added the domain:linear algebra Linear algebra label Dec 19, 2017
@andreasnoack andreasnoack force-pushed the anj/getprops branch 2 times, most recently from badc334 to 7da8002 Compare December 19, 2017 13:03
Copy link
Sponsor Member

@StefanKarpinski StefanKarpinski left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Gorgeous! Why are all the getproperty methods marked as @inline?

3×3 Array{Float64,2}:
0.0 0.0 0.0
0.0 0.0 0.0
0.0 0.0 0.0
```
"""
function getindex(B::BunchKaufman{T}, d::Symbol) where {T<:BlasFloat}
@inline function getproperty(B::BunchKaufman{T}, d::Symbol) where {T<:BlasFloat}
Copy link
Sponsor Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this seems like a lot to inline

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think inlining is necessary to get the constant propagation working sufficiently well to make the factors inferred.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hm. Maybe that wasn't the issue. It seems to work now. I'll try to remove the annotations and see what happens.

3×3 Array{Float64,2}:
0.0 0.0 0.0
0.0 0.0 0.0
0.0 0.0 0.0

julia> F = bkfact(Symmetric(A));

julia> F[:U]*F[:D]*F[:U]' - F[:P]*A*F[:P]'
julia> F.U*F.D*F.U' - F.P*A*F.P'
Copy link
Sponsor Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

oh boy, that is soooo much prettier

d == :L && return LowerTriangular(Symbol(C.uplo) == d ? C.factors : C.factors')
d == :p && return C.piv
if d == :P
@inline function getproperty(C::Cholesky, d::Symbol)
Copy link
Sponsor Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

again, a lot to inline

Copy link
Member

@Sacha0 Sacha0 left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Beautiful! :) I will try to work in a proper review this evening.

@c42f
Copy link
Member

c42f commented Dec 19, 2017

This is really nice. Previously I wrapped SVD for StaticArrays and it wasn't really possible to have performance while sharing the Base interface.

@@ -3423,6 +3423,12 @@ workspace() = error("workspace() is discontinued, check out Revise.jl for an alt
# PR #25113
@deprecate_binding CartesianRange CartesianIndices

# Use getproperty instead of getindex for Factorizations
function getindex(F::Factorization, s::Symbol)
depwarn("F[:$s] is deprecated, use F.$s instead.", :getindex)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Given that this method's F may have any name in the caller's code, F[:...] might be a bit obscure. Perhaps expand this depwarn a bit? A sketch:

Extracting factorization components via getindex(F::Factorization, s::Symbol) methods, usually written F[:s] where F is the <:Factorization and :s identifies the factorization component (for example aQRfact[:Q]), has been deprecated in favor of F.s.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Missed or disagree? :)

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Decided to think a little longer and then forgot. I've decided that I prefer the current shorter version. If users later report that it is unclear we can adjust but I wouldn't expect that.

size(B::BunchKaufman) = size(B.LD)
size(B::BunchKaufman, d::Integer) = size(B.LD, d)
size(B::BunchKaufman) = size(getfield(B, :LD))
size(B::BunchKaufman, d::Integer) = size(getfield(B, :LD), d)
Copy link
Member

@Sacha0 Sacha0 Dec 20, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why switch to explicitly calling getfield here?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

To avoid an infinite cycle between size and getproperty. The reason is that getproperty always calls size.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Cheers, and thanks for the explanation! :)

Cfactors = getfield(C, :factors)
Cuplo = getfield(C, :uplo)
if d == :U
return UpperTriangular(Symbol(Cuplo) == d ? Cfactors : Cfactors')
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Assuming you would like an eager adjoint of Cfactors, Cfactors' -> adjoint(Cfactors)?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I didn't consider the possibility. I'm not sure. Should it really be eager when it could be lazy?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If we want this type stable it needs to be eager, right?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

True. I thought it would be known at compile time but that is of course not the case for Cuplo. I'll change it.

if d == :U
return UpperTriangular(Symbol(Cuplo) == d ? Cfactors : Cfactors')
elseif d == :L
return LowerTriangular(Symbol(Cuplo) == d ? Cfactors : Cfactors')
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Likewise here re. eager versus lazy adjoint :).

Cfactors = getfield(C, :factors)
Cuplo = getfield(C, :uplo)
if d == :U
return UpperTriangular(Symbol(Cuplo) == d ? Cfactors : Cfactors')
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

And likewise here re. eager versus lazy adjoint :).

if d == :U
return UpperTriangular(Symbol(Cuplo) == d ? Cfactors : Cfactors')
elseif d == :L
return LowerTriangular(Symbol(Cuplo) == d ? Cfactors : Cfactors')
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

And likewise here re. eager versus lazy adjoint :).

d == :values && return A.values
d == :vectors && return A.vectors
throw(KeyError(d))
end
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

👍

end
end

getindex(A::LQPackedQ, i::Integer, j::Integer) =
mul!(A, setindex!(zeros(eltype(A), size(A, 2)), 1, j))[i]

getq(A::LQ) = LQPackedQ(A.factors, A.τ)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

💯! Perhaps getq deserves a deprecation?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done

sQf1 = size(Q.factors, 1)
return (!full ? Array(Q) : mul!(Q, Matrix{eltype(Q)}(I, sQf1, sQf1))), R
end
function _qr(A::Union{Number, AbstractMatrix}, ::Val{true}; full::Bool = false)
F = qrfact(A, Val(true))
Q, R, p = getq(F), F[:R]::Matrix{eltype(F)}, F[:p]::Vector{BlasInt}
Q, R, p = F.Q, F.R, F.p
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Fantastic.

return F.Vt'
function getproperty(F::SVD, d::Symbol)
if d == :V
return getfield(F, :Vt)'
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

adjoint(getfield(F, :Vt)) for eager adjoint?

finalizer(free!, A)
return A
end
end
Sparse(p::Ptr{C_Sparse{Tv}}) where {Tv<:VTypes} = Sparse{Tv}(p)

Base.unsafe_convert(::Type{Ptr{Tv}}, A::Sparse{Tv}) where {Tv} = A.p
Base.unsafe_convert(::Type{Ptr{Tv}}, A::Sparse{Tv}) where {Tv} = getfield(A, :ptr)
Copy link
Member

@Sacha0 Sacha0 Dec 20, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Switch necessary? (Edit: I'm guessing so given you perform the same rewrite in the unsafe_convert def below as well?)

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I had to change the name of the field to ptr because F.p is now syntax for extracting a permutation vector.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Nice!

@test_throws DimensionMismatch LinAlg.lowrankupdate(F, ones(eltype($v), length($v)+1))
@test LinAlg.lowrankdowndate(G, $v).$uplo ≈ F.$uplo
@test_throws DimensionMismatch LinAlg.lowrankdowndate(G, ones(eltype($v), length($v)+1))
end
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Perhaps explicit getproperty calls (if possible?) would be simpler than metafication?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The idea was to exercise the compile time version but that is probably not worth the metafication here.

@test convert(Array, usv) ≈ a
@test usv[:Vt]' ≈ usv[:V]
@test_throws KeyError usv[:Z]
@test usv.Vt' ≈ usv.V
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Eager adjoint? :)

Copy link
Member

@Sacha0 Sacha0 left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Wonderful. Simply wonderful :).

@andreasnoack andreasnoack force-pushed the anj/getprops branch 4 times, most recently from 4da58d2 to 894cb01 Compare December 21, 2017 08:14
@Sacha0
Copy link
Member

Sacha0 commented Dec 21, 2017

Rebase and merge? :)

@JeffBezanson
Copy link
Sponsor Member

🎉
Amazing. This is one of those things we've been dreaming of for years :)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

6 participants