From 32112bed7f0e1c665764b81e478c4df3cfbf98b0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Sun, 17 Sep 2023 01:09:09 +0200 Subject: [PATCH 01/10] Polynomial to QCQP meta-solver --- Project.toml | 2 ++ src/PolyJuMP.jl | 1 + src/QCQP/QCQP.jl | 38 ++++++++++++++++++++++++++++++++++++++ test/qcqp.jl | 47 +++++++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 88 insertions(+) create mode 100644 src/QCQP/QCQP.jl create mode 100644 test/qcqp.jl diff --git a/Project.toml b/Project.toml index 1bea461..487064b 100644 --- a/Project.toml +++ b/Project.toml @@ -4,6 +4,7 @@ repo = "https://github.com/jump-dev/PolyJuMP.jl.git" version = "0.7.0" [deps] +DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" DynamicPolynomials = "7c1d4256-1411-5781-91ec-d7bc3513ac07" JuMP = "4076af6c-e467-56ae-b986-b466b2749572" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" @@ -15,6 +16,7 @@ MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0" SemialgebraicSets = "8e049039-38e8-557d-ae3a-bc521ccf6204" [compat] +DataStructures = "0.18" DynamicPolynomials = "0.5" JuMP = "1" MathOptInterface = "1" diff --git a/src/PolyJuMP.jl b/src/PolyJuMP.jl index 31bb5bb..8b5452e 100644 --- a/src/PolyJuMP.jl +++ b/src/PolyJuMP.jl @@ -33,5 +33,6 @@ include("default.jl") include("model.jl") include("KKT/KKT.jl") +include("QCQP/QCQP.jl") end # module diff --git a/src/QCQP/QCQP.jl b/src/QCQP/QCQP.jl new file mode 100644 index 0000000..a49e35e --- /dev/null +++ b/src/QCQP/QCQP.jl @@ -0,0 +1,38 @@ +module QCQP + +import MultivariatePolynomials as MP +import DataStructures + +function decompose(monos::AbstractVector{M}) where {M<:MP.AbstractMonomial} + vars = MP.variables(monos) + quad = DataStructures.OrderedDict{M,Union{Nothing,M}}(var => var for var in vars) + for mono in monos + quad[mono] = nothing + end + candidates = DataStructures.PriorityQueue{eltype(monos),Int}(Base.Order.Reverse) + while any(mono -> isnothing(quad[mono]), keys(quad)) + empty!(candidates) + for mono in keys(quad) + if !isnothing(quad[mono]) + continue + end + for a in keys(quad) + if a != mono && MP.divides(a, mono) + b = MP.div_multiple(mono, a) + if b in keys(quad) + quad[mono] = a + else + candidates[b] = get(candidates, b, 0) + 1 + end + end + end + end + if !isempty(candidates) + next, _ = first(candidates) + quad[next] = nothing + end + end + return quad +end + +end diff --git a/test/qcqp.jl b/test/qcqp.jl new file mode 100644 index 0000000..2921976 --- /dev/null +++ b/test/qcqp.jl @@ -0,0 +1,47 @@ +module TestQCQP + +using Test + +import MultivariatePolynomials as MP +import PolyJuMP + +function _test_decompose(monos, exps) + vars = MP.variables(monos) + M = eltype(monos) + expected = PolyJuMP.QCQP.DataStructures.OrderedDict{M,M}(var => var for var in vars) + for exp in exps + expected[exp[1]] = exp[2] + end + quad = PolyJuMP.QCQP.decompose(monos) + @test quad == expected +end + +function test_decompose(x, y, _) + _test_decompose([x * y], [x * y => y]) + _test_decompose([x^2, y^3, x^2 * y^3], [x^2 => x, y^2 => y, x^2 * y^3 => y^3, y^3 => y^2]) +end + +function runtests(x, y, T) + for name in names(@__MODULE__; all = true) + if startswith("$name", "test_") + @testset "$(name)" begin + getfield(@__MODULE__, name)(x, y, T) + end + end + end +end + +end # module + +using Test + +import DynamicPolynomials +@testset "DynamicPolynomials" begin + DynamicPolynomials.@polyvar(x, y) + TestQCQP.runtests(x, y, Float64) +end +import TypedPolynomials +@testset "TypedPolynomials" begin + TypedPolynomials.@polyvar(x, y) + TestQCQP.runtests(x, y, Float64) +end From 9d3338a5289be950bbc4a30d0bd8a802ab0f9254 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Sun, 17 Sep 2023 21:50:53 +0200 Subject: [PATCH 02/10] Fix format --- src/QCQP/QCQP.jl | 8 ++++++-- test/qcqp.jl | 9 +++++++-- 2 files changed, 13 insertions(+), 4 deletions(-) diff --git a/src/QCQP/QCQP.jl b/src/QCQP/QCQP.jl index a49e35e..f6ffa60 100644 --- a/src/QCQP/QCQP.jl +++ b/src/QCQP/QCQP.jl @@ -3,13 +3,17 @@ module QCQP import MultivariatePolynomials as MP import DataStructures +# TODO special case for squares ? function decompose(monos::AbstractVector{M}) where {M<:MP.AbstractMonomial} vars = MP.variables(monos) - quad = DataStructures.OrderedDict{M,Union{Nothing,M}}(var => var for var in vars) + quad = DataStructures.OrderedDict{M,Union{Nothing,M}}( + var => var for var in vars + ) for mono in monos quad[mono] = nothing end - candidates = DataStructures.PriorityQueue{eltype(monos),Int}(Base.Order.Reverse) + candidates = + DataStructures.PriorityQueue{eltype(monos),Int}(Base.Order.Reverse) while any(mono -> isnothing(quad[mono]), keys(quad)) empty!(candidates) for mono in keys(quad) diff --git a/test/qcqp.jl b/test/qcqp.jl index 2921976..7d7a252 100644 --- a/test/qcqp.jl +++ b/test/qcqp.jl @@ -8,7 +8,9 @@ import PolyJuMP function _test_decompose(monos, exps) vars = MP.variables(monos) M = eltype(monos) - expected = PolyJuMP.QCQP.DataStructures.OrderedDict{M,M}(var => var for var in vars) + expected = PolyJuMP.QCQP.DataStructures.OrderedDict{M,M}( + var => var for var in vars + ) for exp in exps expected[exp[1]] = exp[2] end @@ -18,7 +20,10 @@ end function test_decompose(x, y, _) _test_decompose([x * y], [x * y => y]) - _test_decompose([x^2, y^3, x^2 * y^3], [x^2 => x, y^2 => y, x^2 * y^3 => y^3, y^3 => y^2]) + return _test_decompose( + [x^2, y^3, x^2 * y^3], + [x^2 => x, y^2 => y, x^2 * y^3 => y^3, y^3 => y^2], + ) end function runtests(x, y, T) From 3773423c62fda2e80895d403faa35db8fc97d1b5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Tue, 19 Sep 2023 10:27:16 +0200 Subject: [PATCH 03/10] Add MOI wrapper --- src/QCQP/MOI_wrapper.jl | 55 +++++++++++++++++++++++++++++++++++++++++ src/QCQP/QCQP.jl | 3 +++ 2 files changed, 58 insertions(+) create mode 100644 src/QCQP/MOI_wrapper.jl diff --git a/src/QCQP/MOI_wrapper.jl b/src/QCQP/MOI_wrapper.jl new file mode 100644 index 0000000..e180dbe --- /dev/null +++ b/src/QCQP/MOI_wrapper.jl @@ -0,0 +1,55 @@ +import MathOptInterface as MOI + +struct Optimizer{T,O<:MOI.ModelLike} <: MOI.AbstractOptimizer + model::O +end +function Optimizer(model::MOI.ModelLike) + return Optimizer{Float64,typeof(model)}(model) +end + +function MOI.get( + model::Optimizer{T}, + attr::MOI.Bridges.ListOfNonstandardBridges, +) where {T} + list = copy(MOI.get(model.model, attr)) + push!(list, PolyJuMP.Bridges.Constraint.ToPolynomialBridge{T}) + push!(list, PolyJuMP.Bridges.Objective.ToPolynomialBridge{T}) + return list +end + +MOI.is_empty(model::Optimizer) = MOI.is_empty(model.model) +MOI.empty!(model::Optimizer) = MOI.empty!(model.model) + +MOI.add_variable(model::Optimizer) = MOI.add_variable(model.model) + +function MOI.supports_add_constrained_variable(model::Optimizer, ::Type{S}) where {S<:MOI.AbstractScalarSet} + return MOI.supports_add_constrained_variable(model.model, S) +end + +function MOI.supports_add_constrained_variables(model::Optimizer, ::Type{MOI.Reals}) + return MOI.supports_add_constrained_variables(model.model, MOI.Reals) +end + +function MOI.supports_add_constrained_variables(model::Optimizer, ::Type{S}) where {S<:MOI.AbstractVectorSet} + return MOI.supports_add_constrained_variables(model.model, S) +end + +function MOI.supports(model::Optimizer, attr::MOI.ObjectiveFunction) + return MOI.supports(model.model, attr) +end + +function MOI.supports_constraint(model::Optimizer, ::Type{F}, ::Type{S}) where {F<:MOI.AbstractFunction,S<:MOI.AbstractSet} + return MOI.supports_constraint(model.model, F, S) +end + +function MOI.supports_incremental_interface(model::Optimizer) + return MOI.supports_incremental_interface(model.model) +end +function MOI.copy_to(dest::Optimizer, src::MOI.ModelLike) + return MOI.Utilities.default_copy_to(dest, src) +end + +function MOI.get(model::Optimizer, attr::MOI.SolverName) + name = MOI.get(model.model, attr) + return "PolyJuMP.QCQP with $name" +end diff --git a/src/QCQP/QCQP.jl b/src/QCQP/QCQP.jl index f6ffa60..37519f0 100644 --- a/src/QCQP/QCQP.jl +++ b/src/QCQP/QCQP.jl @@ -2,6 +2,7 @@ module QCQP import MultivariatePolynomials as MP import DataStructures +import PolyJuMP # TODO special case for squares ? function decompose(monos::AbstractVector{M}) where {M<:MP.AbstractMonomial} @@ -39,4 +40,6 @@ function decompose(monos::AbstractVector{M}) where {M<:MP.AbstractMonomial} return quad end +include("MOI_wrapper.jl") + end From 272d2bea44acb8b6340090d03bf591731c6e8e2f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Tue, 19 Sep 2023 10:30:04 +0200 Subject: [PATCH 04/10] Fix format --- src/QCQP/MOI_wrapper.jl | 21 +++++++++++++++++---- 1 file changed, 17 insertions(+), 4 deletions(-) diff --git a/src/QCQP/MOI_wrapper.jl b/src/QCQP/MOI_wrapper.jl index e180dbe..2ee8caf 100644 --- a/src/QCQP/MOI_wrapper.jl +++ b/src/QCQP/MOI_wrapper.jl @@ -22,15 +22,24 @@ MOI.empty!(model::Optimizer) = MOI.empty!(model.model) MOI.add_variable(model::Optimizer) = MOI.add_variable(model.model) -function MOI.supports_add_constrained_variable(model::Optimizer, ::Type{S}) where {S<:MOI.AbstractScalarSet} +function MOI.supports_add_constrained_variable( + model::Optimizer, + ::Type{S}, +) where {S<:MOI.AbstractScalarSet} return MOI.supports_add_constrained_variable(model.model, S) end -function MOI.supports_add_constrained_variables(model::Optimizer, ::Type{MOI.Reals}) +function MOI.supports_add_constrained_variables( + model::Optimizer, + ::Type{MOI.Reals}, +) return MOI.supports_add_constrained_variables(model.model, MOI.Reals) end -function MOI.supports_add_constrained_variables(model::Optimizer, ::Type{S}) where {S<:MOI.AbstractVectorSet} +function MOI.supports_add_constrained_variables( + model::Optimizer, + ::Type{S}, +) where {S<:MOI.AbstractVectorSet} return MOI.supports_add_constrained_variables(model.model, S) end @@ -38,7 +47,11 @@ function MOI.supports(model::Optimizer, attr::MOI.ObjectiveFunction) return MOI.supports(model.model, attr) end -function MOI.supports_constraint(model::Optimizer, ::Type{F}, ::Type{S}) where {F<:MOI.AbstractFunction,S<:MOI.AbstractSet} +function MOI.supports_constraint( + model::Optimizer, + ::Type{F}, + ::Type{S}, +) where {F<:MOI.AbstractFunction,S<:MOI.AbstractSet} return MOI.supports_constraint(model.model, F, S) end From 787ef7cee9a8a14009d54da52bbb02e7325472c8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Sat, 23 Sep 2023 12:21:14 +0200 Subject: [PATCH 05/10] Fixes --- src/QCQP/MOI_wrapper.jl | 19 ++++++++++++++++++- src/functions.jl | 40 ++++++++++++++++++++++++++++++---------- src/nl_to_polynomial.jl | 7 ++++++- 3 files changed, 54 insertions(+), 12 deletions(-) diff --git a/src/QCQP/MOI_wrapper.jl b/src/QCQP/MOI_wrapper.jl index 2ee8caf..b9362ce 100644 --- a/src/QCQP/MOI_wrapper.jl +++ b/src/QCQP/MOI_wrapper.jl @@ -43,10 +43,19 @@ function MOI.supports_add_constrained_variables( return MOI.supports_add_constrained_variables(model.model, S) end -function MOI.supports(model::Optimizer, attr::MOI.ObjectiveFunction) +function MOI.supports(model::Optimizer, attr::MOI.AbstractModelAttribute) return MOI.supports(model.model, attr) end +function MOI.set(model::Optimizer, attr::MOI.AbstractModelAttribute, value) + MOI.set(model.model, attr, value) +end + +function MOI.get(model::Optimizer, attr::MOI.AbstractModelAttribute) + return MOI.get(model.model, attr) +end + + function MOI.supports_constraint( model::Optimizer, ::Type{F}, @@ -55,6 +64,14 @@ function MOI.supports_constraint( return MOI.supports_constraint(model.model, F, S) end +function MOI.add_constraint( + model::Optimizer, + func::MOI.AbstractFunction, + set::MOI.AbstractSet, +) + return MOI.add_constraint(model.model, func, set) +end + function MOI.supports_incremental_interface(model::Optimizer) return MOI.supports_incremental_interface(model.model) end diff --git a/src/functions.jl b/src/functions.jl index 7fa3245..71350ff 100644 --- a/src/functions.jl +++ b/src/functions.jl @@ -1,6 +1,6 @@ """ struct ScalarPolynomialFunction{P<:MP.AbstractPolynomialLike} <: MOI.AbstractScalarFunction - p::P + polynomial::P variables::Vector{MOI.VariableIndex} end @@ -64,6 +64,22 @@ function Base.convert( return ScalarPolynomialFunction{T,P}(poly, variables) end +function _to_polynomial!( + d::Dict{K,V}, + ::Type{T}, + func::MOI.ScalarQuadraticFunction{T}, +) where {K,V,T} + terms = MP.term_type(V, T)[MOI.constant(func)] + for t in func.affine_terms + push!(terms, MP.term(t.coefficient, _to_polynomial!(d, T, t.variable))) + end + for t in func.quadratic_terms + coef = t.variable_1 == t.variable_2 ? t.coefficient / 2 : t.coefficient + push!(terms, MP.term(coef, _to_polynomial!(d, T, t.variable_1) * _to_polynomial!(d, T, t.variable_2))) + end + return MP.polynomial(terms) +end + function Base.convert( ::Type{ScalarPolynomialFunction{T,P}}, func::MOI.ScalarQuadraticFunction{T}, @@ -73,15 +89,8 @@ function Base.convert( quad_variables_2 = [t.variable_2 for t in func.quadratic_terms] variables = [linear_variables; quad_variables_1; quad_variables_2] _, d = _polynomial_variables!(P, variables) - terms = MP.term_type(P)[MOI.constant(func)] - for t in func.affine_terms - push!(terms, MP.term(t.coefficient, d[t.variable])) - end - for t in func.quadratic_terms - coef = t.variable_1 == t.variable_2 ? t.coefficient / 2 : t.coefficient - push!(terms, MP.term(coef, d[t.variable_1] * d[t.variable_2])) - end - return ScalarPolynomialFunction{T,P}(MP.polynomial(terms), variables) + poly = _to_polynomial!(d, T, func) + return ScalarPolynomialFunction{T,P}(poly, variables) end function Base.convert( @@ -149,3 +158,14 @@ function MOI.Utilities.promote_operation( ) where {T,P} return MOI.VectorQuadraticFunction{T} end + +function MOI.Utilities.operate( + op::Union{typeof(+),typeof(-)}, + ::Type{T}, + p::ScalarPolynomialFunction{T,P}, + f::Union{T,MOI.AbstractScalarFunction}, +) where {T,P} + d = Dict(vi => v for (vi, v) in zip(p.variables, MP.variables(p.polynomial))) + poly = _to_polynomial!(d, T, f) + return _scalar_polynomial(d, T, op(p.polynomial, poly)) +end diff --git a/src/nl_to_polynomial.jl b/src/nl_to_polynomial.jl index 84c0095..22e81eb 100644 --- a/src/nl_to_polynomial.jl +++ b/src/nl_to_polynomial.jl @@ -97,10 +97,15 @@ end function _to_polynomial(expr, ::Type{T}) where {T} d = Dict{MOI.VariableIndex,VarType}() poly = _to_polynomial!(d, T, expr) + return _scalar_polynomial(d, T, poly) +end + +function _scalar_polynomial(d::Dict{K,V}, ::Type{T}, poly) where {T,K,V} variable_map = collect(d) sort!(variable_map, by = x -> x[2]) variables = [x[1] for x in variable_map] - return FuncType{T}(poly, variables) + P = MP.polynomial_type(V, T) + return ScalarPolynomialFunction{T,P}(poly, variables) end function _to_polynomial(model::NLToPolynomial{T}, expr) where {T} From 0ece9cbdb2284f500e902b7861bfe4fe2f6ebdd8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Sun, 24 Sep 2023 12:08:41 +0200 Subject: [PATCH 06/10] Fixes --- src/QCQP/MOI_wrapper.jl | 93 ++++++++++++++++++++++++++++++++++++++++- src/QCQP/QCQP.jl | 1 + src/functions.jl | 29 +++++++++++-- 3 files changed, 117 insertions(+), 6 deletions(-) diff --git a/src/QCQP/MOI_wrapper.jl b/src/QCQP/MOI_wrapper.jl index b9362ce..91639ad 100644 --- a/src/QCQP/MOI_wrapper.jl +++ b/src/QCQP/MOI_wrapper.jl @@ -1,10 +1,12 @@ import MathOptInterface as MOI -struct Optimizer{T,O<:MOI.ModelLike} <: MOI.AbstractOptimizer +mutable struct Optimizer{T,O<:MOI.ModelLike} <: MOI.AbstractOptimizer model::O + objective::Union{Nothing,PolyJuMP.ScalarPolynomialFunction{T}} end + function Optimizer(model::MOI.ModelLike) - return Optimizer{Float64,typeof(model)}(model) + return Optimizer{Float64,typeof(model)}(model, nothing) end function MOI.get( @@ -47,6 +49,22 @@ function MOI.supports(model::Optimizer, attr::MOI.AbstractModelAttribute) return MOI.supports(model.model, attr) end +function MOI.supports( + ::Optimizer, + ::MOI.ObjectiveFunction{<:PolyJuMP.ScalarPolynomialFunction}, +) + return true +end + +function MOI.set( + model::Optimizer{T}, + ::MOI.ObjectiveFunction{F}, + f::F, +) where {T,F<:PolyJuMP.ScalarPolynomialFunction{T}} + model.objective = f + return +end + function MOI.set(model::Optimizer, attr::MOI.AbstractModelAttribute, value) MOI.set(model.model, attr, value) end @@ -75,10 +93,81 @@ end function MOI.supports_incremental_interface(model::Optimizer) return MOI.supports_incremental_interface(model.model) end + function MOI.copy_to(dest::Optimizer, src::MOI.ModelLike) return MOI.Utilities.default_copy_to(dest, src) end +MOI.optimize!(model::Optimizer) = MOI.optimize!(model.model) + +function _quad_convert(index, p::MP.AbstractPolynomialLike{T}, div) where {T} + q = zero(MOI.ScalarQuadraticFunction{T}) + for t in MP.terms(p) + α = MP.coefficient(t) + mono = MP.monomial(t) + if MP.degree(mono) == 0 + MA.operate!(+, q, α) + else + vars = MP.effective_variables(mono) + if MP.degree(mono) == 1 + @assert length(vars) == 1 + MA.operate!(MA.add_mul, q, α, index(first(vars))) + elseif MP.degree(mono) == 2 + x = first(vars) + if length(vars) == 1 + y = x + else + @assert length(vars) == 2 + y = vars[2] + end + MA.operate!(MA.add_mul, q, α, index(x), index(y)) + else + x = div[mono] + y = MP.div_multiple(mono, x) + MA.operate!(MA.add_mul, q, α, index(x), index(y)) + end + end + end + return q +end + +function _variable_to_index_map(p::PolyJuMP.ScalarPolynomialFunction{T,P}) where {T,P} + M = MP.monomial_type(P) + return Dict{M,MOI.VariableIndex}( + v => vi for (v, vi) in zip(MP.variables(p.polynomial), p.variables) + ) +end + +function monomial_variable_index(model::Optimizer{T}, d::Dict, div, mono::MP.AbstractMonomialLike) where {T} + if !haskey(d, mono) + d[mono] = MOI.add_variable(model.model) + x = div[mono] + vx = monomial_variable_index(model, d, div, x) + y = MP.div_multiple(mono, x) + vy = monomial_variable_index(model, d, div, y) + MOI.add_constraint(model, + MA.@rewrite(one(T) * d[mono] - one(T) * vx * vy), + MOI.EqualTo(zero(T)), + ) + end + return d[mono] +end + +function MOI.Utilities.final_touch(model::Optimizer{T}, _) where {T} + if !isnothing(model.objective) + d = _variable_to_index_map(model.objective) + p = model.objective.polynomial + monos = MP.monomials(p) + div = decompose(monos) + function index(mono) + return monomial_variable_index(model, d, div, mono) + end + obj = _quad_convert(index, p, div) + MOI.set(model.model, MOI.ObjectiveFunction{typeof(obj)}(), obj) + end + return +end + function MOI.get(model::Optimizer, attr::MOI.SolverName) name = MOI.get(model.model, attr) return "PolyJuMP.QCQP with $name" diff --git a/src/QCQP/QCQP.jl b/src/QCQP/QCQP.jl index 37519f0..4e3bf03 100644 --- a/src/QCQP/QCQP.jl +++ b/src/QCQP/QCQP.jl @@ -1,5 +1,6 @@ module QCQP +import MutableArithmetics as MA import MultivariatePolynomials as MP import DataStructures import PolyJuMP diff --git a/src/functions.jl b/src/functions.jl index 71350ff..83ef50f 100644 --- a/src/functions.jl +++ b/src/functions.jl @@ -1,5 +1,5 @@ """ - struct ScalarPolynomialFunction{P<:MP.AbstractPolynomialLike} <: MOI.AbstractScalarFunction + struct ScalarPolynomialFunction{T,P<:MP.AbstractPolynomial{T}} <: MOI.AbstractScalarFunction polynomial::P variables::Vector{MOI.VariableIndex} end @@ -136,7 +136,7 @@ end function MOI.Utilities.promote_operation( ::typeof(-), ::Type{T}, - F::Type{ScalarPolynomialFunction{T,P}}, + F::Type{<:Union{ScalarPolynomialFunction{T,P},VectorPolynomialFunction{T,P}}}, ) where {T,P} return F end @@ -150,13 +150,34 @@ function MOI.Utilities.promote_operation( return F end -# FIXME +# Placeholder for `promote_operation` +struct VectorPolynomialFunction{T,P<:MP.AbstractPolynomial{T}} <: MOI.AbstractVectorFunction +end + +function MOI.Utilities.scalar_type(::Type{VectorPolynomialFunction{T,P}}) where {T,P} + return PolyJuMP.ScalarPolynomialFunction{T,P} +end + +function MOI.Utilities.is_coefficient_type( + ::Type{<:VectorPolynomialFunction{T}}, + ::Type{T}, +) where {T} + return true +end + +function MOI.Utilities.is_coefficient_type( + ::Type{<:VectorPolynomialFunction}, + ::Type, +) + return false +end + function MOI.Utilities.promote_operation( ::typeof(vcat), ::Type{T}, ::Type{ScalarPolynomialFunction{T,P}}, ) where {T,P} - return MOI.VectorQuadraticFunction{T} + return VectorPolynomialFunction{T,P} end function MOI.Utilities.operate( From 123185d42073f1906dcf2d73182bd4467df4afe7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Sun, 24 Sep 2023 13:35:37 +0200 Subject: [PATCH 07/10] Add test --- src/functions.jl | 43 ++++++++++++++++++++++++++----------------- test/qcqp.jl | 37 +++++++++++++++++++++++++++++++++++++ 2 files changed, 63 insertions(+), 17 deletions(-) diff --git a/src/functions.jl b/src/functions.jl index 83ef50f..975191a 100644 --- a/src/functions.jl +++ b/src/functions.jl @@ -133,23 +133,6 @@ function MOI.Utilities.is_coefficient_type( return S === T end -function MOI.Utilities.promote_operation( - ::typeof(-), - ::Type{T}, - F::Type{<:Union{ScalarPolynomialFunction{T,P},VectorPolynomialFunction{T,P}}}, -) where {T,P} - return F -end - -function MOI.Utilities.promote_operation( - ::typeof(-), - ::Type{T}, - F::Type{ScalarPolynomialFunction{T,P}}, - ::Type{<:Union{T,MOI.Utilities.ScalarLike{T}}}, -) where {T,P} - return F -end - # Placeholder for `promote_operation` struct VectorPolynomialFunction{T,P<:MP.AbstractPolynomial{T}} <: MOI.AbstractVectorFunction end @@ -172,6 +155,32 @@ function MOI.Utilities.is_coefficient_type( return false end +function MOI.Utilities.promote_operation( + ::typeof(-), + ::Type{T}, + F::Type{<:Union{ScalarPolynomialFunction{T,P},VectorPolynomialFunction{T,P}}}, +) where {T,P} + return F +end + +function MOI.Utilities.promote_operation( + ::typeof(-), + ::Type{T}, + F::Type{ScalarPolynomialFunction{T,P}}, + ::Type{<:Union{T,MOI.Utilities.ScalarLike{T}}}, +) where {T,P} + return F +end + +function MOI.Utilities.promote_operation( + ::typeof(-), + ::Type{T}, + F::Type{VectorPolynomialFunction{T,P}}, + ::Type{<:Union{AbstractVector{T},MOI.Utilities.VectorLike{T}}}, +) where {T,P} + return F +end + function MOI.Utilities.promote_operation( ::typeof(vcat), ::Type{T}, diff --git a/test/qcqp.jl b/test/qcqp.jl index 7d7a252..0e066d7 100644 --- a/test/qcqp.jl +++ b/test/qcqp.jl @@ -2,6 +2,7 @@ module TestQCQP using Test +import MathOptInterface as MOI import MultivariatePolynomials as MP import PolyJuMP @@ -26,6 +27,42 @@ function test_decompose(x, y, _) ) end +MOI.Utilities.@model( + Model, + (), + (MOI.LessThan, MOI.GreaterThan, MOI.EqualTo), + (), + (), + (), + (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction), + (), + (), +) + +function MOI.supports( + ::Model, + ::MOI.ObjectiveFunction{MOI.ScalarNonlinearFunction}, +) + return false +end + +function test_objective(x, y, T) + inner = Model{T}() + optimizer = MOI.Utilities.MockOptimizer(inner) + model = PolyJuMP.JuMP.GenericModel{T}(() -> PolyJuMP.QCQP.Optimizer(optimizer)) + PolyJuMP.@variable(model, 0 <= a) + PolyJuMP.@variable(model, 0 <= b) + PolyJuMP.@constraint(model, a + b >= 1) + PolyJuMP.@objective(model, Min, a^3 - a^2 + 2a*b - b^2 + b^3) + PolyJuMP.optimize!(model) + vis = MOI.get(inner, MOI.ListOfVariableIndices()) + @test length(vis) == 4 + F = MOI.ScalarQuadraticFunction{T} + S = MOI.EqualTo{T} + cis = MOI.get(inner, MOI.ListOfConstraintIndices{F,S}()) + @test length(cis) == 2 +end + function runtests(x, y, T) for name in names(@__MODULE__; all = true) if startswith("$name", "test_") From 4bf75beb0864c313e1af8c3674e76cdcf921e245 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Sun, 24 Sep 2023 20:59:41 +0200 Subject: [PATCH 08/10] Add bounds --- src/QCQP/MOI_wrapper.jl | 51 ++++++++++++++++++++++++++--------------- src/nl_to_polynomial.jl | 2 +- test/qcqp.jl | 17 ++++++++++++-- 3 files changed, 48 insertions(+), 22 deletions(-) diff --git a/src/QCQP/MOI_wrapper.jl b/src/QCQP/MOI_wrapper.jl index 91639ad..669e1ee 100644 --- a/src/QCQP/MOI_wrapper.jl +++ b/src/QCQP/MOI_wrapper.jl @@ -22,6 +22,16 @@ end MOI.is_empty(model::Optimizer) = MOI.is_empty(model.model) MOI.empty!(model::Optimizer) = MOI.empty!(model.model) +MOI.is_valid(model::Optimizer, i::MOI.Index) = MOI.is_valid(model.model, i) + +function MOI.get( + model::Optimizer, + attr::MOI.AbstractConstraintAttribute, + ci::MOI.ConstraintIndex, +) + return MOI.get(model.model, attr, ci) +end + MOI.add_variable(model::Optimizer) = MOI.add_variable(model.model) function MOI.supports_add_constrained_variable( @@ -100,7 +110,7 @@ end MOI.optimize!(model::Optimizer) = MOI.optimize!(model.model) -function _quad_convert(index, p::MP.AbstractPolynomialLike{T}, div) where {T} +function _quad_convert(p::MP.AbstractPolynomialLike{T}, index, div) where {T} q = zero(MOI.ScalarQuadraticFunction{T}) for t in MP.terms(p) α = MP.coefficient(t) @@ -108,23 +118,12 @@ function _quad_convert(index, p::MP.AbstractPolynomialLike{T}, div) where {T} if MP.degree(mono) == 0 MA.operate!(+, q, α) else - vars = MP.effective_variables(mono) - if MP.degree(mono) == 1 - @assert length(vars) == 1 - MA.operate!(MA.add_mul, q, α, index(first(vars))) - elseif MP.degree(mono) == 2 - x = first(vars) - if length(vars) == 1 - y = x - else - @assert length(vars) == 2 - y = vars[2] - end - MA.operate!(MA.add_mul, q, α, index(x), index(y)) + if haskey(index, mono) + MA.operate!(MA.add_mul, q, α, index[mono]) else x = div[mono] y = MP.div_multiple(mono, x) - MA.operate!(MA.add_mul, q, α, index(x), index(y)) + MA.operate!(MA.add_mul, q, α, index[x], index[y]) end end end @@ -140,11 +139,19 @@ end function monomial_variable_index(model::Optimizer{T}, d::Dict, div, mono::MP.AbstractMonomialLike) where {T} if !haskey(d, mono) - d[mono] = MOI.add_variable(model.model) x = div[mono] vx = monomial_variable_index(model, d, div, x) y = MP.div_multiple(mono, x) vy = monomial_variable_index(model, d, div, y) + lx, ux = MOI.Utilities.get_bounds(model, T, vx) + ly, uy = MOI.Utilities.get_bounds(model, T, vy) + bounds = (lx * ly, lx * uy, ux * ly, ux * uy) + l = min(bounds...) + if vx == vy + l = max(l, zero(T)) + end + u = max(bounds...) + d[mono], _ = MOI.add_constrained_variable(model.model, MOI.Interval(l, u)) MOI.add_constraint(model, MA.@rewrite(one(T) * d[mono] - one(T) * vx * vy), MOI.EqualTo(zero(T)), @@ -159,10 +166,16 @@ function MOI.Utilities.final_touch(model::Optimizer{T}, _) where {T} p = model.objective.polynomial monos = MP.monomials(p) div = decompose(monos) - function index(mono) - return monomial_variable_index(model, d, div, mono) + for mono in sort(collect(keys(div))) + if haskey(d, mono) + continue + end + a = div[mono] + monomial_variable_index(model, d, div, a) + b = MP.div_multiple(mono, a) + monomial_variable_index(model, d, div, b) end - obj = _quad_convert(index, p, div) + obj = _quad_convert(p, d, div) MOI.set(model.model, MOI.ObjectiveFunction{typeof(obj)}(), obj) end return diff --git a/src/nl_to_polynomial.jl b/src/nl_to_polynomial.jl index 22e81eb..bb799e6 100644 --- a/src/nl_to_polynomial.jl +++ b/src/nl_to_polynomial.jl @@ -102,7 +102,7 @@ end function _scalar_polynomial(d::Dict{K,V}, ::Type{T}, poly) where {T,K,V} variable_map = collect(d) - sort!(variable_map, by = x -> x[2]) + sort!(variable_map, by = x -> x[2], rev = true) variables = [x[1] for x in variable_map] P = MP.polynomial_type(V, T) return ScalarPolynomialFunction{T,P}(poly, variables) diff --git a/test/qcqp.jl b/test/qcqp.jl index 0e066d7..04eab8c 100644 --- a/test/qcqp.jl +++ b/test/qcqp.jl @@ -50,17 +50,30 @@ function test_objective(x, y, T) inner = Model{T}() optimizer = MOI.Utilities.MockOptimizer(inner) model = PolyJuMP.JuMP.GenericModel{T}(() -> PolyJuMP.QCQP.Optimizer(optimizer)) - PolyJuMP.@variable(model, 0 <= a) - PolyJuMP.@variable(model, 0 <= b) + PolyJuMP.@variable(model, 1 <= a <= 2) + PolyJuMP.@variable(model, -5 <= b <= 3) PolyJuMP.@constraint(model, a + b >= 1) PolyJuMP.@objective(model, Min, a^3 - a^2 + 2a*b - b^2 + b^3) PolyJuMP.optimize!(model) vis = MOI.get(inner, MOI.ListOfVariableIndices()) @test length(vis) == 4 + a, b, b2, a2 = vis + @test MOI.Utilities.get_bounds(inner, T, a) == (1, 2) + @test MOI.Utilities.get_bounds(inner, T, b) == (-5, 3) + @test MOI.Utilities.get_bounds(inner, T, a2) == (1, 4) + @test MOI.Utilities.get_bounds(inner, T, b2) == (0, 25) F = MOI.ScalarQuadraticFunction{T} S = MOI.EqualTo{T} cis = MOI.get(inner, MOI.ListOfConstraintIndices{F,S}()) @test length(cis) == 2 + o = one(T) + @test MOI.get(inner, MOI.ConstraintFunction(), cis[1]) ≈ b2 - o * b * b + @test MOI.get(inner, MOI.ConstraintFunction(), cis[2]) ≈ a2 - o * a * a + for ci in cis + @test MOI.get(inner, MOI.ConstraintSet(), ci) == MOI.EqualTo(zero(T)) + end + @test MOI.get(inner, MOI.ObjectiveFunction{F}()) ≈ + -o * a2 + 2o * a * b - o * b2 + o * a * a2 + o * b * b2 end function runtests(x, y, T) From 67981407cf108fc4a27c296638f0bac964bc6946 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Sun, 1 Oct 2023 21:17:31 +0200 Subject: [PATCH 09/10] Add suppor to constraints --- src/QCQP/MOI_wrapper.jl | 186 +++++++++++++++++++++++++++++++++++----- src/QCQP/QCQP.jl | 4 +- test/qcqp.jl | 83 +++++++++++++++--- 3 files changed, 242 insertions(+), 31 deletions(-) diff --git a/src/QCQP/MOI_wrapper.jl b/src/QCQP/MOI_wrapper.jl index 669e1ee..30aefbd 100644 --- a/src/QCQP/MOI_wrapper.jl +++ b/src/QCQP/MOI_wrapper.jl @@ -3,12 +3,19 @@ import MathOptInterface as MOI mutable struct Optimizer{T,O<:MOI.ModelLike} <: MOI.AbstractOptimizer model::O objective::Union{Nothing,PolyJuMP.ScalarPolynomialFunction{T}} + constraints::DataStructures.OrderedDict{Type,Tuple{Type,MOI.Utilities.VectorOfConstraints}} end -function Optimizer(model::MOI.ModelLike) - return Optimizer{Float64,typeof(model)}(model, nothing) +function Optimizer{T}(model::MOI.ModelLike) where {T} + return Optimizer{T,typeof(model)}( + model, + nothing, + DataStructures.OrderedDict{Type,MOI.Utilities.VectorOfConstraints}(), + ) end +Optimizer(model::MOI.ModelLike) = Optimizer{Float64}(model) + function MOI.get( model::Optimizer{T}, attr::MOI.Bridges.ListOfNonstandardBridges, @@ -20,9 +27,24 @@ function MOI.get( end MOI.is_empty(model::Optimizer) = MOI.is_empty(model.model) -MOI.empty!(model::Optimizer) = MOI.empty!(model.model) +function MOI.empty!(model::Optimizer) + MOI.empty!(model.model) + model.objective = nothing + empty!(model.constraints) + return +end MOI.is_valid(model::Optimizer, i::MOI.Index) = MOI.is_valid(model.model, i) +function MOI.is_valid( + model::Optimizer{T}, + ::MOI.ConstraintIndex{ + PolyJuMP.ScalarPolynomialFunction{T}, + S + }, +) where {T,S<:MOI.AbstractScalarSet} + return haskey(model.constraints, S) && + MOI.is_valid(model.constraints[S][2], ci) +end function MOI.get( model::Optimizer, @@ -83,7 +105,6 @@ function MOI.get(model::Optimizer, attr::MOI.AbstractModelAttribute) return MOI.get(model.model, attr) end - function MOI.supports_constraint( model::Optimizer, ::Type{F}, @@ -91,6 +112,17 @@ function MOI.supports_constraint( ) where {F<:MOI.AbstractFunction,S<:MOI.AbstractSet} return MOI.supports_constraint(model.model, F, S) end +function MOI.supports_constraint( + model::Optimizer{T}, + ::Type{<:PolyJuMP.ScalarPolynomialFunction{T}}, + ::Type{S}, +) where {T,S<:MOI.AbstractScalarSet} + return MOI.supports_constraint( + model.model, + MOI.ScalarQuadraticFunction{T}, + S, + ) +end function MOI.add_constraint( model::Optimizer, @@ -99,6 +131,43 @@ function MOI.add_constraint( ) return MOI.add_constraint(model.model, func, set) end +function MOI.add_constraint( + model::Optimizer{T}, + func::PolyJuMP.ScalarPolynomialFunction{T,P}, + set::MOI.AbstractScalarSet, +) where {T,P} + F = typeof(func) + S = typeof(set) + if !haskey(model.constraints, S) + con = MOI.Utilities.VectorOfConstraints{F,S}() + model.constraints[S] = (P,con) + end + return MOI.add_constraint(model.constraints[S][2], func, set) +end + +function MOI.get( + model::Optimizer{T}, + attr::Union{ + MOI.ConstraintFunction, + MOI.ConstraintSet, + }, + ci::MOI.ConstraintIndex{ + <:PolyJuMP.ScalarPolynomialFunction{T}, + S, + }, +) where {T,S} + return MOI.get(model.constraints[S][2], attr, ci) +end + +function MOI.get( + model::Optimizer{T}, + attr::MOI.ListOfConstraintIndices{ + <:PolyJuMP.ScalarPolynomialFunction{T}, + S, + } +) where {T,S<:MOI.AbstractScalarSet} + return MOI.get(model.constraints[S][2], attr) +end function MOI.supports_incremental_interface(model::Optimizer) return MOI.supports_incremental_interface(model.model) @@ -130,13 +199,56 @@ function _quad_convert(p::MP.AbstractPolynomialLike{T}, index, div) where {T} return q end -function _variable_to_index_map(p::PolyJuMP.ScalarPolynomialFunction{T,P}) where {T,P} - M = MP.monomial_type(P) - return Dict{M,MOI.VariableIndex}( - v => vi for (v, vi) in zip(MP.variables(p.polynomial), p.variables) +function _add_monomials!(p::PolyJuMP.ScalarPolynomialFunction, monos1) + monos2 = MP.monomials(p.polynomial) + if isnothing(monos1) + return monos2 + else + return MP.merge_monomial_vectors([monos1, monos2]) + end +end + +function _subs!(p::PolyJuMP.ScalarPolynomialFunction{T,P}, ::Nothing) where {T,P} + return p, Dict{MOI.VariableIndex,MP.variable_union_type(P)}( + vi => var for (vi, var) in zip(p.variables, MP.variables(p.polynomial)) ) end +function _subs!(p::PolyJuMP.ScalarPolynomialFunction, index_to_var::Dict{K,V}) where {K,V} + old_var = V[] + new_var = V[] + for (vi, var) in zip(p.variables, MP.variables(p.polynomial)) + if haskey(index_to_var, vi) + if var != index_to_var[vi] + push!(old_var, var) + push!(new_var, index_to_var[vi]) + end + else + index_to_var[vi] = var + end + end + if !isempty(old_var) + poly = MP.subs(p.polynomial, old_var => new_var) + p = PolyJuMP.ScalarPolynomialFunction(poly, p.variables) + end + return p, index_to_var +end + +function _add_variables!(p::PolyJuMP.ScalarPolynomialFunction{T,P}, d) where {T,P} + if isnothing(d) + d = Dict{MP.monomial_type(P),MOI.VariableIndex}() + else + M = promote_type(keytype(d), MP.monomial_type(P)) + if keytype(d) !== M + d = convert(Dict{M,MOI.VariableIndex}, d) + end + end + for (v, vi) in zip(MP.variables(p.polynomial), p.variables) + d[v] = vi + end + return d +end + function monomial_variable_index(model::Optimizer{T}, d::Dict, div, mono::MP.AbstractMonomialLike) where {T} if !haskey(d, mono) x = div[mono] @@ -160,24 +272,58 @@ function monomial_variable_index(model::Optimizer{T}, d::Dict, div, mono::MP.Abs return d[mono] end +function _add_constraints(model::Optimizer, cis, index_to_var, d, div) + for ci in cis + func = MOI.get(model, MOI.ConstraintFunction(), ci) + set = MOI.get(model, MOI.ConstraintSet(), ci) + func, index_to_var = _subs!(func, index_to_var) + quad = _quad_convert(func.polynomial, d, div) + MOI.add_constraint(model, quad, set) + end +end + function MOI.Utilities.final_touch(model::Optimizer{T}, _) where {T} + index_to_var = nothing + vars = nothing + monos = nothing if !isnothing(model.objective) - d = _variable_to_index_map(model.objective) - p = model.objective.polynomial - monos = MP.monomials(p) - div = decompose(monos) - for mono in sort(collect(keys(div))) - if haskey(d, mono) - continue + func, index_to_var = _subs!(model.objective, index_to_var) + vars = _add_variables!(func, vars) + monos = _add_monomials!(func, monos) + end + if !isempty(model.constraints) + for S in keys(model.constraints) + for ci in MOI.get(model, MOI.ListOfConstraintIndices{ + PolyJuMP.ScalarPolynomialFunction{T,model.constraints[S][1]}, + S + }()) + func = MOI.get(model, MOI.ConstraintFunction(), ci) + func, index_to_var = _subs!(func, index_to_var) + vars = _add_variables!(func, vars) + monos = _add_monomials!(func, monos) end - a = div[mono] - monomial_variable_index(model, d, div, a) - b = MP.div_multiple(mono, a) - monomial_variable_index(model, d, div, b) end - obj = _quad_convert(p, d, div) + end + div = decompose(monos) + for mono in sort(collect(keys(div))) + if haskey(vars, mono) + continue + end + a = div[mono] + monomial_variable_index(model, vars, div, a) + b = MP.div_multiple(mono, a) + monomial_variable_index(model, vars, div, b) + end + if !isnothing(model.objective) + func, index_to_var = _subs!(model.objective, index_to_var) + obj = _quad_convert(func.polynomial, vars, div) MOI.set(model.model, MOI.ObjectiveFunction{typeof(obj)}(), obj) end + for S in keys(model.constraints) + F = PolyJuMP.ScalarPolynomialFunction{T,model.constraints[S][1]} + cis = MOI.get(model, MOI.ListOfConstraintIndices{F,S}()) + _add_constraints(model, cis, index_to_var, vars, div) + end return end diff --git a/src/QCQP/QCQP.jl b/src/QCQP/QCQP.jl index 4e3bf03..47fbc67 100644 --- a/src/QCQP/QCQP.jl +++ b/src/QCQP/QCQP.jl @@ -12,7 +12,9 @@ function decompose(monos::AbstractVector{M}) where {M<:MP.AbstractMonomial} var => var for var in vars ) for mono in monos - quad[mono] = nothing + if MP.degree(mono) > 1 + quad[mono] = nothing + end end candidates = DataStructures.PriorityQueue{eltype(monos),Int}(Base.Order.Reverse) diff --git a/test/qcqp.jl b/test/qcqp.jl index 04eab8c..a3afbbd 100644 --- a/test/qcqp.jl +++ b/test/qcqp.jl @@ -30,7 +30,7 @@ end MOI.Utilities.@model( Model, (), - (MOI.LessThan, MOI.GreaterThan, MOI.EqualTo), + (MOI.LessThan, MOI.GreaterThan, MOI.EqualTo, MOI.Interval), (), (), (), @@ -46,14 +46,18 @@ function MOI.supports( return false end -function test_objective(x, y, T) +function _test_objective_or_constraint(x, y, T, obj::Bool) inner = Model{T}() optimizer = MOI.Utilities.MockOptimizer(inner) - model = PolyJuMP.JuMP.GenericModel{T}(() -> PolyJuMP.QCQP.Optimizer(optimizer)) + model = PolyJuMP.JuMP.GenericModel{T}(() -> PolyJuMP.QCQP.Optimizer{T}(optimizer)) PolyJuMP.@variable(model, 1 <= a <= 2) PolyJuMP.@variable(model, -5 <= b <= 3) PolyJuMP.@constraint(model, a + b >= 1) - PolyJuMP.@objective(model, Min, a^3 - a^2 + 2a*b - b^2 + b^3) + if obj + PolyJuMP.@objective(model, Min, a^3 - a^2 + 2a*b - b^2 + b^3) + else + PolyJuMP.@constraint(model, 0 <= a^3 - a^2 + 2a*b - b^2 + b^3 <= 1) + end PolyJuMP.optimize!(model) vis = MOI.get(inner, MOI.ListOfVariableIndices()) @test length(vis) == 4 @@ -72,14 +76,73 @@ function test_objective(x, y, T) for ci in cis @test MOI.get(inner, MOI.ConstraintSet(), ci) == MOI.EqualTo(zero(T)) end - @test MOI.get(inner, MOI.ObjectiveFunction{F}()) ≈ - -o * a2 + 2o * a * b - o * b2 + o * a * a2 + o * b * b2 + exp = -o * a2 + 2o * a * b - o * b2 + o * a * a2 + o * b * b2 + if obj + @test MOI.get(inner, MOI.ObjectiveFunction{F}()) ≈ exp + else + S = MOI.Interval{T} + cis = MOI.get(inner, MOI.ListOfConstraintIndices{F,S}()) + @test length(cis) == 1 + @test MOI.get(inner, MOI.ConstraintFunction(), cis[1]) ≈ exp + end +end + +function test_objective(x, y, T) + return _test_objective_or_constraint(x, y, T, true) +end + +function test_constraint(x, y, T) + return _test_objective_or_constraint(x, y, T, false) +end + +function test_objective_and_constraint(x, y, T) + inner = Model{T}() + optimizer = MOI.Utilities.MockOptimizer(inner) + model = PolyJuMP.JuMP.GenericModel{T}(() -> PolyJuMP.QCQP.Optimizer{T}(optimizer)) + PolyJuMP.@variable(model, -2 <= a <= 3) + PolyJuMP.@variable(model, 5 <= b <= 7) + PolyJuMP.@constraint(model, 0 <= a^3 <= 1) + PolyJuMP.@constraint(model, 0 <= b^3 <= 1) + PolyJuMP.@constraint(model, 0 <= a^3 * b^3 + a^6 <= 1) + PolyJuMP.@objective(model, Max, a^6 * b^3) + PolyJuMP.optimize!(model) + vis = MOI.get(inner, MOI.ListOfVariableIndices()) + @test length(vis) == 7 + a, b, b2, a2, a3, b3, a6 = vis + @test MOI.Utilities.get_bounds(inner, T, a) == (-2, 3) + @test MOI.Utilities.get_bounds(inner, T, b) == (5, 7) + @test MOI.Utilities.get_bounds(inner, T, b2) == (25, 49) + @test MOI.Utilities.get_bounds(inner, T, a2) == (0, 9) + @test MOI.Utilities.get_bounds(inner, T, a3) == (-18, 27) + @test MOI.Utilities.get_bounds(inner, T, b3) == (125, 343) + @test MOI.Utilities.get_bounds(inner, T, a6) == (0, 729) + F = MOI.ScalarQuadraticFunction{T} + S = MOI.EqualTo{T} + cis = MOI.get(inner, MOI.ListOfConstraintIndices{F,S}()) + @test length(cis) == 5 + o = one(T) + z = zero(T) + @test MOI.get(inner, MOI.ConstraintFunction(), cis[1]) ≈ b2 - o * b * b + @test MOI.get(inner, MOI.ConstraintFunction(), cis[2]) ≈ a2 - o * a * a + @test MOI.get(inner, MOI.ConstraintFunction(), cis[3]) ≈ a3 - o * a2 * a + @test MOI.get(inner, MOI.ConstraintFunction(), cis[4]) ≈ b3 - o * b2 * b + @test MOI.get(inner, MOI.ConstraintFunction(), cis[5]) ≈ a6 - o * a3 * a3 + for ci in cis + @test MOI.get(inner, MOI.ConstraintSet(), ci) == MOI.EqualTo(zero(T)) + end + @test MOI.get(inner, MOI.ObjectiveFunction{F}()) ≈ o * a6 * b3 + S = MOI.Interval{T} + cis = MOI.get(inner, MOI.ListOfConstraintIndices{F,S}()) + @test length(cis) == 3 + @test MOI.get(inner, MOI.ConstraintFunction(), cis[1]) ≈ o * a3 + z * a * b + @test MOI.get(inner, MOI.ConstraintFunction(), cis[2]) ≈ o * b3 + z * a * b + @test MOI.get(inner, MOI.ConstraintFunction(), cis[3]) ≈ o * a3 * b3 + o * a6 end -function runtests(x, y, T) +function runtests(x, y) for name in names(@__MODULE__; all = true) if startswith("$name", "test_") - @testset "$(name)" begin + @testset "$(name) $T" for T in [Int, Float64] getfield(@__MODULE__, name)(x, y, T) end end @@ -93,10 +156,10 @@ using Test import DynamicPolynomials @testset "DynamicPolynomials" begin DynamicPolynomials.@polyvar(x, y) - TestQCQP.runtests(x, y, Float64) + TestQCQP.runtests(x, y) end import TypedPolynomials @testset "TypedPolynomials" begin TypedPolynomials.@polyvar(x, y) - TestQCQP.runtests(x, y, Float64) + TestQCQP.runtests(x, y) end From 1ea12ac054733fbbc70e16e9736e870023365caa Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Sun, 1 Oct 2023 21:19:55 +0200 Subject: [PATCH 10/10] Fix format --- src/QCQP/MOI_wrapper.jl | 76 ++++++++++++++++++++++++----------------- src/functions.jl | 25 ++++++++++---- test/qcqp.jl | 15 +++++--- 3 files changed, 74 insertions(+), 42 deletions(-) diff --git a/src/QCQP/MOI_wrapper.jl b/src/QCQP/MOI_wrapper.jl index 30aefbd..9e6d25d 100644 --- a/src/QCQP/MOI_wrapper.jl +++ b/src/QCQP/MOI_wrapper.jl @@ -3,7 +3,10 @@ import MathOptInterface as MOI mutable struct Optimizer{T,O<:MOI.ModelLike} <: MOI.AbstractOptimizer model::O objective::Union{Nothing,PolyJuMP.ScalarPolynomialFunction{T}} - constraints::DataStructures.OrderedDict{Type,Tuple{Type,MOI.Utilities.VectorOfConstraints}} + constraints::DataStructures.OrderedDict{ + Type, + Tuple{Type,MOI.Utilities.VectorOfConstraints}, + } end function Optimizer{T}(model::MOI.ModelLike) where {T} @@ -37,13 +40,10 @@ end MOI.is_valid(model::Optimizer, i::MOI.Index) = MOI.is_valid(model.model, i) function MOI.is_valid( model::Optimizer{T}, - ::MOI.ConstraintIndex{ - PolyJuMP.ScalarPolynomialFunction{T}, - S - }, + ::MOI.ConstraintIndex{PolyJuMP.ScalarPolynomialFunction{T},S}, ) where {T,S<:MOI.AbstractScalarSet} return haskey(model.constraints, S) && - MOI.is_valid(model.constraints[S][2], ci) + MOI.is_valid(model.constraints[S][2], ci) end function MOI.get( @@ -98,7 +98,7 @@ function MOI.set( end function MOI.set(model::Optimizer, attr::MOI.AbstractModelAttribute, value) - MOI.set(model.model, attr, value) + return MOI.set(model.model, attr, value) end function MOI.get(model::Optimizer, attr::MOI.AbstractModelAttribute) @@ -140,31 +140,22 @@ function MOI.add_constraint( S = typeof(set) if !haskey(model.constraints, S) con = MOI.Utilities.VectorOfConstraints{F,S}() - model.constraints[S] = (P,con) + model.constraints[S] = (P, con) end return MOI.add_constraint(model.constraints[S][2], func, set) end function MOI.get( model::Optimizer{T}, - attr::Union{ - MOI.ConstraintFunction, - MOI.ConstraintSet, - }, - ci::MOI.ConstraintIndex{ - <:PolyJuMP.ScalarPolynomialFunction{T}, - S, - }, + attr::Union{MOI.ConstraintFunction,MOI.ConstraintSet}, + ci::MOI.ConstraintIndex{<:PolyJuMP.ScalarPolynomialFunction{T},S}, ) where {T,S} return MOI.get(model.constraints[S][2], attr, ci) end function MOI.get( model::Optimizer{T}, - attr::MOI.ListOfConstraintIndices{ - <:PolyJuMP.ScalarPolynomialFunction{T}, - S, - } + attr::MOI.ListOfConstraintIndices{<:PolyJuMP.ScalarPolynomialFunction{T},S}, ) where {T,S<:MOI.AbstractScalarSet} return MOI.get(model.constraints[S][2], attr) end @@ -208,13 +199,20 @@ function _add_monomials!(p::PolyJuMP.ScalarPolynomialFunction, monos1) end end -function _subs!(p::PolyJuMP.ScalarPolynomialFunction{T,P}, ::Nothing) where {T,P} - return p, Dict{MOI.VariableIndex,MP.variable_union_type(P)}( +function _subs!( + p::PolyJuMP.ScalarPolynomialFunction{T,P}, + ::Nothing, +) where {T,P} + return p, + Dict{MOI.VariableIndex,MP.variable_union_type(P)}( vi => var for (vi, var) in zip(p.variables, MP.variables(p.polynomial)) ) end -function _subs!(p::PolyJuMP.ScalarPolynomialFunction, index_to_var::Dict{K,V}) where {K,V} +function _subs!( + p::PolyJuMP.ScalarPolynomialFunction, + index_to_var::Dict{K,V}, +) where {K,V} old_var = V[] new_var = V[] for (vi, var) in zip(p.variables, MP.variables(p.polynomial)) @@ -234,7 +232,10 @@ function _subs!(p::PolyJuMP.ScalarPolynomialFunction, index_to_var::Dict{K,V}) w return p, index_to_var end -function _add_variables!(p::PolyJuMP.ScalarPolynomialFunction{T,P}, d) where {T,P} +function _add_variables!( + p::PolyJuMP.ScalarPolynomialFunction{T,P}, + d, +) where {T,P} if isnothing(d) d = Dict{MP.monomial_type(P),MOI.VariableIndex}() else @@ -249,7 +250,12 @@ function _add_variables!(p::PolyJuMP.ScalarPolynomialFunction{T,P}, d) where {T, return d end -function monomial_variable_index(model::Optimizer{T}, d::Dict, div, mono::MP.AbstractMonomialLike) where {T} +function monomial_variable_index( + model::Optimizer{T}, + d::Dict, + div, + mono::MP.AbstractMonomialLike, +) where {T} if !haskey(d, mono) x = div[mono] vx = monomial_variable_index(model, d, div, x) @@ -263,8 +269,10 @@ function monomial_variable_index(model::Optimizer{T}, d::Dict, div, mono::MP.Abs l = max(l, zero(T)) end u = max(bounds...) - d[mono], _ = MOI.add_constrained_variable(model.model, MOI.Interval(l, u)) - MOI.add_constraint(model, + d[mono], _ = + MOI.add_constrained_variable(model.model, MOI.Interval(l, u)) + MOI.add_constraint( + model, MA.@rewrite(one(T) * d[mono] - one(T) * vx * vy), MOI.EqualTo(zero(T)), ) @@ -293,10 +301,16 @@ function MOI.Utilities.final_touch(model::Optimizer{T}, _) where {T} end if !isempty(model.constraints) for S in keys(model.constraints) - for ci in MOI.get(model, MOI.ListOfConstraintIndices{ - PolyJuMP.ScalarPolynomialFunction{T,model.constraints[S][1]}, - S - }()) + for ci in MOI.get( + model, + MOI.ListOfConstraintIndices{ + PolyJuMP.ScalarPolynomialFunction{ + T, + model.constraints[S][1], + }, + S, + }(), + ) func = MOI.get(model, MOI.ConstraintFunction(), ci) func, index_to_var = _subs!(func, index_to_var) vars = _add_variables!(func, vars) diff --git a/src/functions.jl b/src/functions.jl index 975191a..f0d2a2a 100644 --- a/src/functions.jl +++ b/src/functions.jl @@ -75,7 +75,14 @@ function _to_polynomial!( end for t in func.quadratic_terms coef = t.variable_1 == t.variable_2 ? t.coefficient / 2 : t.coefficient - push!(terms, MP.term(coef, _to_polynomial!(d, T, t.variable_1) * _to_polynomial!(d, T, t.variable_2))) + push!( + terms, + MP.term( + coef, + _to_polynomial!(d, T, t.variable_1) * + _to_polynomial!(d, T, t.variable_2), + ), + ) end return MP.polynomial(terms) end @@ -134,10 +141,12 @@ function MOI.Utilities.is_coefficient_type( end # Placeholder for `promote_operation` -struct VectorPolynomialFunction{T,P<:MP.AbstractPolynomial{T}} <: MOI.AbstractVectorFunction -end +struct VectorPolynomialFunction{T,P<:MP.AbstractPolynomial{T}} <: + MOI.AbstractVectorFunction end -function MOI.Utilities.scalar_type(::Type{VectorPolynomialFunction{T,P}}) where {T,P} +function MOI.Utilities.scalar_type( + ::Type{VectorPolynomialFunction{T,P}}, +) where {T,P} return PolyJuMP.ScalarPolynomialFunction{T,P} end @@ -158,7 +167,9 @@ end function MOI.Utilities.promote_operation( ::typeof(-), ::Type{T}, - F::Type{<:Union{ScalarPolynomialFunction{T,P},VectorPolynomialFunction{T,P}}}, + F::Type{ + <:Union{ScalarPolynomialFunction{T,P},VectorPolynomialFunction{T,P}}, + }, ) where {T,P} return F end @@ -195,7 +206,9 @@ function MOI.Utilities.operate( p::ScalarPolynomialFunction{T,P}, f::Union{T,MOI.AbstractScalarFunction}, ) where {T,P} - d = Dict(vi => v for (vi, v) in zip(p.variables, MP.variables(p.polynomial))) + d = Dict( + vi => v for (vi, v) in zip(p.variables, MP.variables(p.polynomial)) + ) poly = _to_polynomial!(d, T, f) return _scalar_polynomial(d, T, op(p.polynomial, poly)) end diff --git a/test/qcqp.jl b/test/qcqp.jl index a3afbbd..427111c 100644 --- a/test/qcqp.jl +++ b/test/qcqp.jl @@ -49,14 +49,16 @@ end function _test_objective_or_constraint(x, y, T, obj::Bool) inner = Model{T}() optimizer = MOI.Utilities.MockOptimizer(inner) - model = PolyJuMP.JuMP.GenericModel{T}(() -> PolyJuMP.QCQP.Optimizer{T}(optimizer)) + model = PolyJuMP.JuMP.GenericModel{T}( + () -> PolyJuMP.QCQP.Optimizer{T}(optimizer), + ) PolyJuMP.@variable(model, 1 <= a <= 2) PolyJuMP.@variable(model, -5 <= b <= 3) PolyJuMP.@constraint(model, a + b >= 1) if obj - PolyJuMP.@objective(model, Min, a^3 - a^2 + 2a*b - b^2 + b^3) + PolyJuMP.@objective(model, Min, a^3 - a^2 + 2a * b - b^2 + b^3) else - PolyJuMP.@constraint(model, 0 <= a^3 - a^2 + 2a*b - b^2 + b^3 <= 1) + PolyJuMP.@constraint(model, 0 <= a^3 - a^2 + 2a * b - b^2 + b^3 <= 1) end PolyJuMP.optimize!(model) vis = MOI.get(inner, MOI.ListOfVariableIndices()) @@ -98,7 +100,9 @@ end function test_objective_and_constraint(x, y, T) inner = Model{T}() optimizer = MOI.Utilities.MockOptimizer(inner) - model = PolyJuMP.JuMP.GenericModel{T}(() -> PolyJuMP.QCQP.Optimizer{T}(optimizer)) + model = PolyJuMP.JuMP.GenericModel{T}( + () -> PolyJuMP.QCQP.Optimizer{T}(optimizer), + ) PolyJuMP.@variable(model, -2 <= a <= 3) PolyJuMP.@variable(model, 5 <= b <= 7) PolyJuMP.@constraint(model, 0 <= a^3 <= 1) @@ -136,7 +140,8 @@ function test_objective_and_constraint(x, y, T) @test length(cis) == 3 @test MOI.get(inner, MOI.ConstraintFunction(), cis[1]) ≈ o * a3 + z * a * b @test MOI.get(inner, MOI.ConstraintFunction(), cis[2]) ≈ o * b3 + z * a * b - @test MOI.get(inner, MOI.ConstraintFunction(), cis[3]) ≈ o * a3 * b3 + o * a6 + @test MOI.get(inner, MOI.ConstraintFunction(), cis[3]) ≈ + o * a3 * b3 + o * a6 end function runtests(x, y)