diff --git a/docs/src/tutorials/conic/ellipse_approx.jl b/docs/src/tutorials/conic/ellipse_approx.jl index 81fbbf5082d..0d78f1f2722 100644 --- a/docs/src/tutorials/conic/ellipse_approx.jl +++ b/docs/src/tutorials/conic/ellipse_approx.jl @@ -109,7 +109,7 @@ m, n = size(S) [i in 1:m], S[i, :]' * Z * S[i, :] - 2 * S[i, :]' * z + s <= 1, ) -@constraint(model, [t; vec(Z)] in MOI.RootDetConeSquare(n)) +@constraint(model, (t, Z) in MOI.RootDetConeSquare(n)) @objective(model, Max, t) optimize!(model) Test.@test termination_status(model) == OPTIMAL #src @@ -208,9 +208,11 @@ set_silent(model) ## The former constraint S[i, :]' * Z * S[i, :] - 2 * S[i, :]' * z + s <= 1 f = [1 - S[i, :]' * Z * S[i, :] + 2 * S[i, :]' * z - s for i in 1:m] @constraint(model, f in MOI.Nonnegatives(m)) -## The former constraint [t; vec(Z)] in MOI.RootDetConeSquare(n) -Z_upper = [Z[i, j] for j in 1:n for i in 1:j] -@constraint(model, 1 * vcat(t, Z_upper) .+ 0 in MOI.RootDetConeTriangle(n)) +## The former constraint (t, Z) in MOI.RootDetConeSquare(n) +@constraint( + model, + (1.0 * t, LinearAlgebra.Symmetric(1.0 .* Z)) in MOI.RootDetConeTriangle(n), +) ## The former @objective(model, Max, t) @objective(model, Max, 1 * t + 0) optimize!(model) diff --git a/docs/src/tutorials/conic/experiment_design.jl b/docs/src/tutorials/conic/experiment_design.jl index 88edd583dc4..f4d61b05e8b 100644 --- a/docs/src/tutorials/conic/experiment_design.jl +++ b/docs/src/tutorials/conic/experiment_design.jl @@ -210,11 +210,8 @@ set_silent(dOpt) @variable(dOpt, t) @objective(dOpt, Max, t) @constraint(dOpt, sum(np) <= n) -E = V * LinearAlgebra.diagm(0 => np ./ n) * V' -@constraint( - dOpt, - [t, 1, (E[i, j] for i in 1:q for j in 1:i)...] in MOI.LogDetConeTriangle(q) -) +E = LinearAlgebra.Symmetric(V * LinearAlgebra.diagm(0 => np ./ n) * V') +@constraint(dOpt, (t, 1.0, E) in MOI.LogDetConeTriangle(q)) optimize!(dOpt) objective_value(dOpt) #- diff --git a/docs/src/tutorials/conic/min_ellipse.jl b/docs/src/tutorials/conic/min_ellipse.jl index 1fc33027a85..c934cced579 100644 --- a/docs/src/tutorials/conic/min_ellipse.jl +++ b/docs/src/tutorials/conic/min_ellipse.jl @@ -119,7 +119,7 @@ end # the conic reformulation: @variable(model, log_det_P) -@constraint(model, [log_det_P; 1; vec(P²)] in MOI.LogDetConeSquare(n)) +@constraint(model, (log_det_P, 1, P²) in MOI.LogDetConeSquare(n)) @objective(model, Max, log_det_P) # Now, solve the program: diff --git a/docs/src/tutorials/conic/tips_and_tricks.jl b/docs/src/tutorials/conic/tips_and_tricks.jl index 2471a19029b..2d7f7b46bdf 100644 --- a/docs/src/tutorials/conic/tips_and_tricks.jl +++ b/docs/src/tutorials/conic/tips_and_tricks.jl @@ -95,7 +95,7 @@ set_silent(model) @variable(model, x[1:3]) @variable(model, t) @constraint(model, sum(x) == 1) -@constraint(model, [t; x] in SecondOrderCone()) +@constraint(model, (t, x) in SecondOrderCone()) @objective(model, Min, t) optimize!(model) value(t), value.(x) @@ -117,7 +117,7 @@ set_silent(model) @variable(model, θ) @variable(model, t) @expression(model, residuals, θ * data .- target) -@constraint(model, [t; 0.5; residuals] in RotatedSecondOrderCone()) +@constraint(model, (t, 0.5, residuals) in RotatedSecondOrderCone()) @objective(model, Min, t) optimize!(model) value(θ), value(t) @@ -140,7 +140,7 @@ set_silent(model) @variable(model, x == 1.5) @variable(model, z) @objective(model, Min, z) -@constraint(model, [x, 1, z] in MOI.ExponentialCone()) +@constraint(model, (x, 1, z) in MOI.ExponentialCone()) optimize!(model) value(z), exp(1.5) @@ -153,7 +153,7 @@ set_silent(model) @variable(model, x) @variable(model, z == 1.5) @objective(model, Max, x) -@constraint(model, [x, 1, z] in MOI.ExponentialCone()) +@constraint(model, (x, 1, z) in MOI.ExponentialCone()) optimize!(model) value(x), log(1.5) @@ -170,7 +170,7 @@ set_silent(model) @objective(model, Min, t) @variable(model, u[1:N]) @constraint(model, sum(u) <= 1) -@constraint(model, [i = 1:N], [x[i] - t, 1, u[i]] in MOI.ExponentialCone()) +@constraint(model, [i = 1:N], (x[i] - t, 1, u[i]) in MOI.ExponentialCone()) optimize!(model) value(t), log(sum(exp.(x0))) @@ -214,7 +214,7 @@ set_silent(model) @objective(model, Max, sum(t)) @constraint(model, sum(x) == 1) @constraint(model, A * x .<= b) -@constraint(model, [i = 1:n], [t[i], x[i], 1] in MOI.ExponentialCone()) +@constraint(model, [i = 1:n], (t[i], x[i], 1) in MOI.ExponentialCone()) optimize!(model) objective_value(model) @@ -232,7 +232,7 @@ set_silent(model) @objective(model, Max, -t) @constraint(model, sum(x) == 1) @constraint(model, A * x .<= b) -@constraint(model, [t; ones(n); x] in MOI.RelativeEntropyCone(2n + 1)) +@constraint(model, (t, ones(n), x) in MOI.RelativeEntropyCone(2n + 1)) optimize!(model) objective_value(model) @@ -252,7 +252,7 @@ model = Model(SCS.Optimizer) set_silent(model) @variable(model, t) @variable(model, x >= 1.5) -@constraint(model, [t, 1, x] in MOI.PowerCone(1 / 3)) +@constraint(model, (t, 1, x) in MOI.PowerCone(1 / 3)) @objective(model, Min, t) optimize!(model) value(t), value(x) @@ -273,7 +273,7 @@ function p_norm(x::Vector, p) set_silent(model) @variable(model, r[1:N]) @variable(model, t) - @constraint(model, [i = 1:N], [r[i], t, x[i]] in MOI.PowerCone(1 / p)) + @constraint(model, [i = 1:N], (r[i], t, x[i]) in MOI.PowerCone(1 / p)) @constraint(model, sum(r) == t) @objective(model, Min, t) optimize!(model) @@ -334,10 +334,64 @@ set_silent(model) @variable(model, x[1:4]) @variable(model, t) @constraint(model, sum(x) == 1) -@constraint(model, [t; x] in MOI.GeometricMeanCone(5)) +@constraint(model, (t, x) in MOI.GeometricMeanCone(5)) optimize!(model) value(t), value.(x) +# ## RootDetCone + +# The [`MOI.RootDetConeSquare`](@ref) is a cone of the form: +# ```math +# K = \{ (t, X) \in \mathbb{R}^{1+d^2} : t \le \det(X)^{\frac{1}{d}} \} +# ``` + +model = Model(SCS.Optimizer) +set_silent(model) +@variable(model, t) +@variable(model, X[1:2, 1:2]) +@objective(model, Max, t) +@constraint(model, (t, X) in MOI.RootDetConeSquare(2)) +@constraint(model, X .== [2 1; 1 3]) +optimize!(model) +value(t), sqrt(LinearAlgebra.det(value.(X))) + +# If `X` is symmetric, then you can use [`MOI.RootDetConeTriangle`](@ref) +# instead. This can be more efficient because the solver does not need to add +# additional constraints to ensure `X` is symmetric. + +@constraint( + model, + (t, LinearAlgebra.Symmetric(X)) in MOI.RootDetConeTriangle(2), +) + +# ## LogDetCone + +# The [`MOI.LogDetConeSquare`](@ref) is a cone of the form: +# ```math +# K = \{ (t, u, X) \in \mathbb{R}^{2+d^2} : t \le u \log(\det(X / u)) \} +# ``` + +model = Model(SCS.Optimizer) +set_silent(model) +@variable(model, t) +@variable(model, u) +@variable(model, X[1:2, 1:2]) +@objective(model, Max, t) +@constraint(model, (t, u, X) in MOI.LogDetConeSquare(2)) +@constraint(model, X .== [2 1; 1 3]) +@constraint(model, u == 0.5) +optimize!(model) +value(t), 0.5 * log(LinearAlgebra.det(value.(X) ./ 0.5)) + +# If `X` is symmetric, then you can use [`MOI.LogDetConeTriangle`](@ref) +# instead. This can be more efficient because the solver does not need to add +# additional constraints to ensure `X` is symmetric. + +@constraint( + model, + (t, u, LinearAlgebra.Symmetric(X)) in MOI.LogDetConeTriangle(2), +) + # ## Other Cones and Functions # For other cones supported by JuMP, check out the diff --git a/src/print.jl b/src/print.jl index ba03e5dccd7..fd97da8f1c1 100644 --- a/src/print.jl +++ b/src/print.jl @@ -696,6 +696,10 @@ function function_string( return str * "\\end{bmatrix}" end +function function_string(mime::MIME, f::Tuple) + return string("(", join([function_string(mime, fi) for fi in f], ", "), ")") +end + function function_string(mode, constraint::AbstractConstraint) f = reshape_vector(jump_function(constraint), shape(constraint)) return function_string(mode, f) diff --git a/src/sd.jl b/src/sd.jl index d409771d324..dae7be1a2c0 100644 --- a/src/sd.jl +++ b/src/sd.jl @@ -768,3 +768,120 @@ function build_variable( x = _vectorize_variables(_error, variables) return VariablesConstrainedOnCreation(x, set, SymmetricMatrixShape(n)) end + +""" + TupleShape(args::Pair{<:AbstractShape,Int}...) <: AbstractShape + +A shape for representing tuple-valued functions, where each element of the tuple +is a different [`AbstractShape`](@ref) mapped to the number of scalar elements +in that shape. + +## Example + +```jldoctest +julia> shape = TupleShape(ScalarShape() => 1, SymmetricMatrixShape(2) => 3); +``` +""" +struct TupleShape{T} <: AbstractShape + shapes::T + TupleShape(args::Pair{<:AbstractShape,Int}...) = new{typeof(args)}(args) +end + +function vectorize(f::Tuple, shape::TupleShape) + vectors = [vectorize(fi, si[1]) for (fi, si) in zip(f, shape.shapes)] + return reduce(vcat, vectors) +end + +reshape_set(set::MOI.AbstractVectorSet, ::TupleShape) = set + +function reshape_vector(v::Vector{T}, shape::TupleShape) where {T} + out = Any[] + offset = 0 + for (s, dimension) in shape.shapes + push!(out, reshape_vector(v[offset.+(1:dimension)], s)) + offset += dimension + end + return tuple(out...) +end + +_shape_from_function(::Union{Real,AbstractJuMPScalar}) = ScalarShape() => 1 +_shape_from_function(x::AbstractVector) = VectorShape() => length(x) + +function build_constraint( + error_fn::Function, + f::Tuple, + set::Union{MOI.AbstractScalarSet,MOI.AbstractVectorSet}, +) + return error_fn( + "The tuple function $(typeof(f)) is not supported for a set of type " * + "$(typeof(set)).", + ) +end + +function build_constraint( + error_fn::Function, + f::Tuple{Vararg{Union{Real,AbstractJuMPScalar,AbstractVector}}}, + set::MOI.AbstractVectorSet, +) + shape = TupleShape(_shape_from_function.(f)...) + return VectorConstraint(vectorize(f, shape), set, shape) +end + +function build_constraint( + error_fn::Function, + f::Tuple{<:Union{Real,AbstractJuMPScalar},<:LinearAlgebra.Symmetric}, + set::MOI.RootDetConeTriangle, +) + n = LinearAlgebra.checksquare(f[2]) + shape = TupleShape( + ScalarShape() => 1, + SymmetricMatrixShape(n) => div(n * (n + 1), 2), + ) + return VectorConstraint(vectorize(f, shape), set, shape) +end + +function build_constraint( + error_fn::Function, + f::Tuple{<:Union{Real,AbstractJuMPScalar},<:AbstractMatrix}, + set::MOI.RootDetConeSquare, +) + n = LinearAlgebra.checksquare(f[2]) + shape = TupleShape(ScalarShape() => 1, SquareMatrixShape(n) => n^2) + return VectorConstraint(vectorize(f, shape), set, shape) +end + +function build_constraint( + error_fn::Function, + f::Tuple{ + <:Union{Real,AbstractJuMPScalar}, + <:Union{Real,AbstractJuMPScalar}, + <:LinearAlgebra.Symmetric, + }, + set::MOI.LogDetConeTriangle, +) + n = LinearAlgebra.checksquare(f[3]) + shape = TupleShape( + ScalarShape() => 1, + ScalarShape() => 1, + SymmetricMatrixShape(n) => div(n * (n + 1), 2), + ) + return VectorConstraint(vectorize(f, shape), set, shape) +end + +function build_constraint( + error_fn::Function, + f::Tuple{ + <:Union{Real,AbstractJuMPScalar}, + <:Union{Real,AbstractJuMPScalar}, + <:AbstractMatrix, + }, + set::MOI.LogDetConeSquare, +) + n = LinearAlgebra.checksquare(f[3]) + shape = TupleShape( + ScalarShape() => 1, + ScalarShape() => 1, + SquareMatrixShape(n) => n^2, + ) + return VectorConstraint(vectorize(f, shape), set, shape) +end diff --git a/src/sets.jl b/src/sets.jl index e5651338fdf..a283cd7a782 100644 --- a/src/sets.jl +++ b/src/sets.jl @@ -41,6 +41,16 @@ function build_constraint( return build_constraint(_error, func, moi_set(set, length(func))) end +function build_constraint( + error_fn::Function, + f::Tuple{Vararg{Union{Real,AbstractJuMPScalar,AbstractVector},N}}, + set::AbstractVectorSet, +) where {N} + shape = TupleShape(_shape_from_function.(f)...) + args = vectorize(f, shape) + return VectorConstraint(args, moi_set(set, length(args)), shape) +end + """ build_constraint( _error::Function, diff --git a/src/shapes.jl b/src/shapes.jl index 62d1472db70..1e0693f04ad 100644 --- a/src/shapes.jl +++ b/src/shapes.jl @@ -112,6 +112,8 @@ Shape of scalar constraints. """ struct ScalarShape <: AbstractShape end reshape_vector(α, ::ScalarShape) = α +reshape_vector(x::Vector, ::ScalarShape) = first(x) +vectorize(x, ::ScalarShape) = x """ VectorShape diff --git a/test/test_constraint.jl b/test/test_constraint.jl index 16f256baade..44b2d3dc3c2 100644 --- a/test/test_constraint.jl +++ b/test/test_constraint.jl @@ -1710,4 +1710,93 @@ function test_eltype_for_constraint_primal_complex_float64() return end +function test_tuple_shape_vector() + model = Model() + @variable(model, t) + @variable(model, x[1:3]) + @variable(model, u) + @variable(model, v) + @variable(model, w) + for set in ( + MOI.ExponentialCone(), + MOI.DualExponentialCone(), + MOI.PowerCone(0.5), + MOI.DualPowerCone(0.5), + ) + c = @constraint(model, (u, v, w) in set) + obj = constraint_object(c) + @test obj.func == [u, v, w] + @test reshape_vector(obj.func, obj.shape) == (u, v, w) + @test reshape_set(obj.set, obj.shape) == obj.set + end + for set in + (MOI.SecondOrderCone(4), SecondOrderCone(), MOI.GeometricMeanCone(4)) + c = @constraint(model, (t, x) in set) + obj = constraint_object(c) + @test obj.func == [t; x] + @test reshape_vector(obj.func, obj.shape) == (t, x) + @test reshape_set(obj.set, obj.shape) == obj.set + end + for set in (MOI.RotatedSecondOrderCone(5), RotatedSecondOrderCone()) + c = @constraint(model, (t, u, x) in set) + obj = constraint_object(c) + @test obj.func == [t; u; x] + @test reshape_vector(obj.func, obj.shape) == (t, u, x) + @test reshape_set(obj.set, obj.shape) == obj.set + end + return +end + +function test_tuple_shape_det_cone() + model = Model() + @variable(model, t) + @variable(model, u) + @variable(model, X[1:2, 1:2]) + Y = LinearAlgebra.Symmetric(X) + U = [X[i, j] for j in 1:2 for i in 1:j] + c = @constraint(model, (t, X) in MOI.RootDetConeSquare(2)) + obj = constraint_object(c) + @test obj.func == [t; vec(X)] + @test reshape_vector(obj.func, obj.shape) == (t, X) + @test reshape_set(obj.set, obj.shape) == obj.set + c = @constraint(model, (t, Y) in MOI.RootDetConeTriangle(2)) + obj = constraint_object(c) + @test obj.func == [t; U] + @test reshape_vector(obj.func, obj.shape) == (t, Y) + @test reshape_set(obj.set, obj.shape) == obj.set + c = @constraint(model, (t, u, X) in MOI.LogDetConeSquare(2)) + obj = constraint_object(c) + @test obj.func == [t; u; vec(X)] + @test reshape_vector(obj.func, obj.shape) == (t, u, X) + @test reshape_set(obj.set, obj.shape) == obj.set + c = @constraint(model, (t, u, Y) in MOI.LogDetConeTriangle(2)) + obj = constraint_object(c) + @test obj.func == [t; u; U] + @test reshape_vector(obj.func, obj.shape) == (t, u, Y) + @test reshape_set(obj.set, obj.shape) == obj.set + return +end + +function test_tuple_shape_unsupported_error() + model = Model() + @variable(model, t) + @variable(model, x[1:2]) + err = ErrorException( + "In `@constraint(model, (t,) in MOI.ZeroOne())`: " * + "The tuple function $(Tuple{VariableRef}) is not supported for " * + "a set of type $(MOI.ZeroOne).", + ) + @test_throws_strip err @constraint(model, (t,) in MOI.ZeroOne()) + err = ErrorException( + "In `@constraint(model, ((x[1], x[2]),) in MOI.SOS1([1.0, 2.0]))`: " * + "The tuple function $(Tuple{Tuple{VariableRef, VariableRef}}) is not " * + "supported for a set of type $(MOI.SOS1{Float64}).", + ) + @test_throws_strip( + err, + @constraint(model, ((x[1], x[2]),) in MOI.SOS1([1.0, 2.0])), + ) + return +end + end diff --git a/test/test_print.jl b/test/test_print.jl index 790db8a1c40..77849b44fec 100644 --- a/test/test_print.jl +++ b/test/test_print.jl @@ -877,6 +877,20 @@ function test_print_hermitian_psd_cone() return end +function test_print_rootdetcone() + model = Model() + @variable(model, t) + @variable(model, X[1:2, 1:2], Symmetric) + in_sym = JuMP._math_symbol(MIME("text/plain"), :in) + set = MOI.RootDetConeTriangle(2) + c = @constraint(model, (t, X) in set) + @test sprint(io -> show(io, MIME("text/plain"), c)) == + "(t, [X[1,1] X[1,2];\n X[1,2] X[2,2]]) $in_sym $set" + @test sprint(io -> show(io, MIME("text/latex"), c)) == + "\$\$ (t, \\begin{bmatrix}\nX_{1,1} & X_{1,2}\\\\\n\\cdot & X_{2,2}\\\\\n\\end{bmatrix}) \\in \\text{$set} \$\$" + return +end + function test_print_complex_string_round() @test JuMP._string_round(1.0 + 0.0im) == "1" @test JuMP._string_round(-1.0 + 0.0im) == "-1"