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

Replace all full(X) calls with convert(Array, X) and migrate tests #18850

Closed
wants to merge 6 commits into from

Conversation

Sacha0
Copy link
Member

@Sacha0 Sacha0 commented Oct 8, 2016

This PR reconstitutes and rebases #17079. #17079 was reverted in #17564 due to JuliaStats/DataArrays.jl#208; the series of PRs associated with #17660 should prevent this PR from causing JuliaStats/DataArrays.jl#208 to reoccur (see also JuliaLang/LinearAlgebra.jl#229). Best!

@nalimilan
Copy link
Member

Any reason not to use the compact Array(A) form instead?

@Sacha0
Copy link
Member Author

Sacha0 commented Oct 9, 2016

Any reason not to use the compact Array(A) form instead?

Good catch! No reason. I wrote convert(Array, by force of habit. Array( seems much better. Will fix. Thanks!

@Sacha0 Sacha0 force-pushed the repfull branch 3 times, most recently from 29dcbd8 to 51f2e80 Compare October 9, 2016 23:08
@Sacha0
Copy link
Member Author

Sacha0 commented Oct 10, 2016

Replaced convert(Array, with Array( everywhere relevant. Also fixed a rebase mistake, so this PR now passes CI. (The lone AV x86_64 failure was a git-related error in the pkg tests.) Best!

@kshyatt kshyatt added linear algebra Linear algebra test This change adds or pertains to unit tests labels Oct 10, 2016
@@ -38,7 +38,7 @@ Bidiagonal(dv::AbstractVector, ev::AbstractVector) = throw(ArgumentError("did yo
Constructs an upper (`uplo='U'`) or lower (`uplo='L'`) bidiagonal matrix using the
given diagonal (`dv`) and off-diagonal (`ev`) vectors. The result is of type `Bidiagonal`
and provides efficient specialized linear solvers, but may be converted into a regular
matrix with [`full`](:func:`full`). `ev`'s length must be one less than the length of `dv`.
matrix with [`Array(_)`](:func:`convert`). `ev`'s length must be one less than the length of `dv`.
Copy link
Contributor

Choose a reason for hiding this comment

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

which should these be linking to?

Copy link
Member Author

@Sacha0 Sacha0 Oct 15, 2016

Choose a reason for hiding this comment

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

convert? Would ...may be converted into a regular matrix with [convert(Array, )](:func:convert) (or Array() for short)... be better? Thoughts? Thanks!

A_mul_B!{T<:BlasFloat}(A::LQ{T}, B::QR{T}) = A[:L]*LAPACK.ormlq!('L','N',A.factors,A.τ,full(B))
A_mul_B!{T<:BlasFloat}(A::QR{T}, B::LQ{T}) = A_mul_B!(zeros(full(A)), full(A), full(B))
A_mul_B!{T<:BlasFloat}(A::LQ{T}, B::QR{T}) = A[:L]*LAPACK.ormlq!('L','N',A.factors,A.τ,Array(B))
A_mul_B!{T<:BlasFloat}(A::QR{T}, B::LQ{T}) = A_mul_B!(zeros(Array(A)), Array(A), Array(B))
Copy link
Contributor

Choose a reason for hiding this comment

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

wouldn't zeros(size(A)) be equivalent and allocate less?

Copy link
Member Author

Choose a reason for hiding this comment

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

Good call. Fixed on push. Thanks!

(./)(A::Number, B::SparseMatrixCSC) = (./)(A, Array(B))
(./)(A::SparseMatrixCSC, B::Array) = (./)(Array(A), B)
(./)(A::Array, B::SparseMatrixCSC) = (./)(A, Array(B))
(./)(A::SparseMatrixCSC, B::SparseMatrixCSC) = (./)(Array( A), Array(B))
Copy link
Contributor

Choose a reason for hiding this comment

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

excess space in Array( A)

Copy link
Member Author

Choose a reason for hiding this comment

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

Good catch. Fixed on push. Thanks!

@@ -119,7 +119,7 @@ for elty in (Float32, Float64, Complex64, Complex128, Int)
invFsv = Fs\vv
x = Ts\vv
@test x ≈ invFsv
@test full(full(Tldlt)) ≈ Fs
@test Array(Array(Tldlt)) ≈ Fs
Copy link
Contributor

Choose a reason for hiding this comment

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

what's going on here?

Copy link
Contributor

Choose a reason for hiding this comment

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

Are you asking why we have two full calls? It's because the first full goes LDLt -> Tridiagonal, and the second goes Tridiagonal -> Matrix

Copy link
Contributor

Choose a reason for hiding this comment

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

ah right my "use a different name for..." issue. what type does Array(::LDLt) return?

Copy link
Member Author

Choose a reason for hiding this comment

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

Derp is what's going on here. Fixed on push. (Array(::LDLt) returns an Array.) Thanks!

@kshyatt
Copy link
Contributor

kshyatt commented Oct 12, 2016

If it's consistent with the old API it should be one of Tridiagonal or
SymTridiagonal

On Wed, Oct 12, 2016, 11:41 AM Tony Kelman notifications@github.com wrote:

@tkelman commented on this pull request.

In test/linalg/tridiag.jl #18850:

@@ -119,7 +119,7 @@ for elty in (Float32, Float64, Complex64, Complex128, Int)
invFsv = Fs\vv
x = Ts\vv
@test x ≈ invFsv

  •        @test full(full(Tldlt)) ≈ Fs
    
  •        @test Array(Array(Tldlt)) ≈ Fs
    

ah right my "use a different name for..." issue. what type does
Array(::LDLt) return?


You are receiving this because you commented.
Reply to this email directly, view it on GitHub
#18850, or mute the thread
https://github.com/notifications/unsubscribe-auth/AAyk43tKILsHBdhRWBxSP5EPGdJuabNuks5qzSnegaJpZM4KR2SL
.

@tkelman
Copy link
Contributor

tkelman commented Oct 12, 2016

the return type of full(::LDLt{eltype, arraytype}) on 0.4 was the arraytype. Need to check on 0.5. Converting to Array should always return an Array (and full on 0.6-dev as of a few weeks ago does this), so this operation that returns the arraytype from a factorization doesn't have its own name right now.

@kshyatt
Copy link
Contributor

kshyatt commented Oct 12, 2016

@tkelman probably better to open a new issue/find that old one for this, but... unfactorize?

@andreasnoack
Copy link
Member

...or full 😄 . In one of the issues, it was concluded that using convert(Type{AbstractMatrix}, x::Factorization) would work for factorizations.

@tkelman
Copy link
Contributor

tkelman commented Oct 13, 2016

Are the factorizations themselves AbstractArrays?

There's a bug here, the PR that reimplemented full in terms of convert did not do so correctly for the factorization types, and now this test change is not testing the original intent either.

@Sacha0
Copy link
Member Author

Sacha0 commented Oct 15, 2016

In JuliaLang/LinearAlgebra.jl#229 consensus was to replace full(factorization)'s behavior with convert(AbstractMatrix/AbstractArray, factorization). #17066 realized that change,

julia> A = SymTridiagonal(2*ones(4), ones(3));

julia> FA = factorize(A)
Base.LinAlg.LDLt{Float64,SymTridiagonal{Float64}}([2.0 0.5 0.0 0.0; 0.5 1.5 0.666667 0.0; 0.0 0.666667 1.33333 0.75; 0.0 0.0 0.75 1.25])

julia> convert(AbstractMatrix, FA)
4×4 SymTridiagonal{Float64}:
 2.0  1.0       
 1.0  2.0  1.0   
     1.0  2.0  1.0
         1.0  2.0

julia> convert(Array, FA)
4×4 Array{Float64,2}:
 2.0  1.0  0.0  0.0
 1.0  2.0  1.0  0.0
 0.0  1.0  2.0  1.0
 0.0  0.0  1.0  2.0

but @tkelman is correct that something is amiss with full(factorization)'s reimplementation as a child of convert(AbstractMatrix, factorization):

julia> full(FA)
4×4 Array{Float64,2}:
 2.0  1.0  0.0  0.0
 1.0  2.0  1.0  0.0
 0.0  1.0  2.0  1.0
 0.0  0.0  1.0  2.0

Having a look at that now (edit: found bug, PR forthcoming), but should be orthogonal to this PR. Thanks and best!

@Sacha0
Copy link
Member Author

Sacha0 commented Oct 15, 2016

#18938 indicates that some of the tests touched in this PR should be more strict. Similarly, some of the full(X) calls presently replaced by Array(X) might rather mean AbstractArray(X) where X is a factorization.

Additionally, we never resolved the full(Triangular/Symmetric/Hermitian) question discussed in #17066, so some full(X) calls where X is one of those wrapper types might not mean Array(X).

Perhaps breaking this PR up into a series of smaller PRs would make this transition safer / smoother? (To decouple the simple, safe cases from the tricky ones?) Thoughts? Thanks and best!

@andreasnoack
Copy link
Member

Additionally, we never resolved the full(Triangular/Symmetric/Hermitian) question discussed in #17066, so some full(X) calls where X is one of those wrapper types might not mean Array(X).

We are looking for a generic function that returns the wrapped array type with elements equal to the wrapper. That function used to be called full so the conclusion to use explicit type conversion instead of a generic function might have been too hasty. The name full might have been wrong but the functionality seems right.

An alternative solution would be to have a mattype function similar to the eltype. Then we could convert with mattype(A)(A) which should work for Symmetric etc. as well as Factorizations.

@Sacha0
Copy link
Member Author

Sacha0 commented Oct 15, 2016

That function used to be called full so the conclusion to use explicit type conversion instead of a generic function might have been too hasty. The name full might have been wrong but the functionality seems right.

So far we've disambiguated two of full's meanings, convert(Array, X) and convert(AbstractArray, X). That we have not yet disambiguated full's third meaning does not imply that disambiguating the two preceding meanings was unwise. To the contrary, only through this disambiguation process have we plumbed the depths of full's depravity :). And the number and variety of ambiguities / challenges we've encountered only emphasizes full's shortcomings / the value of displacing it with better-defined entities.

An alternative solution would be to have a mattype function similar to the eltype. Then we could convert with mattype(A)(A) which should work for Symmetric etc. as well as Factorizations.

Yes, in #17066 we agreed that this approach to handling full(Triangular/Symmetric/...) might work well (but discussion faltered). Annotated matrices and factorizations are fundamentally different objects though, and full's semantics for factorizations (almost) uniformly conform to convert(AbstractArray, X) or convert(Array, X).

Seems we have consensus that something like convert(mattype(A), A) for full(Triangular/Symmetric/...) is the path forward (independent of the name mattype)?

Best!

@andreasnoack
Copy link
Member

I don't question that a renaming was a good idea but I think there is a common pattern among the annotation types and the Factorizations that full was able to capture with some success. That I decided to use the same name as we're already using for densifying sparse matrices was probably not a good idea but I think dedicated generic function instead of convert was useful.

Thanks for reminding of #17066 (comment) in which I conclude that parenttype/mattype is problematic. A solution could be to define the mattype function recursively such that mattype(A::Symmetric)=mattype(A.data) and mattype{T}(::Matrix{T})=Matrix{T}, mattype(A:SubArray)=mattype(parent(A)). But again, this would lead to the behavior that full used to have for the annotation types and Factorizations so maybe renaming full for those types is the simpler solution.

@Sacha0
Copy link
Member Author

Sacha0 commented Oct 26, 2016

The series of pull requests listed just above should cover all files in which all instances of full translate cleanly to Array/AbstractArray. In the remaining files, replacing some instances of full requires consensus on the point discussed above (what full for Symmetric/Hermitian/AbstractTriangular should translate to). Thoughts on how to resolve that issue so that we can move this effort forward? Thanks!

@andreasnoack
Copy link
Member

Could a way forward be to try to define a storagetype/mattype/parenttype function for the annotation types? Then we can punt a bit on what this function should do for Factorizations. I think this xtype function should call itself for arrays with an underlying array (XTriangular, Symmetric, and SubArray) and return the array type for e.g. Array, SparseMatrixCSC, DistributedArray, and SArray. Would you be able to try out such a PR?

@Sacha0
Copy link
Member Author

Sacha0 commented Oct 28, 2016

Sounds good, and shall do! Thanks!

@Sacha0
Copy link
Member Author

Sacha0 commented Jun 29, 2017

With JuliaLang/LinearAlgebra.jl#229 on the 1.0 milestone, perhaps this PR should receive a 1.0 tag as well?

@tkelman tkelman added this to the 1.0 milestone Jun 29, 2017
@JeffBezanson
Copy link
Member

In some cases would it make sense to convert to DenseArray? E.g. if you had a distributed sparse array and wanted to get a distributed dense array.

@Sacha0
Copy link
Member Author

Sacha0 commented Jul 6, 2017

In some cases would it make sense to convert to DenseArray? E.g. if you had a distributed sparse array and wanted to get a distributed dense array.

Absolutely! :) Replacing full uniformly with Array as in this (outdated) pull request is incorrect. full should instead transform to Array, AbstractArray, DenseArray, or one of a number of other possibilities depending on context. (The later, smaller pull requests that reference this pull request (see above) transform full more appropriately / generally, though I imagine some of the transformations to Array in those pull requests could instead have been DenseArray or some other broader type.) Best!

@Sacha0
Copy link
Member Author

Sacha0 commented Oct 24, 2017

Closed by #24250 and the series of preceding pull requests. Best!

@Sacha0 Sacha0 closed this Oct 24, 2017
@Sacha0 Sacha0 deleted the repfull branch October 24, 2017 16:10
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
arrays [a, r, r, a, y, s] design Design of APIs or of the language itself linear algebra Linear algebra needs decision A decision on this change is needed test This change adds or pertains to unit tests
Projects
None yet
Development

Successfully merging this pull request may close these issues.

7 participants