diff --git a/src/operators.jl b/src/operators.jl index 937548c7414..e98fad9b9dc 100644 --- a/src/operators.jl +++ b/src/operators.jl @@ -360,10 +360,21 @@ function _multiplyt!(ret::Array{T}, lhs::Matrix, rhs::SparseMatrixCSC) where T<: ret end - +# See https://github.com/JuliaLang/julia/issues/27015 +function Base.Matrix(S::SparseMatrixCSC{VariableRef}) + A = zeros(AffExpr, S.m, S.n) + for Sj in 1:S.n + for Sk in nzrange(S, Sj) + Si = S.rowval[Sk] + Sv = S.nzval[Sk] + A[Si, Sj] = Sv + end + end + return A +end # TODO: implement sparse * sparse code as in base/sparse/linalg.jl (spmatmul) -_multiply!(ret::AbstractArray{T}, lhs::SparseMatrixCSC, rhs::SparseMatrixCSC) where {T<:JuMPTypes} = _multiply!(ret, lhs, full(rhs)) -_multiplyt!(ret::AbstractArray{T}, lhs::SparseMatrixCSC, rhs::SparseMatrixCSC) where {T<:JuMPTypes} = _multiplyt!(ret, lhs, full(rhs)) +_multiply!(ret::AbstractArray{T}, lhs::SparseMatrixCSC, rhs::SparseMatrixCSC) where {T<:JuMPTypes} = _multiply!(ret, lhs, Matrix(rhs)) +_multiplyt!(ret::AbstractArray{T}, lhs::SparseMatrixCSC, rhs::SparseMatrixCSC) where {T<:JuMPTypes} = _multiplyt!(ret, lhs, Matrix(rhs)) _multiply!(ret, lhs, rhs) = A_mul_B!(ret, lhs, rhs) _multiplyt!(ret, lhs, rhs) = At_mul_B!(ret, lhs, rhs) diff --git a/src/quadexpr.jl b/src/quadexpr.jl index 8f70e8b86a1..70b80025074 100644 --- a/src/quadexpr.jl +++ b/src/quadexpr.jl @@ -69,6 +69,11 @@ function isequal_canonical(quad::GenericQuadExpr{CoefType,VarType}, other::Gener vset = Set((v1,v2)) d[vset] = c + get(d, vset, zero(CoefType)) end + for k in keys(d) + if iszero(d[k]) + delete!(d, k) + end + end return d end d1 = canonicalize(quad) diff --git a/test/old/operator.jl b/test/old/operator.jl index afde85d38c5..4a67a11c340 100644 --- a/test/old/operator.jl +++ b/test/old/operator.jl @@ -72,276 +72,6 @@ Base.transpose(t::MySumType) = MySumType(t.a) return true end - @testset "Vectorized operations" begin - - @testset "Transpose" begin - m = Model() - @variable(m, x[1:3]) - @variable(m, y[1:2,1:3]) - @variable(m, z[2:5]) - @test vec_eq(x', [x[1] x[2] x[3]]) - @test vec_eq(transpose(x), [x[1] x[2] x[3]]) - @test vec_eq(y', [y[1,1] y[2,1] - y[1,2] y[2,2] - y[1,3] y[2,3]]) - @test vec_eq(transpose(y), - [y[1,1] y[2,1] - y[1,2] y[2,2] - y[1,3] y[2,3]]) - @test_throws ErrorException z' - @test_throws ErrorException transpose(z) - end - - @testset "Vectorized arithmetic" begin - m = Model() - @variable(m, x[1:3]) - A = [2 1 0 - 1 2 1 - 0 1 2] - B = sparse(A) - @variable(m, X11) - @variable(m, X23) - X = sparse([1, 2], [1, 3], [X11, X23], 3, 3) # for testing Variable - @variable(m, Xd[1:3, 1:3]) - Y = sparse([1, 2], [1, 3], [2X11, 4X23], 3, 3) # for testing GenericAffExpr - Yd = [2X11 0 0 - 0 0 4X23 - 0 0 0] - Z = sparse([1, 2], [1, 3], [X11^2, 2X23^2], 3, 3) # for testing GenericQuadExpr - Zd = [X11^2 0 0 - 0 0 2X23^2 - 0 0 0] - v = [4, 5, 6] - @test vec_eq(A*x, [2x[1] + x[2] - 2x[2] + x[1] + x[3] - x[2] + 2x[3]]) - @test vec_eq(A*x, B*x) - @test vec_eq(A*x, @JuMP.Expression(B*x)) - @test vec_eq(@JuMP.Expression(A*x), @JuMP.Expression(B*x)) - @test vec_eq(x'*A, [2x[1]+x[2]; 2x[2]+x[1]+x[3]; x[2]+2x[3]]') - @test vec_eq(x'*A, x'*B) - @test vec_eq(x'*A, @JuMP.Expression(x'*B)) - @test vec_eq(@JuMP.Expression(x'*A), @JuMP.Expression(x'*B)) - @test vec_eq(x'*A*x, [2x[1]*x[1] + 2x[1]*x[2] + 2x[2]*x[2] + 2x[2]*x[3] + 2x[3]*x[3]]) - @test vec_eq(x'A*x, x'*B*x) - @test vec_eq(x'*A*x, @JuMP.Expression(x'*B*x)) - @test vec_eq(@JuMP.Expression(x'*A*x), @JuMP.Expression(x'*B*x)) - - y = A*x - @test vec_eq(-x, [-x[1], -x[2], -x[3]]) - @test vec_eq(-y, [-2x[1] - x[2] - -x[1] - 2x[2] - x[3] - -x[2] - 2x[3]]) - @test vec_eq(y + 1, [2x[1] + x[2] + 1 - x[1] + 2x[2] + x[3] + 1 - x[2] + 2x[3] + 1]) - @test vec_eq(y - 1, [2x[1] + x[2] - 1 - x[1] + 2x[2] + x[3] - 1 - x[2] + 2x[3] - 1]) - @test vec_eq(y + 2ones(3), [2x[1] + x[2] + 2 - x[1] + 2x[2] + x[3] + 2 - x[2] + 2x[3] + 2]) - @test vec_eq(y - 2ones(3), [2x[1] + x[2] - 2 - x[1] + 2x[2] + x[3] - 2 - x[2] + 2x[3] - 2]) - @test vec_eq(2ones(3) + y, [2x[1] + x[2] + 2 - x[1] + 2x[2] + x[3] + 2 - x[2] + 2x[3] + 2]) - @test vec_eq(2ones(3) - y, [-2x[1] - x[2] + 2 - -x[1] - 2x[2] - x[3] + 2 - -x[2] - 2x[3] + 2]) - @test vec_eq(y + x, [3x[1] + x[2] - x[1] + 3x[2] + x[3] - x[2] + 3x[3]]) - @test vec_eq(x + y, [3x[1] + x[2] - x[1] + 3x[2] + x[3] - x[2] + 3x[3]]) - @test vec_eq(2y + 2x, [6x[1] + 2x[2] - 2x[1] + 6x[2] + 2x[3] - 2x[2] + 6x[3]]) - @test vec_eq(y - x, [ x[1] + x[2] - x[1] + x[2] + x[3] - x[2] + x[3]]) - @test vec_eq(x - y, [-x[1] - x[2] - -x[1] - x[2] - x[3] - -x[2] - x[3]]) - @test vec_eq(y + x[:], [3x[1] + x[2] - x[1] + 3x[2] + x[3] - x[2] + 3x[3]]) - @test vec_eq(x[:] + y, [3x[1] + x[2] - x[1] + 3x[2] + x[3] - x[2] + 3x[3]]) - - @test vec_eq(@JuMP.Expression(A*x/2), A*x/2) - @test vec_eq(X*v, [4X11; 6X23; 0]) - @test vec_eq(v'*X, [4X11 0 5X23]) - @test vec_eq(v.'*X, [4X11 0 5X23]) - @test vec_eq(X'*v, [4X11; 0; 5X23]) - @test vec_eq(X.'*v, [4X11; 0; 5X23]) - @test vec_eq(X*A, [2X11 X11 0 - 0 X23 2X23 - 0 0 0 ]) - @test vec_eq(A*X, [2X11 0 X23 - X11 0 2X23 - 0 0 X23]) - @test vec_eq(A*X', [2X11 0 0 - X11 X23 0 - 0 2X23 0]) - @test vec_eq(X'*A, [2X11 X11 0 - 0 0 0 - X23 2X23 X23]) - @test vec_eq(X.'*A, [2X11 X11 0 - 0 0 0 - X23 2X23 X23]) - @test vec_eq(A'*X, [2X11 0 X23 - X11 0 2X23 - 0 0 X23]) - @test vec_eq(X.'*A, X'*A) - @test vec_eq(A.'*X, A'*X) - @test vec_eq(X*A, X*B) - @test vec_eq(Y'*A, Y.'*A) - @test vec_eq(A*Y', A*Y.') - @test vec_eq(Z'*A, Z.'*A) - @test vec_eq(Xd'*Y, Xd.'*Y) - @test vec_eq(Y'*Xd, Y.'*Xd) - @test vec_eq(Xd'*Xd, Xd.'*Xd) - # @test_broken vec_eq(A*X, B*X) - # @test_broken vec_eq(A*X', B*X') - @test vec_eq(X'*A, X'*B) - # @test_broken(X'*X, X.'*X) # sparse quadratic known to be broken, see #912 - end - - @testset "Dot-ops" begin - m = Model() - @variable(m, x[1:2,1:2]) - A = [1 2; - 3 4] - B = sparse(A) - y = SparseMatrixCSC(2, 2, copy(B.colptr), copy(B.rowval), vec(x)) - @test vec_eq(A.+x, [1+x[1,1] 2+x[1,2]; - 3+x[2,1] 4+x[2,2]]) - @test vec_eq(A.+x, B.+x) - # @test vec_eq(A.+x, A.+y) == true - # @test vec_eq(A.+y, B.+y) == true - @test vec_eq(x.+A, [1+x[1,1] 2+x[1,2]; - 3+x[2,1] 4+x[2,2]]) - @test vec_eq(x.+A, x.+B) == true - @test vec_eq(x.+A, y.+A) - @test vec_eq(x .+ x, [2x[1,1] 2x[1,2]; 2x[2,1] 2x[2,2]]) - # @test vec_eq(y.+A, y.+B) == true - @test vec_eq(A.-x, [1-x[1,1] 2-x[1,2]; - 3-x[2,1] 4-x[2,2]]) - @test vec_eq(A.-x, B.-x) - @test vec_eq(A.-x, A.-y) - @test vec_eq(x .- x, [zero(AffExpr) for _1 in 1:2, _2 in 1:2]) - # @test vec_eq(A.-y, B.-y) == true - @test vec_eq(x.-A, [-1+x[1,1] -2+x[1,2]; - -3+x[2,1] -4+x[2,2]]) - @test vec_eq(x.-A, x.-B) - @test vec_eq(x.-A, y.-A) - # @test vec_eq(y.-A, y.-B) == true - @test vec_eq(A.*x, [1*x[1,1] 2*x[1,2]; - 3*x[2,1] 4*x[2,2]]) - @test vec_eq(A.*x, B.*x) - @test vec_eq(A.*x, A.*y) - # @test vec_eq(A.*y, B.*y) == true - - @test vec_eq(x.*A, [1*x[1,1] 2*x[1,2]; - 3*x[2,1] 4*x[2,2]]) - @test vec_eq(x.*A, x.*B) - @test vec_eq(x.*A, y.*A) - # @test vec_eq(y.*A, y.*B) == true - - @test vec_eq(x .* x, [x[1,1]^2 x[1,2]^2; x[2,1]^2 x[2,2]^2]) - @test_throws ErrorException vec_eq(A./x, [1*x[1,1] 2*x[1,2]; - 3*x[2,1] 4*x[2,2]]) - @test vec_eq(x./A, [1/1*x[1,1] 1/2*x[1,2]; - 1/3*x[2,1] 1/4*x[2,2]]) - @test vec_eq(x./A, x./B) - @test vec_eq(x./A, y./A) - # @test vec_eq(A./y, B./y) == true - - @test vec_eq((2*x) / 3, full((2*y) / 3)) - @test vec_eq(2 * (x/3), full(2 * (y/3))) - @test vec_eq(x[1,1] * A, full(x[1,1] * B)) - end - - @testset "Vectorized comparisons" begin - m = Model() - @variable(m, x[1:3]) - A = [1 2 3 - 0 4 5 - 6 0 7] - B = sparse(A) - if VERSION < v"0.6.0-dev.2074" # julia PR #19670 - @constraint(m, x'*A*x .>= 1) - else - # force vector output - @constraint(m, reshape(x,(1,3))*A*x .>= 1) - end - @test vec_eq(m.quadconstr[1].terms, [x[1]*x[1] + 2x[1]*x[2] + 4x[2]*x[2] + 9x[1]*x[3] + 5x[2]*x[3] + 7x[3]*x[3] - 1]) - @test m.quadconstr[1].sense == :(>=) - if VERSION < v"0.6.0-dev.2074" # julia PR #19670 - @constraint(m, x'*A*x .>= 1) - else - @constraint(m, x'*A*x >= 1) - end - @test vec_eq(m.quadconstr[1].terms, m.quadconstr[2].terms) - - mat = [ 3x[1] + 12x[3] + 4x[2] - 2x[1] + 12x[2] + 10x[3] - 15x[1] + 5x[2] + 21x[3]] - - @constraint(m, (x'A)' + 2A*x .<= 1) - terms = map(v->v.terms, m.linconstr[1:3]) - lbs = map(v->v.lb, m.linconstr[1:3]) - ubs = map(v->v.ub, m.linconstr[1:3]) - @test vec_eq(terms, mat) - @test lbs == fill(-Inf, 3) - @test ubs == fill( 1, 3) - @test vec_eq((x'A)' + 2A*x, (x'A)' + 2B*x) - @test vec_eq((x'A)' + 2A*x, (x'B)' + 2A*x) - @test vec_eq((x'A)' + 2A*x, (x'B)' + 2B*x) - @test vec_eq((x'A)' + 2A*x, @JuMP.Expression((x'A)' + 2A*x)) - @test vec_eq((x'A)' + 2A*x, @JuMP.Expression((x'B)' + 2A*x)) - @test vec_eq((x'A)' + 2A*x, @JuMP.Expression((x'A)' + 2B*x)) - @test vec_eq((x'A)' + 2A*x, @JuMP.Expression((x'B)' + 2B*x)) - - @constraint(m, -1 .<= (x'A)' + 2A*x .<= 1) - terms = map(v->v.terms, m.linconstr[4:6]) - lbs = map(v->v.lb, m.linconstr[4:6]) - ubs = map(v->v.ub, m.linconstr[4:6]) - @test vec_eq(terms, mat) == true - @test lbs == fill(-1, 3) - @test ubs == fill( 1, 3) - - @constraint(m, -[1:3;] .<= (x'A)' + 2A*x .<= 1) - terms = map(v->v.terms, m.linconstr[7:9]) - lbs = map(v->v.lb, m.linconstr[7:9]) - ubs = map(v->v.ub, m.linconstr[7:9]) - @test vec_eq(terms, mat) == true - @test lbs == -[1:3;] - @test ubs == fill( 1, 3) - - @constraint(m, -[1:3;] .<= (x'A)' + 2A*x .<= [3:-1:1;]) - terms = map(v->v.terms, m.linconstr[10:12]) - lbs = map(v->v.lb, m.linconstr[10:12]) - ubs = map(v->v.ub, m.linconstr[10:12]) - @test vec_eq(terms, mat) == true - @test lbs == -[1:3;] - @test ubs == [3:-1:1;] - - @constraint(m, -[1:3;] .<= (x'A)' + 2A*x .<= 3) - terms = map(v->v.terms, m.linconstr[13:15]) - lbs = map(v->v.lb, m.linconstr[13:15]) - ubs = map(v->v.ub, m.linconstr[13:15]) - @test vec_eq(terms, mat) == true - @test lbs == -[1:3;] - @test ubs == fill(3,3) - end - - end - @testset "JuMPArray concatenation" begin m = Model() @variable(m, x[1:3]) diff --git a/test/operator.jl b/test/operator.jl index c50cb157e5a..6b47d2024b8 100644 --- a/test/operator.jl +++ b/test/operator.jl @@ -247,4 +247,283 @@ using Compat.Test end end end + + vec_eq(x,y) = vec_eq([x;], [y;]) + + function vec_eq(x::AbstractArray, y::AbstractArray) + size(x) == size(y) || return false + for i in 1:length(x) + v, w = convert(AffExpr,x[i]), convert(AffExpr,y[i]) + JuMP.isequal_canonical(v, w) || return false + end + return true + end + + function vec_eq(x::Array{QuadExpr}, y::Array{QuadExpr}) + size(x) == size(y) || return false + for i in 1:length(x) + JuMP.isequal_canonical(x[i], y[i]) || return false + end + return true + end + + @testset "Vectorized operations" begin + @testset "Transpose" begin + m = Model() + @variable(m, x[1:3]) + @variable(m, y[1:2,1:3]) + @variable(m, z[2:5]) + @test vec_eq(x', [x[1] x[2] x[3]]) + @test vec_eq(transpose(x), [x[1] x[2] x[3]]) + @test vec_eq(y', [y[1,1] y[2,1] + y[1,2] y[2,2] + y[1,3] y[2,3]]) + @test vec_eq(transpose(y), + [y[1,1] y[2,1] + y[1,2] y[2,2] + y[1,3] y[2,3]]) + @test (z')' == z + @test transpose(transpose(z)) == z + end + + @testset "Vectorized arithmetic" begin + m = Model() + @variable(m, x[1:3]) + A = [2 1 0 + 1 2 1 + 0 1 2] + B = sparse(A) + @variable(m, X11) + @variable(m, X23) + X = sparse([1, 2], [1, 3], [X11, X23], 3, 3) # for testing Variable + @test vec_eq([X11 0. 0.; 0. 0. X23; 0. 0. 0.], @inferred Matrix(X)) + @variable(m, Xd[1:3, 1:3]) + Y = sparse([1, 2], [1, 3], [2X11, 4X23], 3, 3) # for testing GenericAffExpr + Yd = [2X11 0 0 + 0 0 4X23 + 0 0 0] + Z = sparse([1, 2], [1, 3], [X11^2, 2X23^2], 3, 3) # for testing GenericQuadExpr + Zd = [X11^2 0 0 + 0 0 2X23^2 + 0 0 0] + v = [4, 5, 6] + @test vec_eq(A*x, [2x[1] + x[2] + 2x[2] + x[1] + x[3] + x[2] + 2x[3]]) + @test vec_eq(A*x, B*x) + @test vec_eq(A*x, @JuMP.Expression(B*x)) + @test vec_eq(@JuMP.Expression(A*x), @JuMP.Expression(B*x)) + @test vec_eq(x'*A, [2x[1]+x[2]; 2x[2]+x[1]+x[3]; x[2]+2x[3]]') + @test vec_eq(x'*A, x'*B) + @test vec_eq(x'*A, @JuMP.Expression(x'*B)) + @test vec_eq(@JuMP.Expression(x'*A), @JuMP.Expression(x'*B)) + @test JuMP.isequal_canonical(x'*A*x, 2x[1]*x[1] + 2x[1]*x[2] + 2x[2]*x[2] + 2x[2]*x[3] + 2x[3]*x[3]) + @test JuMP.isequal_canonical(x'A*x, x'*B*x) + @test JuMP.isequal_canonical(x'*A*x, @JuMP.Expression(x'*B*x)) + @test JuMP.isequal_canonical(@JuMP.Expression(x'*A*x), @JuMP.Expression(x'*B*x)) + + y = A*x + @test vec_eq(-x, [-x[1], -x[2], -x[3]]) + @test vec_eq(-y, [-2x[1] - x[2] + -x[1] - 2x[2] - x[3] + -x[2] - 2x[3]]) + @test vec_eq(y + 1, [2x[1] + x[2] + 1 + x[1] + 2x[2] + x[3] + 1 + x[2] + 2x[3] + 1]) + @test vec_eq(y - 1, [2x[1] + x[2] - 1 + x[1] + 2x[2] + x[3] - 1 + x[2] + 2x[3] - 1]) + @test vec_eq(y + 2ones(3), [2x[1] + x[2] + 2 + x[1] + 2x[2] + x[3] + 2 + x[2] + 2x[3] + 2]) + @test vec_eq(y - 2ones(3), [2x[1] + x[2] - 2 + x[1] + 2x[2] + x[3] - 2 + x[2] + 2x[3] - 2]) + @test vec_eq(2ones(3) + y, [2x[1] + x[2] + 2 + x[1] + 2x[2] + x[3] + 2 + x[2] + 2x[3] + 2]) + @test vec_eq(2ones(3) - y, [-2x[1] - x[2] + 2 + -x[1] - 2x[2] - x[3] + 2 + -x[2] - 2x[3] + 2]) + @test vec_eq(y + x, [3x[1] + x[2] + x[1] + 3x[2] + x[3] + x[2] + 3x[3]]) + @test vec_eq(x + y, [3x[1] + x[2] + x[1] + 3x[2] + x[3] + x[2] + 3x[3]]) + @test vec_eq(2y + 2x, [6x[1] + 2x[2] + 2x[1] + 6x[2] + 2x[3] + 2x[2] + 6x[3]]) + @test vec_eq(y - x, [ x[1] + x[2] + x[1] + x[2] + x[3] + x[2] + x[3]]) + @test vec_eq(x - y, [-x[1] - x[2] + -x[1] - x[2] - x[3] + -x[2] - x[3]]) + @test vec_eq(y + x[:], [3x[1] + x[2] + x[1] + 3x[2] + x[3] + x[2] + 3x[3]]) + @test vec_eq(x[:] + y, [3x[1] + x[2] + x[1] + 3x[2] + x[3] + x[2] + 3x[3]]) + + @test vec_eq(@JuMP.Expression(A*x/2), A*x/2) + @test vec_eq(X*v, [4X11; 6X23; 0]) + @test vec_eq(v'*X, [4X11 0 5X23]) + @test vec_eq(v.'*X, [4X11 0 5X23]) + @test vec_eq(X'*v, [4X11; 0; 5X23]) + @test vec_eq(X.'*v, [4X11; 0; 5X23]) + @test vec_eq(X*A, [2X11 X11 0 + 0 X23 2X23 + 0 0 0 ]) + @test vec_eq(A*X, [2X11 0 X23 + X11 0 2X23 + 0 0 X23]) + @test vec_eq(A*X', [2X11 0 0 + X11 X23 0 + 0 2X23 0]) + @test vec_eq(X'*A, [2X11 X11 0 + 0 0 0 + X23 2X23 X23]) + @test vec_eq(X.'*A, [2X11 X11 0 + 0 0 0 + X23 2X23 X23]) + @test vec_eq(A'*X, [2X11 0 X23 + X11 0 2X23 + 0 0 X23]) + @test vec_eq(X.'*A, X'*A) + @test vec_eq(A.'*X, A'*X) + @test vec_eq(X*A, X*B) + @test vec_eq(Y'*A, Y.'*A) + @test vec_eq(A*Y', A*Y.') + @test vec_eq(Z'*A, Z.'*A) + @test vec_eq(Xd'*Y, Xd.'*Y) + @test vec_eq(Y'*Xd, Y.'*Xd) + @test vec_eq(Xd'*Xd, Xd.'*Xd) + @test vec_eq(A*X, B*X) + @test_broken vec_eq(A*X', B*X') # See https://github.com/JuliaOpt/JuMP.jl/issues/1276 + @test vec_eq(X'*A, X'*B) + @test vec_eq(X'*X, X.'*X) + end + + @testset "Dot-ops" begin + m = Model() + @variable(m, x[1:2,1:2]) + A = [1 2; + 3 4] + B = sparse(A) + y = SparseMatrixCSC(2, 2, copy(B.colptr), copy(B.rowval), vec(x)) + @test vec_eq(A.+x, [1+x[1,1] 2+x[1,2]; + 3+x[2,1] 4+x[2,2]]) + @test vec_eq(A.+x, B.+x) + @test vec_eq(A.+x, A.+y) + @test vec_eq(A.+y, B.+y) + @test vec_eq(x.+A, [1+x[1,1] 2+x[1,2]; + 3+x[2,1] 4+x[2,2]]) + @test vec_eq(x.+A, x.+B) + @test vec_eq(x.+A, y.+A) + @test vec_eq(x .+ x, [2x[1,1] 2x[1,2]; 2x[2,1] 2x[2,2]]) + @test vec_eq(y.+A, y.+B) + @test vec_eq(A.-x, [1-x[1,1] 2-x[1,2]; + 3-x[2,1] 4-x[2,2]]) + @test vec_eq(A.-x, B.-x) + @test vec_eq(A.-x, A.-y) + @test vec_eq(x .- x, [zero(AffExpr) for _1 in 1:2, _2 in 1:2]) + @test vec_eq(A.-y, B.-y) + @test vec_eq(x.-A, [-1+x[1,1] -2+x[1,2]; + -3+x[2,1] -4+x[2,2]]) + @test vec_eq(x.-A, x.-B) + @test vec_eq(x.-A, y.-A) + @test vec_eq(y.-A, y.-B) + @test vec_eq(A.*x, [1*x[1,1] 2*x[1,2]; + 3*x[2,1] 4*x[2,2]]) + @test vec_eq(A.*x, B.*x) + @test vec_eq(A.*x, A.*y) + @test vec_eq(A.*y, B.*y) + + @test vec_eq(x.*A, [1*x[1,1] 2*x[1,2]; + 3*x[2,1] 4*x[2,2]]) + @test vec_eq(x.*A, x.*B) + @test vec_eq(x.*A, y.*A) + @test vec_eq(y.*A, y.*B) + + @test vec_eq(x .* x, [x[1,1]^2 x[1,2]^2; x[2,1]^2 x[2,2]^2]) + @test_throws ErrorException vec_eq(A./x, [1*x[1,1] 2*x[1,2]; + 3*x[2,1] 4*x[2,2]]) + @test vec_eq(x./A, [1/1*x[1,1] 1/2*x[1,2]; + 1/3*x[2,1] 1/4*x[2,2]]) + @test vec_eq(x./A, x./B) + @test vec_eq(x./A, y./A) + @test_throws ErrorException A./y + @test_throws ErrorException B./y + + @test vec_eq((2*x) / 3, full((2*y) / 3)) + @test vec_eq(2 * (x/3), full(2 * (y/3))) + @test vec_eq(x[1,1] .* A, full(x[1,1] .* B)) + end + + @testset "Vectorized comparisons" begin + m = Model() + @variable(m, x[1:3]) + A = [1 2 3 + 0 4 5 + 6 0 7] + B = sparse(A) + # force vector output + cref1 = @constraint(m, reshape(x,(1,3))*A*x .>= 1) + c1 = JuMP.constraintobject.(cref1, QuadExpr, MOI.GreaterThan) + f1 = map(c -> c.func, c1) + @test vec_eq(f1, [x[1]*x[1] + 2x[1]*x[2] + 4x[2]*x[2] + 9x[1]*x[3] + 5x[2]*x[3] + 7x[3]*x[3]]) + @test all(c -> c.set.lower == 1, c1) + + cref2 = @constraint(m, x'*A*x >= 1) + c2 = JuMP.constraintobject.(cref2, QuadExpr, MOI.GreaterThan) + @test vec_eq(f1[1], c2.func) + + mat = [ 3x[1] + 12x[3] + 4x[2] + 2x[1] + 12x[2] + 10x[3] + 15x[1] + 5x[2] + 21x[3]] + + cref3 = @constraint(m, (x'A)' + 2A*x .<= 1) + c3 = JuMP.constraintobject.(cref3, AffExpr, MOI.LessThan) + f3 = map(c->c.func, c3) + @test vec_eq(f3, mat) + @test all(c -> c.set.upper == 1, c3) + @test vec_eq((x'A)' + 2A*x, (x'A)' + 2B*x) + @test vec_eq((x'A)' + 2A*x, (x'B)' + 2A*x) + @test vec_eq((x'A)' + 2A*x, (x'B)' + 2B*x) + @test vec_eq((x'A)' + 2A*x, @JuMP.Expression((x'A)' + 2A*x)) + @test vec_eq((x'A)' + 2A*x, @JuMP.Expression((x'B)' + 2A*x)) + @test vec_eq((x'A)' + 2A*x, @JuMP.Expression((x'A)' + 2B*x)) + @test vec_eq((x'A)' + 2A*x, @JuMP.Expression((x'B)' + 2B*x)) + + cref4 = @constraint(m, -1 .<= (x'A)' + 2A*x .<= 1) + c4 = JuMP.constraintobject.(cref4, AffExpr, MOI.Interval) + f4 = map(c->c.func, c4) + @test vec_eq(f4, mat) + @test all(c -> c.set.lower == -1, c4) + @test all(c -> c.set.upper == 1, c4) + + cref5 = @constraint(m, -[1:3;] .<= (x'A)' + 2A*x .<= 1) + c5 = JuMP.constraintobject.(cref5, AffExpr, MOI.Interval) + f5 = map(c->c.func, c5) + @test vec_eq(f5, mat) + @test map(c -> c.set.lower, c5) == -[1:3;] + @test all(c -> c.set.upper == 1, c4) + + cref6 = @constraint(m, -[1:3;] .<= (x'A)' + 2A*x .<= [3:-1:1;]) + c6 = JuMP.constraintobject.(cref6, AffExpr, MOI.Interval) + f6 = map(c->c.func, c6) + @test vec_eq(f6, mat) + @test map(c -> c.set.lower, c6) == -[1:3;] + @test map(c -> c.set.upper, c6) == [3:-1:1;] + + cref7 = @constraint(m, -[1:3;] .<= (x'A)' + 2A*x .<= 3) + c7 = JuMP.constraintobject.(cref7, AffExpr, MOI.Interval) + f7 = map(c->c.func, c7) + @test vec_eq(f7, mat) + @test map(c -> c.set.lower, c7) == -[1:3;] + @test all(c -> c.set.upper == 3, c7) + end + end end