From 2fde3bb89893e2d7cd1ae80120fc5d97ee02886c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Mon, 2 Oct 2023 09:03:12 +0200 Subject: [PATCH] Polynomial to QCQP meta-solver (#95) * Polynomial to QCQP meta-solver * Fix format * Add MOI wrapper * Fix format * Fixes * Fixes * Add test * Add bounds * Add suppor to constraints * Fix format --- Project.toml | 2 + src/PolyJuMP.jl | 1 + src/QCQP/MOI_wrapper.jl | 347 ++++++++++++++++++++++++++++++++++++++++ src/QCQP/QCQP.jl | 48 ++++++ src/functions.jl | 91 +++++++++-- src/nl_to_polynomial.jl | 9 +- test/qcqp.jl | 170 ++++++++++++++++++++ 7 files changed, 652 insertions(+), 16 deletions(-) create mode 100644 src/QCQP/MOI_wrapper.jl 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/MOI_wrapper.jl b/src/QCQP/MOI_wrapper.jl new file mode 100644 index 0000000..9e6d25d --- /dev/null +++ b/src/QCQP/MOI_wrapper.jl @@ -0,0 +1,347 @@ +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{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, +) 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) +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, + 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( + 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.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) + return 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}, + ::Type{S}, +) 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, + func::MOI.AbstractFunction, + set::MOI.AbstractSet, +) + 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) +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(p::MP.AbstractPolynomialLike{T}, index, 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 + 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]) + end + end + end + return q +end + +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] + 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)), + ) + end + 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) + 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 + end + 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 + +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 new file mode 100644 index 0000000..47fbc67 --- /dev/null +++ b/src/QCQP/QCQP.jl @@ -0,0 +1,48 @@ +module QCQP + +import MutableArithmetics as MA +import MultivariatePolynomials as MP +import DataStructures +import PolyJuMP + +# 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 + ) + for mono in monos + if MP.degree(mono) > 1 + quad[mono] = nothing + end + 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 + +include("MOI_wrapper.jl") + +end diff --git a/src/functions.jl b/src/functions.jl index 7fa3245..f0d2a2a 100644 --- a/src/functions.jl +++ b/src/functions.jl @@ -1,6 +1,6 @@ """ - struct ScalarPolynomialFunction{P<:MP.AbstractPolynomialLike} <: MOI.AbstractScalarFunction - p::P + struct ScalarPolynomialFunction{T,P<:MP.AbstractPolynomial{T}} <: MOI.AbstractScalarFunction + polynomial::P variables::Vector{MOI.VariableIndex} end @@ -64,6 +64,29 @@ 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 +96,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( @@ -124,10 +140,36 @@ function MOI.Utilities.is_coefficient_type( return S === T end +# 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(-), ::Type{T}, - F::Type{ScalarPolynomialFunction{T,P}}, + F::Type{ + <:Union{ScalarPolynomialFunction{T,P},VectorPolynomialFunction{T,P}}, + }, ) where {T,P} return F end @@ -141,11 +183,32 @@ function MOI.Utilities.promote_operation( return F end -# FIXME +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}, ::Type{ScalarPolynomialFunction{T,P}}, ) where {T,P} - return MOI.VectorQuadraticFunction{T} + return VectorPolynomialFunction{T,P} +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..bb799e6 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]) + sort!(variable_map, by = x -> x[2], rev = true) 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} diff --git a/test/qcqp.jl b/test/qcqp.jl new file mode 100644 index 0000000..427111c --- /dev/null +++ b/test/qcqp.jl @@ -0,0 +1,170 @@ +module TestQCQP + +using Test + +import MathOptInterface as MOI +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]) + 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 + +MOI.Utilities.@model( + Model, + (), + (MOI.LessThan, MOI.GreaterThan, MOI.EqualTo, MOI.Interval), + (), + (), + (), + (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction), + (), + (), +) + +function MOI.supports( + ::Model, + ::MOI.ObjectiveFunction{MOI.ScalarNonlinearFunction}, +) + return false +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), + ) + 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) + 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 + 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 + 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) + for name in names(@__MODULE__; all = true) + if startswith("$name", "test_") + @testset "$(name) $T" for T in [Int, Float64] + 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) +end +import TypedPolynomials +@testset "TypedPolynomials" begin + TypedPolynomials.@polyvar(x, y) + TestQCQP.runtests(x, y) +end