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

clean up some linear algebra primitives for JuMP eltypes #1456

Merged
merged 11 commits into from
Sep 15, 2018
Merged

Conversation

jrevels
Copy link
Contributor

@jrevels jrevels commented Sep 5, 2018

Did this while working on some of the TODOs floating about this code (e.g. properly intercepting sparse matmul on Julia >v0.7). I can just continue that work within this PR, or open a different PR on top of this one later - whatever's easiest for review.

Let's see if tests pass with what I currently have...

EDIT: Tests now pass locally

@jrevels jrevels force-pushed the jr/opsupdate branch 3 times, most recently from 65769b4 to afd0e11 Compare September 5, 2018 20:37
Copy link
Member

@mlubin mlubin left a comment

Choose a reason for hiding this comment

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

Please fix the CI failure on 0.6.

src/operators.jl Outdated

densify_with_jump_eltype(x::AbstractMatrix) = convert(Matrix, x)

function densify_with_jump_eltype(x::SparseMatrixCSC{V}) where V<:AbstractVariableRef
Copy link
Member

Choose a reason for hiding this comment

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

Why have special behavior for SparseMatrixCSC{<:AbstractVariableRef} and not AbstractMatrix{<:AbstractVariableRef}?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Just because (AFAICT, maybe I missed something) that was what the code was doing before, when it was pirate-overloading Matrix. I could change to AbstractMatrix{<:AbstractVariableRef} if you want.

src/operators.jl Outdated
return convert(Matrix{GenericAffExpr{Float64,V}}, x)
end

# see https://github.com/JuliaLang/julia/pull/18218
Copy link
Member

Choose a reason for hiding this comment

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

Nit: Proper capitalization and punctuation in comments. "See JuliaLang/julia#18218."

src/operators.jl Outdated
_A_mul_B_ret_dims(A::AbstractMatrix, B::AbstractVector) = (size(A, 1),)
_A_mul_B_ret_dims(A::AbstractMatrix, B::AbstractMatrix) = (size(A, 1), size(B, 2))

# don't do size checks here; defer that to `_A_mul_B(A, B)`
Copy link
Member

Choose a reason for hiding this comment

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

Nit: Proper capitalization and punctuation in comments.

src/operators.jl Outdated
_sizehint_expr!(q, n)
for k ∈ 1:n
tmp = convert(T, lhs[k,i]*rhs[k,j]) # transpose
function fillwithzeros!(A, ::Type{T}) where {T}
Copy link
Member

Choose a reason for hiding this comment

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

Style nit: fill_with_zeros! or fill_with_zeros.

end

###############################################################################
# `_A_mul_B!(ret, A, B)` adds the result of `A*B` into the buffer `ret`. We roll our own
Copy link
Member

Choose a reason for hiding this comment

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

Nit: It's weird for the ############################################################################### separator not to cover the width of the text. Could we remove these? (I know they're used above, but they're pretty old..)

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Sure, just used them since they already seemed to be here.

Is there a preferred "code block comment separator" format for JuMP? I don't necessarily care what format it takes, but breaking the code up with these kind of things vastly improves readability.

Copy link
Member

Choose a reason for hiding this comment

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

"code block comment separator" format for JuMP?

No, not really. I thought about it briefly but don't have a clear answer. Maybe more consistent use of empty lines would make it easier to separate sections without large horizontal comment bars.

src/operators.jl Outdated
function Base.:-(x::AbstractArray{T}) where T<:JuMPTypes
ret = similar(x, typeof(-one(T)))
for I in eachindex(ret)
ret[I] = -x[I]
end
ret
end
Base.:*(x::AbstractArray{T}) where {T<:JuMPTypes} = x

# TODO This will interact poorly with other packages that overload broadcast on Julia >v0.7;
Copy link
Member

Choose a reason for hiding this comment

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

What breaks if we don't define this code on 0.7+? It's not clear why we need to override this.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

good point, I'll try wrapping this in a VERSION conditional

src/operators.jl Outdated

Base.:*(x::AbstractArray{T}) where {T<:JuMPTypes} = x

Base.:/(A::SparseMatrixCSC{T}, B::Number) where {T<:JuMPTypes} = SparseMatrixCSC(A.m, A.n, copy(A.colptr), copy(A.rowval), A.nzval ./ B)
Copy link
Member

Choose a reason for hiding this comment

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

Nit: Use one-line functions only when the function fits on one line (i.e., 80 characters).

src/operators.jl Outdated
end

Base.:*(A::Number, B::SparseMatrixCSC{T}) where {T<:JuMPTypes} = SparseMatrixCSC(B.m, B.n, copy(B.colptr), copy(B.rowval), A .* B.nzval)
Copy link
Member

Choose a reason for hiding this comment

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

Why aren't these covered by the built-ins? Because JuMPTypes aren't Numbers? This merits a comment.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

That's a good question. This code was here before I got here, though :p I'll try to read deeper and figure out an answer

src/operators.jl Show resolved Hide resolved
src/operators.jl Outdated
# ensured that the eltype of `ret` is appropriate, and has zeroed the elements of `ret` (if
# desired).

function _A_mul_B!(ret::AbstractArray{T}, A, B) where {T<:JuMPTypes}
Copy link
Member

Choose a reason for hiding this comment

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

Nit: Spaces between binary operators, so T <: JuMPTypes.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Hmmm. No spaces surrounding <: seems to be the norm in the rest of the codebase. Seems best to change them all at once or not all to avoid inconsistent style - should I do that in this PR or leave until later?

Copy link
Member

@mlubin mlubin Sep 10, 2018

Choose a reason for hiding this comment

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

We're fixing the style issues as we touch lines of code. It's pretty painful to update the whole code base at once, but you're welcome to if you'd like. Otherwise it's ok to have a mix of styles as long as we're moving in the right direction.

Copy link
Contributor Author

@jrevels jrevels Sep 11, 2018

Choose a reason for hiding this comment

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

This seems like an easy one, I'll just go for it

EDIT: I overestimated my regex abilities and underestimated the formatting issues caused by naive find/replace here...so I only made this change to operators.jl.

@jrevels
Copy link
Contributor Author

jrevels commented Sep 6, 2018

On the style nits - sorry, I didn't see http://www.juliaopt.org/JuMP.jl/latest/style.html until just now :) I'll get those fixes in soon.

@codecov
Copy link

codecov bot commented Sep 6, 2018

Codecov Report

Merging #1456 into master will increase coverage by 0.5%.
The diff coverage is 87.7%.

Impacted file tree graph

@@            Coverage Diff            @@
##           master    #1456     +/-   ##
=========================================
+ Coverage   89.22%   89.72%   +0.5%     
=========================================
  Files          26       26             
  Lines        3723     3700     -23     
=========================================
- Hits         3322     3320      -2     
+ Misses        401      380     -21
Impacted Files Coverage Δ
src/operators.jl 93.66% <87.7%> (+5.35%) ⬆️
src/JuMP.jl 78.67% <0%> (+0.31%) ⬆️
src/parseexpr.jl 94.88% <0%> (+0.78%) ⬆️

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 4496a54...33e8359. Read the comment docs.

@jrevels
Copy link
Contributor Author

jrevels commented Sep 10, 2018

Bump. This is ready (pending responses to my previous questions from the initial review).


Base.transpose(x::AbstractJuMPScalar) = x

# Can remove the following code once == overloading is removed

function Compat.LinearAlgebra.issymmetric(x::Matrix{T}) where T<:JuMPTypes
function Compat.LinearAlgebra.issymmetric(x::Matrix{T}) where {T <: JuMPTypes}
Copy link
Member

Choose a reason for hiding this comment

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

Why use { and } when it is not a one line function ?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

stylistic consistency (this is also more readable IMO, but that's subjective)

There probably should just be a rule in the style guide about the desired style here

_sizehint_expr!(q, n)
for k ∈ 1:n
tmp = convert(T, lhs[k,i]*rhs[k,j]) # transpose
function fill_with_zeros!(A, ::Type{T}) where {T}
Copy link
Member

Choose a reason for hiding this comment

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

Isn't A an AbstractArray{T} ? Why is the second argument needed ?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

It's bad practice to dispatch on element type in cases like these; the eltype of A and T only need be related by a subtype relation, not equality.

# Don't do size checks here; defer that to `_A_mul_B(A, B)`.
function _A_mul_B_ret(A, B, dims...)
T = _A_mul_B_eltype(eltype(A), eltype(B))
ret = Array{T}(undef, _A_mul_B_ret_dims(A, B))
Copy link
Member

Choose a reason for hiding this comment

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

Why not do zeros(T, ...) instead of Array... followed by fill_with_zeros! ?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

I just carried this pattern over from the old code. It does leave a nice path to generalize to non-Array rets later, but that's probably premature.

I could replace it with zeros if we want to make that change.

Copy link
Member

Choose a reason for hiding this comment

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

zeros is broken for JuMP types because zero returns a mutable value.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Huh, TIL

julia> mutable struct Foo
       end

julia> Base.zero(::Type{Foo}) = Foo()

julia> x = zeros(Foo, 2)
2-element Array{Foo,1}:
 Foo()
 Foo()

julia> x[1] === x[2]
true

That's arguably a bug in Base's implementation. If I wanted a single instance, I would've just used fill(Foo(), 2)...

end

###############################################################################
# `_A_mul_B!(ret, A, B)` adds the result of `A*B` into the buffer `ret`. We roll our own
Copy link
Member

Choose a reason for hiding this comment

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

"code block comment separator" format for JuMP?

No, not really. I thought about it briefly but don't have a clear answer. Maybe more consistent use of empty lines would make it easier to separate sections without large horizontal comment bars.

test/operator.jl Outdated
@test JuMP.isequal_canonical((2 .* x) ./ 3, Matrix((2 .* y) ./ 3))
@test JuMP.isequal_canonical(2 .* (x ./ 3), Matrix(2 * (y ./ 3)))
@test JuMP.isequal_canonical((x[1,1],) .* A, Matrix((x[1,1],) .* B))
z = JuMP.densify_with_jump_eltype((2 .* y) ./ 3)
Copy link
Member

Choose a reason for hiding this comment

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

Using an internal JuMP function to construct the expected answer seems unnecessary and makes it harder to verify that the test is correct.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Yeah - before this PR, these tests' reliance on internal constructors was obscured by the fact that the constructors were pirated overloads of Matrix. Now that it's easily visible where this is happening, I've added a TODO to refactor the tests

test/operator.jl Outdated
@test JuMP.isequal_canonical((2 .* x) ./ 3, Matrix((2 .* y) ./ 3))
@test JuMP.isequal_canonical(2 .* (x ./ 3), Matrix(2 * (y ./ 3)))
@test JuMP.isequal_canonical((x[1,1],) .* A, Matrix((x[1,1],) .* B))
z = JuMP.densify_with_jump_eltype((2 .* y) ./ 3)
Copy link
Member

Choose a reason for hiding this comment

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

These tests are really poorly structured but I'll let this go for now. Unit tests should test only one thing and not refer to variables like y that are defined far away.

@mlubin mlubin merged commit 7ada455 into master Sep 15, 2018
@jrevels jrevels deleted the jr/opsupdate branch September 15, 2018 15:55
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Development

Successfully merging this pull request may close these issues.

3 participants