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

sparse vectors join the higher order function party #19690

Merged
merged 3 commits into from
Dec 31, 2016

Conversation

Sacha0
Copy link
Member

@Sacha0 Sacha0 commented Dec 23, 2016

This pull request extends sparse map[!]/broadcast[!] to sparse vectors and sparse vector/matrix combinations.

Specifically, this pull request's first commit: (1) introduces a common interface to SparseVector and SparseMatrixCSC sufficient for the purposes of map[!]/broadcast[!]; (2) rewrites the existing sparse map[!]/broadcast[!] code against that interface; and (3) places the relevant code in the new SparseArrays.HigherOrderFns submodule which lives in base/sparse/higherorderfns.jl, loaded after definition of both SparseVector and SparseMatrixCSC.

This pull request's second commit: (1) condenses and systematizes existing tests for generic sparse map[!]/broadcast[!]; (2) extends those tests to sparse vectors and sparse vector/matrix combinations; and (3) collects those tests, as well older tests of sparse broadcast[!], into the new test file test/sparse/higherorderfns.jl corresponding to base/sparse/higherorderfns.jl.

Performance is more or less the same (previous generic sparse map[!]/broadcast[!] methods versus those in this pull request), faster in some cases where I fixed a type stability regression.

Best!

(Tangentially, this code (and likely an appreciable fraction of sparsevector.jl/sparsematrix.jl) could be further unified and simplified by revising SparseVector/SparseMatrixCSC with a common interface (cleaner than that in this pull request). If that idea sparks interest in this thread, I'll write a mini-julep at some point.)

@Sacha0 Sacha0 added sparse Sparse arrays broadcast Applying a function over a collection labels Dec 23, 2016
@Sacha0 Sacha0 added this to the 0.6.0 milestone Dec 23, 2016
@stevengj
Copy link
Member

@nanosoldier runbenchmarks(ALL, vs = ":master")

@nanosoldier
Copy link
Collaborator

Your benchmark job has completed - possible performance regressions were detected. A full report can be found here. cc @jrevels

@test all(A[1,:] .* 3 .== AF[1,:] .* 3)
#@test A[:,1] .* 3 == AF[:,1] .* 3
@test all(A[:,1] .* 3 .== AF[:,1] .* 3)
#TODO: simple comparation with == returns false because the left side is a (two-dimensional) SparseMatrixCSC
Copy link
Contributor

Choose a reason for hiding this comment

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

this todo is really old, predates SparseVector being in base

Copy link
Member Author

Choose a reason for hiding this comment

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

(Is this a request to clear this TODO in this pull request? Or a request to clear it later? Thanks!)

Copy link
Contributor

Choose a reason for hiding this comment

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

I was about to do this in a local branch, but it would cause a conflict 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.

Happy to sort that however you think best --- let me know. Thanks!

Copy link
Contributor

@tkelman tkelman Dec 23, 2016

Choose a reason for hiding this comment

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

it's just uncommenting, and deleting the all and todo lines

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. Thanks!

@stevengj
Copy link
Member

Benchmarks look good; the floatexp benchmarks seem noisy as usual.

@stevengj
Copy link
Member

Why do we still have broadcast(::typeof(+), ...) methods for SparseVector?

@stevengj
Copy link
Member

Do we still need the methods specialized on a pair of sparse matrices?

@stevengj
Copy link
Member

It still doesn't handle mixtures of sparse vectors/matrices and scalars, correct?

@Sacha0
Copy link
Member Author

Sacha0 commented Dec 23, 2016

Why do we still have broadcast(::typeof(+), ...) methods for SparseVector?

base/sparse/sparsevectors.jl contains a few not-exported map-like methods and several children [1] (e.g. the methods you mention). This pull request should enable removal of some of that code without regressions and more with potentially negligible regressions. That cleanup seems a substantial project in itself (with potentially controversial points, and largely orthogonal to this work in practice). So in the interest of keeping this pull request confined / uncontroversial / readily mergeable, I left that work separate.

[1] ... mostly defined against AbstractSparseVector, though most of those methods are SparseVector-specific?

Do we still need the methods specialized on a pair of sparse matrices?

Yes: The generic methods still appreciably lag the specialized methods in some cases. (A few performance improvements I have in mind might close the gap enough to nix some of the specialized methods. Post feature freeze work though.)

It still doesn't handle mixtures of sparse vectors/matrices and scalars, correct?

Correct: #19641 / merging #19667 blocks my solution for handling combinations including broadcast scalars. Best!

@test broadcast(+, X, Y) == sparse(broadcast(+, fX, fY))
@test broadcast(*, X, Y) == sparse(broadcast(*, fX, fY))
@test broadcast(f, X, Y) == sparse(broadcast(f, fX, fY))
try Base.Broadcast.broadcast_indices(spzeros((shapeX .- 1)...), Y)
Copy link
Contributor

Choose a reason for hiding this comment

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

split this over multiple lines. if the intention is that the try should always throw, then add a @test false after so it doesn't silently pass if it happens not to throw

Copy link
Member Author

@Sacha0 Sacha0 Dec 24, 2016

Choose a reason for hiding this comment

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

The intention is for the @test_throws in the catch body to only evaluate when the statement in the try body throws. Is there a better way to express this? Thanks!

(In other words, the test logic is: Check whether the arguments are shape-broadcast-incompatible. If so, make certain the method being @test_throws'd throws.)

Copy link
Contributor

Choose a reason for hiding this comment

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

Is there another way of determining shape compatibility than checking for something throwing? Calling some internal boolean function maybe? Otherwise seems like you might not be checking carefully enough which inputs do or do not trigger the test

Copy link
Member Author

@Sacha0 Sacha0 Dec 24, 2016

Choose a reason for hiding this comment

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

Is there another way of determining shape compatibility than checking for something throwing? Calling some internal boolean function maybe?

To the best of my knowledge none such exists. (I should instead call check_broadcast_indices for clarity, but that similarly throws or returns nothing. On second thought no, broadcast_indices is correct as check_broadcast_indices requires the broadcast shape (which you retrieve from broadcast_indices) as a argument.) Thoughts? Thanks!

Copy link
Member Author

Choose a reason for hiding this comment

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

(Split over multiple lines. Thanks!)

Copy link
Contributor

Choose a reason for hiding this comment

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

hm maybe we should make that check function return a boolean instead of nothing-or-throw?

Copy link
Member Author

Choose a reason for hiding this comment

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

Seems @timholy designed that mechanism --- perhaps ping him? (Orthogonal to this pull request in any case?) Thanks!

@@ -1674,168 +1596,9 @@ end
# 19304
@inferred hcat(sparse(rand(2,1)), eye(2,2))

# Test that broadcast[!](f, [C::SparseMatrixCSC], A::SparseMatrixCSC, B::SparseMatrixCSC)
# returns the correct (densely populated) result when f(zero(eltype(A)), zero(eltype(B))) != 0
Copy link
Contributor

Choose a reason for hiding this comment

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

which set of the new tests replaces this?

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 broadcast[!] implementation specialized for pairs of (input) sparse vectors/matrices testset replaces those tests.

X = sparse(broadcast(f, Array(A), Array(B), Array(C)))
broadcast!(f, copy(X), A, B, C) # warmup for @allocated
@test_broken (@allocated broadcast!(f, X, A, B, C)) == 0
# this last test allocates 16 bytes in the entry point for broadcast!, but none of the
Copy link
Contributor

@tkelman tkelman Dec 24, 2016

Choose a reason for hiding this comment

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

figure out what was happening here? nevermind

@Sacha0
Copy link
Member Author

Sacha0 commented Dec 26, 2016

Absent objections or requests for time, I plan to merge this Wednesday morning PDT. Best!

@stevengj
Copy link
Member

stevengj commented Dec 26, 2016

You can now remove the workaround in this test, yes, and just do @test exact_equal(x .* x, abs2(x)), or abs2.(x)?

@Sacha0
Copy link
Member Author

Sacha0 commented Dec 26, 2016

You can now remove the workaround in this test, yes, and just do @test exact_equal(x .* x, abs2(x)), or abs2.(x)?

Happily yes, and done. (The abs2.(x) change is in #18564.) Thanks!

(If you or @tkelman would be comfortable merging this PR prior to Wednesday, please feel welcome.)

@tkelman
Copy link
Contributor

tkelman commented Dec 26, 2016

I have not finished reviewing this yet, but would like to do so today or tomorrow.

edit: got through it all, locally undoing some of the rearrangements as https://github.com/JuliaLang/julia/compare/d8a57182f49787055b3a864cc86401988be4bcdf...tkelman:tk/sparsevechof?w=1 just for the purposes of reviewing made it a little easier to compare

@test broadcast(*, X, Y, Z) == sparse(broadcast(*, fX, fY, fZ))
@test broadcast(f, X, Y, Z) == sparse(broadcast(f, fX, fY, fZ))
try Base.Broadcast.broadcast_indices(spzeros((shapeX .- 1)...), Y, Z)
catch @test_throws DimensionMismatch broadcast(+, spzeros((shapeX .- 1)...), Y, Z) end
Copy link
Contributor

Choose a reason for hiding this comment

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

this try-catch (and the one a few lines lower) would also read more clearly if on multiple lines

Copy link
Contributor

Choose a reason for hiding this comment

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

though I still think it would be preferable to have this serve as an independent verification of how broadcast_indices is expected to behave for the differently-sized inputs, rather than relying on it throwing - if it were to be refactored to no longer throw, then this test would silently stop doing anything

Copy link
Member Author

Choose a reason for hiding this comment

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

this try-catch (and the one a few lines lower) would also read more clearly if on multiple lines

Thanks for catching these other two instances! Split into multiple lines on push.

though I still think it would be preferable to have this serve as an independent verification of how broadcast_indices is expected to behave for the differently-sized inputs, rather than relying on it throwing - if it were to be refactored to no longer throw, then this test would silently stop doing anything

Agreed. Not being certain how best to achieve that strengthening at this time, can we revisit that point later? Thanks!

Copy link
Contributor

Choose a reason for hiding this comment

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

leave a todo comment

# entry point for broadcast! with --track-allocation=all, but that first line
# almost certainly should not allocate. so not certain what's going on.
@test broadcast!(f, Q, X, Y, Z) == sparse(broadcast!(f, fQ, fX, fY, fZ))
# --> test shape checks for both braodcast and broadcast! entry points
Copy link
Contributor

Choose a reason for hiding this comment

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

"braodcast" typo

Copy link
Member Author

Choose a reason for hiding this comment

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

Fixed on push. Thanks!

# separate horizontal and vertical expansion
@test broadcast(op, A, B, C) == sparse(broadcast(op, fA, fB, fC))
@test broadcast!(op, X, A, B, C) == sparse(broadcast!(op, fX, fA, fB, fC))
# simultaneous horizontal and vertical expansion
Copy link
Contributor

Choose a reason for hiding this comment

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

does this require 4 inputs or is it covered by the new loop with only 3?

Copy link
Member Author

Choose a reason for hiding this comment

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

Three input arguments suffice to test the relevant code paths, and indeed the testset titled "broadcast[!] implementation capable of handling >2 (input) sparse vectors/matrices" covers the relevant code. Thanks!

# (5) Define _map_[not]zeropres! specialized for a pair of (input) sparse vectors/matrices.
# (6) Define general _map_[not]zeropres! capable of handling >2 (input) sparse vectors/matrices.
# (7) Define _broadcast_[not]zeropres! specialized for a pair of (input) sparse vectors/matrices.
# (8) Define general _broadcast_[not]zeropres! capabel of handling >2 (input) sparse vectors/matrices.
Copy link
Contributor

Choose a reason for hiding this comment

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

"capabel" typo

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!

end
# helper functions for these methods and some of those below
@inline _densecoloffsets(A::SparseVector) = 0
@inline _densecoloffsets(A::SparseMatrixCSC) = 0:A.m:(A.m*A.n - 1)
Copy link
Contributor

Choose a reason for hiding this comment

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

why isn't the end point (A.m * (A.n-1)) 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.

Facepalm-worthy unintentional code obfuscation is why. Your suggestion is much better. Fixed on push. Thanks!

Copy link
Contributor

Choose a reason for hiding this comment

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

A.m * (0:A.n-1) might be even better, not sure

@inbounds for j in 1:C.n
C.colptr[j] = Ck
ks = _startindforbccol_all(j, expandshorzs, As)
stopks = _stopindforbccol_all(j + 1, expandshorzs, As)
Copy link
Contributor

Choose a reason for hiding this comment

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

was the j+1 here a mistake or did you change the definition of this helper function?

Copy link
Member Author

Choose a reason for hiding this comment

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

Changed the helper function's definition IIRC, but my memory is vague on that point. Thanks!

…e vector/matrix combinations.

Extend generic sparse map[!]/broadcast[!] to sparse vectors and sparse vector/matrix combinations. Do so by introducing a common interface to SparseVector and SparseMatrixCSC for the purposes of map[!]/broadcast[!], and rewriting sparse map[!]/broadcast[!] against that interface. Relocate that code to a separate file/module base/sparse/higherorderfns.jl/SparseArrays.HigherOrderFns, loaded after definition of both SparseVector and SparseMatrixCSC.
Condense and systematize existing tests for generic sparse map[!]/broadcast[!], and extend to sparse vectors and vector/matrix combinations. Relocate new test code to a separate file test/sparse/higherorderfns.jl corresponding to base/sparse/higherorderfns.jl. (Test/sparse/sparsevector.jl is hypothetically confined to SparseVectors, and test/sparse/sparse.jl mostly dedicated to SparseMatrixCSCs.) Move older tests of sparse broadcast[!] into that new file as well.
@Sacha0
Copy link
Member Author

Sacha0 commented Dec 30, 2016

Absent objections I plan to merge this once CI approves. Best!

@Sacha0 Sacha0 merged commit 1ed4545 into JuliaLang:master Dec 31, 2016
@Sacha0 Sacha0 deleted the sparsevechof branch December 31, 2016 03:43
@Sacha0
Copy link
Member Author

Sacha0 commented Dec 31, 2016

(AV i686 failure unrelated.)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
broadcast Applying a function over a collection sparse Sparse arrays
Projects
None yet
Development

Successfully merging this pull request may close these issues.

4 participants