-
Notifications
You must be signed in to change notification settings - Fork 151
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
Porting SDiagonal from Bridge.jl #240
Conversation
Originally contributed by D. Getz (https://github.com/getzdan), M. Schauer at https://github.com/mschauer/Bridge.jl under MIT License
I add tests the following days if nobody beats me. |
Codecov Report
@@ Coverage Diff @@
## master #240 +/- ##
==========================================
- Coverage 88.11% 87.81% -0.31%
==========================================
Files 33 34 +1
Lines 2247 2306 +59
==========================================
+ Hits 1980 2025 +45
- Misses 267 281 +14
Continue to review full report at Codecov.
|
I put some more work into this, this is probably halfway done, help is most welcome. |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Thanks for all of this - appreciated very much.
Don't be scared by the number of comments... they are each pretty simple. :)
src/SDiagonal.jl
Outdated
|
||
Base.@propagate_inbounds function getindex{N,T}(D::SDiagonal{N,T}, i::Int, j::Int) | ||
@boundscheck checkbounds(D, i, j) | ||
if i == j |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This is a perfect place to use ifelse
instead of branching (will produce faster code)
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
also looks nicer
src/SDiagonal.jl
Outdated
end | ||
end | ||
|
||
# avoid linear indexing? |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Would be nice. Unfortunately most the static arrays internal code assumes a linear indexing style - we should change this, but it won't be easy...
src/SDiagonal.jl
Outdated
convert{N,T}(::Type{SDiagonal{N,T}}, D::SDiagonal{N,T}) = D | ||
convert{N,T}(::Type{SDiagonal{N,T}}, D::SDiagonal) = SDiagonal{N,T}(convert(SVector{N,T}, D.diag)) | ||
|
||
size{N}(D::SDiagonal{N}) = (N,N) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Is this necessary? There should be a fallback already.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
ok
src/SDiagonal.jl
Outdated
|
||
|
||
convert{N,T}(::Type{SDiagonal{N,T}}, D::SDiagonal{N,T}) = D | ||
convert{N,T}(::Type{SDiagonal{N,T}}, D::SDiagonal) = SDiagonal{N,T}(convert(SVector{N,T}, D.diag)) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
D::SDiagonal{N}
maybe?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
ok
src/SDiagonal.jl
Outdated
# Originally contributed by D. Getz (https://github.com/getzdan), M. Schauer | ||
# at https://github.com/mschauer/Bridge.jl under MIT License | ||
|
||
import Base: getindex,setindex!,==,-,+,*,/,\,transpose,ctranspose,convert, size, abs, real, imag, conj, eye, inv |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
merge these imports with the global ones?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Still this
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
ok
src/SDiagonal.jl
Outdated
/{T<:Number}(D::SDiagonal, x::T) = SDiagonal(D.diag / x) | ||
*(Da::SDiagonal, Db::SDiagonal) = SDiagonal(Da.diag .* Db.diag) | ||
*(D::SDiagonal, V::SVector) = D.diag .* V | ||
*(V::SVector, D::SDiagonal) = D.diag .* V |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Should V
be a RowVector
?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
yes
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
or rather: this is already a fallback https://github.com/JuliaLang/julia/blob/master/base/linalg/rowvector.jl#L181 which looks like it is no overhead
src/SDiagonal.jl
Outdated
*(Da::SDiagonal, Db::SDiagonal) = SDiagonal(Da.diag .* Db.diag) | ||
*(D::SDiagonal, V::SVector) = D.diag .* V | ||
*(V::SVector, D::SDiagonal) = D.diag .* V | ||
*(A::SMatrix, D::SDiagonal) = scalem(A,D.diag) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
A::AbstractMatrix
(and below)
src/SDiagonal.jl
Outdated
*(V::SVector, D::SDiagonal) = D.diag .* V | ||
*(A::SMatrix, D::SDiagonal) = scalem(A,D.diag) | ||
*(D::SDiagonal, A::SMatrix) = scalem(D.diag,A) | ||
\(D::SDiagonal, b::SVector) = D.diag .\ b |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Can we make this work for matrices as well?
also /(vector, diagonal)
?
src/SDiagonal.jl
Outdated
logm(D::SDiagonal) = SDiagonal(log.(D.diag)) | ||
sqrtm(D::SDiagonal) = SDiagonal(sqrt.(D.diag)) | ||
|
||
\(D::SDiagonal, B::SMatrix) = scalem(1 ./ D.diag, B) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
B::AbstractMatrix
(and below)
src/SDiagonal.jl
Outdated
/(Da::SDiagonal, Db::SDiagonal) = SDiagonal(Da.diag ./ Db.diag ) | ||
|
||
function inv{N,T}(D::SDiagonal{N,T}) | ||
for i = 1:N |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Possibly @inbounds
?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Also, we generally unroll such loops elsewhere...
It's also slightly annoying that the code to check if we should throw a singular exception will be slower than taking the inverse itself... but I guess Base.Diagonal
does this too?
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Yes, but I guess the divisions will still be the most expensive operation here.
Not at all, thank you for the careful reading! |
src/SDiagonal.jl
Outdated
|
||
@generated function scalem{M, N}(a::StaticMatrix{M,N}, b::StaticVector{N}) | ||
expr = vec([:(a[$j,$i]*b[$i]) for j=1:M, i=1:N]) | ||
:(let val1 = ($(expr[1])); similar_type(SMatrix{M,N},typeof(val1))(val1, $(expr[2:end]...)); end) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Perhaps these should be inlined?
:(@_inline_meta; let ... )
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
ok
src/SDiagonal.jl
Outdated
convert{N,T}(::Type{SDiagonal{N,T}}, D::SDiagonal{N,T}) = D | ||
convert{N,T}(::Type{SDiagonal{N,T}}, D::SDiagonal{N}) = SDiagonal{N,T}(convert(SVector{N,T}, D.diag)) | ||
|
||
Base.@propagate_inbounds function getindex{N,T}(D::SDiagonal{N,T}, i::Int, j::Int) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
@propagate_inbounds
is needed if you want to defer bounds checking to an inner function. This should be simply @inline
since you are doing your own @boundscheck
and @inbounds
here.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
ok
A final style note is that we've been moving to using |
src/SDiagonal.jl
Outdated
end | ||
|
||
# avoid linear indexing? | ||
Base.@propagate_inbounds function getindex{N,T}(D::SDiagonal{N,T}, k::Int) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This is a good use of @propagate_inbounds
.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
However, it has already been imported so you don't need the Base.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
ok
@andyferris I hope I addressed all your comments. |
It's looking good. There are some more imports we could move to the main package file, that's about it that I see. |
I left them there because they are imported locally elsewhere: https://github.com/JuliaArrays/StaticArrays.jl/search?utf8=✓&q=import&type= We could centralize these imports independently. |
Fine - no worries then :) |
Thanks again |
@@ -0,0 +1,106 @@ | |||
# Originally contributed by D. Getz (https://github.com/getzdan), M. Schauer |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
By the way, if you're porting someone else's code, and you want to attribute them, it's also quite possible to make them the author in git using something like git commit --author="Some Body <[email protected]>"
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
When I do this kind of thing, I'd be inclined to add the original chunk of code (in probably non-working form) under the other author's name as a single commit (@mention
ing them to make sure they're happy with that). Then add any necessary changes under your own name in further commits.
The nice thing about doing it this way is you avoid baking authorship into comments which people feel they can't remove (like the one above), but you also get to do proper attribution which I feel is quite important.
Just some thoughts, I'm happy this was merged already.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
When you port entire files (like in this case) it is also possible to preserve the full git history http://gbayer.com/development/moving-files-from-one-git-repository-to-another-preserving-history/
Ref JuliaLang/julia#22718 for a potential future rewrite to utilize |
|
||
# Bad input | ||
@test_throws Exception SMatrix{1,Int}() | ||
@test_throws Exception SMatrix{2,Int}((1,)) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
@test_throws Exception
is generally a very bad idea. It will pass if the code has a typo.
* Porting SDiagonal from Bridge.jl Originally contributed by D. Getz (https://github.com/getzdan), M. Schauer at https://github.com/mschauer/Bridge.jl under MIT License * Inherit from AbstractMatrix * Fixes and tests for constructors * More tests, more fixes * Address andyferris' comments * SDiagonal: next round of comments * Remaining tests for SDiagonal * Where-notation in SDiagonal * Test for inv(::SDiagonal)
This was originally written for FixedSizeArrays and ported to StaticArrays by @getzdan. As I think we both still aren't very familiar with the internals of StaticArrays please review carefully.
Tag #235