diff --git a/src/QCQP/MOI_wrapper.jl b/src/QCQP/MOI_wrapper.jl index 54bdbbc..4651972 100644 --- a/src/QCQP/MOI_wrapper.jl +++ b/src/QCQP/MOI_wrapper.jl @@ -3,11 +3,31 @@ import MathOptInterface as MOI const _Model{F,S} = MOI.Utilities.UniversalFallback{MOI.Utilities.VectorOfConstraints{F,S}} +""" + QuadratizationMethod() <: MOI.AbstractOptimizerAttribute + +Quadratization method used by the optimizer, an +[`AbstractQuadratizationMethod`](@ref); [`CompositeQuadratization()`](@ref) +by default. +""" +struct QuadratizationMethod <: MOI.AbstractOptimizerAttribute end + +""" + EnforceConsistencyConstraints() <: MOI.AbstractOptimizerAttribute + +If `true`, adds the redundant constraints `v1 * v2 == v3 * v4` for every +pair of products of variables representing the same monomial. +`false` by default since it can create many constraints. +""" +struct EnforceConsistencyConstraints <: MOI.AbstractOptimizerAttribute end + mutable struct Optimizer{T,O<:MOI.ModelLike} <: MOI.AbstractOptimizer model::O objective::Union{Nothing,PolyJuMP.ScalarPolynomialFunction{T}} constraints::DataStructures.OrderedDict{Type,Tuple{Type,_Model}} index_map::Dict{MOI.ConstraintIndex,MOI.ConstraintIndex} + quadratization_method::AbstractQuadratizationMethod + consistency_constraints::Bool end function Optimizer{T}(optimizer::MOI.ModelLike) where {T} @@ -23,6 +43,8 @@ function Optimizer{T}(optimizer::MOI.ModelLike) where {T} nothing, DataStructures.OrderedDict{Type,MOI.Utilities.VectorOfConstraints}(), Dict{MOI.ConstraintIndex,MOI.ConstraintIndex}(), + CompositeQuadratization(), + false ) end @@ -177,6 +199,21 @@ function MOI.supports( return true end +function MOI.supports( + ::Optimizer, + ::MOI.ObjectiveFunction{MOI.ScalarNonlinearFunction}, +) + return false +end + +function MOI.supports_constraint( + ::Optimizer, + ::Type{MOI.ScalarNonlinearFunction}, + ::Type{<:MOI.AbstractScalarSet}, +) + return false +end + function MOI.set( model::Optimizer{T}, ::MOI.ObjectiveFunction{F}, @@ -194,6 +231,25 @@ function MOI.get(model::Optimizer, attr::MOI.AbstractModelAttribute) return MOI.get(model.model, attr) end +function MOI.supports(model::Optimizer, attr::MOI.AbstractOptimizerAttribute) + return MOI.supports(model.model, attr) +end + +function MOI.supports( + ::Optimizer, + ::Union{QuadratizationMethod,EnforceConsistencyConstraints}, +) + return true +end + +function MOI.set(model::Optimizer, attr::MOI.AbstractOptimizerAttribute, value) + return MOI.set(model.model, attr, value) +end + +function MOI.get(model::Optimizer, attr::MOI.AbstractOptimizerAttribute) + return MOI.get(model.model, attr) +end + function MOI.supports_constraint( model::Optimizer, ::Type{F}, @@ -273,6 +329,39 @@ end MOI.optimize!(model::Optimizer) = MOI.optimize!(model.model) + +MOI.Utilities.map_indices(::Any, ::QuadratizationMethod, x) = x + +function MOI.set( + model::Optimizer{T}, + ::QuadratizationMethod, + value::AbstractQuadratizationMethod +) where {T} + model.quadratization_method = value +end + +function MOI.get( + model::Optimizer{T}, + ::QuadratizationMethod, +) where {T} + return model.quadratization_method +end + +function MOI.set( + model::Optimizer{T}, + ::EnforceConsistencyConstraints, + value::Bool +) where {T} + model.consistency_constraints = value +end + +function MOI.get( + model::Optimizer{T}, + ::EnforceConsistencyConstraints, +) where {T} + return model.consistency_constraints +end + function _quad_convert(p::MP.AbstractPolynomialLike{T}, index, div) where {T} q = zero(MOI.ScalarQuadraticFunction{T}) for t in MP.terms(p) @@ -468,6 +557,64 @@ function _add_constraints( return end +# Keys of the products `v1 * v2` and of the pairs of products `v1*v2, v3*v4` +# already considered, independent of the order of the factors. +_product_key(a::MOI.VariableIndex, b::MOI.VariableIndex) = minmax(a.value, b.value) +_products_key(p, q) = minmax(p, q) + +function _add_consistency_constraints!( + model::Optimizer{T}, + vars::Dict, + monos::Vector{Vector}, + cache::Set, + mono::MP.AbstractMonomialLike; + verbose::Bool, +) where {T} + max_deg = length(monos) + v = vars[mono] + for a in keys(vars) + if MP.degree(a) <= 1 + continue + end + va = vars[a] + key = _product_key(va, v) + if key in cache + continue + end + push!(cache, key) + + deg = MP.degree(a)+MP.degree(mono) + for b in reduce(vcat, monos[max(1, deg-max_deg):min(deg-1,max_deg)]) + if b == a + continue + end + vb = vars[b] + for c in monos[deg-MP.degree(b)] + if c == a + continue + end + vc = vars[c] + key = _products_key(_product_key(va, v), _product_key(vb, vc)) + if key in cache + continue + end + push!(cache, key) + if mono*a == b*c + if verbose + println("Consistency constraint: $mono * $a == $b * $c") + end + MOI.Utilities.normalize_and_add_constraint( + model, + MA.@rewrite(one(T)*v*va - one(T)*vb*vc), + MOI.EqualTo(zero(T)); + allow_modify_function = true, + ) + end + end + end + end +end + function MOI.Utilities.final_touch(model::Optimizer{T}, _) where {T} index_to_var = nothing vars = nothing @@ -487,7 +634,9 @@ function MOI.Utilities.final_touch(model::Optimizer{T}, _) where {T} end end if !isnothing(monos) - div = decompose(monos) + max_deg = maximum(MP.degree.(monos)) + verbose = MOI.supports(model, MOI.Silent()) && !MOI.get(model, MOI.Silent()) + div = decompose(monos, model.quadratization_method, verbose=verbose) for mono in sort(collect(keys(div))) if haskey(vars, mono) continue @@ -497,6 +646,28 @@ function MOI.Utilities.final_touch(model::Optimizer{T}, _) where {T} b = MP.div_multiple(mono, a) monomial_variable_index(model, vars, div, b) end + if verbose + println("Done with quadratization!") + end + if model.consistency_constraints + cache = Set{Any}() + sorted_monos = Vector[[] for _ in 1:max_deg] + for mono in keys(vars) + deg = MP.degree(mono) + push!(sorted_monos[deg], mono) + end + i = 1 + l = length(keys(vars)) + for mono in keys(vars) + if MP.degree(mono) > 1 + _add_consistency_constraints!(model, vars, sorted_monos, cache, mono; verbose) + end + if verbose + println("Consistency constraints: finished monomial $i/$l") + end + i += 1 + end + end if !isnothing(model.objective) func, index_to_var = _subs!(model.objective, index_to_var) obj = _quad_convert(func.polynomial, vars, div) diff --git a/src/QCQP/QCQP.jl b/src/QCQP/QCQP.jl index 47fbc67..27895e2 100644 --- a/src/QCQP/QCQP.jl +++ b/src/QCQP/QCQP.jl @@ -2,11 +2,99 @@ module QCQP import MutableArithmetics as MA import MultivariatePolynomials as MP +import DynamicPolynomials as DP import DataStructures +import MathOptInterface as MOI import PolyJuMP +""" + abstract type AbstractQuadratizationMethod end + +Method choosing how each monomial of degree larger than 2 is written as the +product of two monomials of lower degree, see [`decompose`](@ref) and the +[`QuadratizationMethod`](@ref) attribute. +""" +abstract type AbstractQuadratizationMethod end + +""" + GreedyQuadratization() + +Monomials are decomposed by repeatedly adding the most frequent missing factor. +""" +struct GreedyQuadratization <: AbstractQuadratizationMethod end + +""" + @enum Decomposition Greedy Balanced Sequential + +Third phase of [`CompositeQuadratization`](@ref): + +* `Greedy`: substitutes the most frequent pair of variables first; +* `Balanced`: splits multilinear monomials into two halves, requires the smart decomposition; +* `Sequential`: sequential decomposition, separating the first two variables from the rest. +""" +@enum Decomposition begin + Greedy + Balanced + Sequential +end + +""" + CompositeQuadratization(; + decompo::Decomposition = Balanced, + smart::Bool = true, + variable_reuse::Bool = true, + ) + +Multi-level quadratization: if `smart`, the smart decomposition first exploits +even exponents and leaves only multilinear monomials, which are then decomposed +by `decompo`. With `variable_reuse`, a monomial that is the product of two +monomials already known is decomposed as such. +""" +struct CompositeQuadratization <: AbstractQuadratizationMethod + decomposition::Decomposition + smart::Bool + variable_reuse::Bool + function CompositeQuadratization(; + decompo::Decomposition = Balanced, + smart::Bool = true, + variable_reuse::Bool = true, + ) + if decompo == Balanced && !smart + throw(ArgumentError("The balanced decomposition requires `smart = true`.")) + end + return new(decompo, smart, variable_reuse) + end +end + +""" + OptimalQuadratization(optimizer::MOI.AbstractOptimizer) + +Minimizes the number of auxiliary variables by solving a binary program with +`optimizer`. The candidate auxiliary variables are the products of pairs of +variables and products of these candidates, so the result is minimal among these candidates. +""" +struct OptimalQuadratization <: AbstractQuadratizationMethod + optimizer::MOI.AbstractOptimizer +end + + +# Maps the variables to themselves and decomposes the monomials of degree 2 +# of `monos` and of the decompositions of `quad` into their variables. +function _decompose_quadratic!(quad, vars, monos) + for x in vars + quad[x] = x + for y in vars + quad[x*y] = x + end + end + return quad +end + # TODO special case for squares ? -function decompose(monos::AbstractVector{M}) where {M<:MP.AbstractMonomial} +function decompose(monos::AbstractVector{M}, ::GreedyQuadratization; verbose::Bool=true) where {M<:MP.AbstractMonomial} + if verbose + println("Running greedy quadratization scheme...") + end vars = MP.variables(monos) quad = DataStructures.OrderedDict{M,Union{Nothing,M}}( var => var for var in vars @@ -43,6 +131,525 @@ function decompose(monos::AbstractVector{M}) where {M<:MP.AbstractMonomial} return quad end + +""" + kac_phase(monos; variable_reuse::Bool, multilinear::Bool) + +Decomposition inspired by Karia, Adjiman and Chachuat (2022), shared by the +smart decomposition (`multilinear = false`) and the +balanced decomposition (`multilinear = true`): monomials with an +odd exponent are split into their odd part and their even part, monomials +with even exponents into the square of their half and, if `multilinear`, +multilinear monomials into two halves; otherwise multilinear monomials are +left undecomposed (`nothing`) for the next phase. +""" +function kac_phase(monos::AbstractVector{M}; variable_reuse::Bool, multilinear::Bool, verbose::Bool=true) where {M<:MP.AbstractMonomial} + # This function exists because smart decomposition into balanced decomposition + # used to be one method, so the code evolved into having this function, and when + # the method was split this function was kept but parametrized to avoid breaking anything, + # it probably could do with a refactor at some point + vars = MP.variables(monos) + quad = DataStructures.OrderedDict{M,Union{Nothing,M}}() + + stack = DataStructures.Stack{M}() + for mono in monos + if MP.degree(mono) == 0 + continue + end + quad[mono] = nothing + push!(stack, mono) + end + + while !isempty(stack) + mono = pop!(stack) + deg = MP.degree(mono) + if deg <= 2 || get(quad, mono, nothing) !== nothing + continue + end + + if variable_reuse + stop = false + for a in keys(quad) + if a != mono && MP.divides(a, mono) + b = MP.div_multiple(mono, a) + if haskey(quad, b) + quad[mono] = a + stop = true + break + end + end + end + if stop + continue + end + end + + l, r = nothing, nothing + mono_vars = MP.variables(mono) + exps = MP.exponents(mono) + if maximum(exps) == 1 && multilinear + v = filter(x -> MP.degree(mono, x) == 1, mono_vars) + l = prod(v[1:deg÷2]) + r = prod(v[deg÷2+1:deg]) + elseif maximum(exps) == 1 && !multilinear + quad[mono] = nothing + continue + elseif any((exps .% 2) .== 1) + l = MP.div_multiple(mono, mono) + r = copy(mono) + for (v, e) in MP.powers(mono) + if e % 2 == 1 + l *= v + r = MP.div_multiple(r, v) + end + end + else + l = MP.div_multiple(mono, mono) + for (v, e) in MP.powers(mono) + l *= v^(e÷2) + end + end + + if l !== nothing + quad[mono] = l + push!(stack, l) + end + if r !== nothing + push!(stack, r) + end + end + + return quad +end + + +function smart_decompo_phase(monos::AbstractVector{M}; variable_reuse::Bool, verbose::Bool=true) where {M<:MP.AbstractMonomial} + return kac_phase(monos, variable_reuse=variable_reuse, multilinear=false, verbose=verbose) +end + + +function cut_in_half_phase(monos::AbstractVector{M}; variable_reuse::Bool, verbose::Bool=true) where {M<:MP.AbstractMonomial} + return kac_phase(monos, variable_reuse=variable_reuse, multilinear=true, verbose=verbose) +end + + +function greedy_phase(monos::AbstractVector{M}; verbose::Bool=true) where {M<:MP.AbstractMonomial} + vars = MP.variables(monos) + quad = DataStructures.OrderedDict{M,Union{Nothing,M}}( + var => var for var in vars + ) + + var_to_mono = Dict{M, M}(var => var for var in vars) + mono_to_var = Dict{M, M}() + + occurrence = DataStructures.PriorityQueue{M, Tuple{Int, M}}(Base.Order.Reverse) + + monos_list = M[mono for mono in monos] + while any(x -> MP.degree(x) > 2, monos_list) + empty!(occurrence) + for mono in monos_list + if MP.degree(mono) <= 2 + continue + end + vs = MP.variables(mono) + exps = MP.exponents(mono) + nonzeroidx = [i for i in eachindex(exps) if exps[i] > 0] + for k in eachindex(nonzeroidx) + i = nonzeroidx[k] + for j in nonzeroidx[k:end] + if i == j && exps[i] <= 1 + continue + end + q = vs[i]*vs[j] + curr, _ = get(occurrence, q, (0, q)) + occ = i == j ? exps[i]÷2 : min(exps[i], exps[j]) + occurrence[q] = (curr + occ, q) + end + end + end + + for (i, mono) in enumerate(monos_list) + if MP.degree(mono) <= 2 + continue + end + for a in keys(occurrence) + if MP.divides(a, mono) + deg = 0 + b = mono + while MP.divides(a, b) + deg += 1 + b = MP.div_multiple(b, a) + end + if haskey(mono_to_var, a) + v = mono_to_var[a] + else + v = MP.similar_variable(vars[begin], :v) + var_to_mono[v] = a + mono_to_var[a] = v + end + monos_list[i] = b*v^deg + break + end + end + end + end + + _evaluate(x) = begin + mono = var_to_mono[x] + vars = MP.variables(mono) + if MP.degree(mono) == 1 + return x + else + nonzero_vars = filter(v -> MP.degree(mono, v) >= 1, vars) + if length(nonzero_vars) == 1 + nonzero_vars = [nonzero_vars[1], nonzero_vars[1]] + end + return prod(_evaluate.(nonzero_vars), init=one(M)) + end + end + + stack = DataStructures.Stack{M}() + for mono in monos_list + if MP.degree(mono) > 1 + push!(stack, mono) + end + end + + while !isempty(stack) + bilin = pop!(stack) + vars = filter(v -> MP.degree(bilin, v) >= 1, MP.variables(bilin)) + a, b = length(vars) == 2 ? (vars[1], vars[2]) : (vars[1], vars[1]) + a_mono, b_mono = _evaluate(a), _evaluate(b) + mono = a_mono * b_mono + + if haskey(quad, mono) && !isnothing(quad[mono]) + continue + end + + quad[mono] = a_mono + + if !haskey(quad, a_mono) && MP.degree(a_mono) > 2 + quad[a_mono] = nothing + push!(stack, var_to_mono[a]) + end + if !haskey(quad, b_mono) && MP.degree(b_mono) > 2 + quad[b_mono] = nothing + push!(stack, var_to_mono[b]) + end + end + + return quad +end + + +function fallback_phase(monos::AbstractVector{M}; verbose::Bool=true) 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 + deg = MP.degree(mono) + if deg <= 2 + continue + end + + terms = M[] + for (var, n) in MP.powers(mono) + for _ in 1:n + push!(terms, var) + end + end + + a = (deg-1)÷2 + epp = prod(terms[2*a+1:deg]) + quad[mono] = epp + + terms2 = terms[1:2*a] + ep = prod(terms2) + while MP.degree(ep) > 2 + if haskey(quad, ep) + break + end + div = prod(terms2[1:2]) + quad[ep] = div + terms2 = terms2[3:end] + ep = prod(terms2) + end + end + + return quad +end + + +function decompose(monos::AbstractVector{M}, method::OptimalQuadratization; verbose::Bool=true) where {M<:MP.AbstractMonomial} + if verbose + println("Running minimum auxiliary variables quadratization scheme...") + end + + vars = MP.variables(monos) + quad = DataStructures.OrderedDict{M,Union{Nothing,M}}( + var => var for var in vars + ) + deg = maximum(MP.degree.(monos)) + + var_to_mono = DataStructures.OrderedDict{MOI.VariableIndex, M}() + var_divisor = DataStructures.OrderedDict{MOI.VariableIndex, Union{M,MOI.VariableIndex}}() + + optimizer = method.optimizer + all_monos = monos + monos = filter(x -> MP.degree(x) >= 3, monos) + tmp = Set{M}() + # find all pairs of variables + for mono in monos + vs = MP.variables(mono) + exps = MP.exponents(mono) + nonzeroidx = [i for i in eachindex(exps) if exps[i] > 0] + for k in eachindex(nonzeroidx) + i = nonzeroidx[k] + for j in nonzeroidx[k:end] + bilin = vs[i]*vs[j] + if (i == j && exps[i] <= 1) || bilin in tmp + continue + end + new_var, _ = MOI.add_constrained_variable(optimizer, MOI.ZeroOne()) + var_to_mono[new_var] = bilin + var_divisor[new_var] = MP.div_multiple(bilin, vs[j]) + push!(tmp, bilin) + end + end + end + empty!(tmp) + + if verbose + println("found all pairs...") + end + + explored_pairs = Set{Tuple{MOI.VariableIndex,Union{MOI.VariableIndex,eltype(vars)}}}() + curr_deg = 2 + while curr_deg < deg-1 + aux_vars = collect(keys(var_to_mono)) + for (i, va) in enumerate(aux_vars) + a = var_to_mono[va] + for vb in aux_vars[i:end] + b = var_to_mono[vb] + if (va, vb) in explored_pairs || (vb, va) in explored_pairs + continue + end + push!(explored_pairs, (va, vb)) + stop = true + for mono in monos + if MP.divides(a*b, mono) + stop = false + break + end + end + if stop + continue + end + + if MP.degree(a*b) >= deg + continue + end + new_var, _ = MOI.add_constrained_variable(optimizer, MOI.ZeroOne()) + + x = MOI.ScalarAffineTerm(1.0, new_var) + y = MOI.ScalarAffineTerm(-1.0, va) + MOI.add_constraint(optimizer, MOI.ScalarAffineFunction([x, y], 0.0), MOI.LessThan(0.0)) + + x = MOI.ScalarAffineTerm(1.0, new_var) + y = MOI.ScalarAffineTerm(-1.0, vb) + MOI.add_constraint(optimizer, MOI.ScalarAffineFunction([x, y], 0.0), MOI.LessThan(0.0)) + + var_to_mono[new_var] = a*b + var_divisor[new_var] = va + end + for v in vars + if (va, v) in explored_pairs + continue + end + push!(explored_pairs, (va, v)) + stop = true + for mono in monos + if MP.divides(a*v, mono) + stop = false + break + end + end + if stop + continue + end + + new_var, _ = MOI.add_constrained_variable(optimizer, MOI.ZeroOne()) + + x = MOI.ScalarAffineTerm(1.0, new_var) + y = MOI.ScalarAffineTerm(-1.0, va) + MOI.add_constraint(optimizer, MOI.ScalarAffineFunction([x, y], 0.0), MOI.LessThan(0.0)) + + var_to_mono[new_var] = a*v + var_divisor[new_var] = va + end + end + if verbose + println("iteration deg = $curr_deg out of $deg done...") + end + curr_deg *= 2 + end + + if verbose + println("done with variable creation...") + end + + mono_to_vars = Dict{M, Vector{MOI.VariableIndex}}() + for (v, m) in var_to_mono + if !haskey(mono_to_vars, m) + mono_to_vars[m] = MOI.VariableIndex[] + end + push!(mono_to_vars[m], v) + end + + for mono in monos + quad_terms = MOI.ScalarQuadraticTerm{Float64}[] + affine_terms = MOI.ScalarAffineTerm{Float64}[] + for (a, vs) in mono_to_vars + if MP.divides(a, mono) + b = MP.div_multiple(mono, a) + if b == a || MP.degree(b) == 1 + for v in vs + push!(affine_terms, MOI.ScalarAffineTerm(1.0, v)) + end + elseif MP.degree(b) > 1 + for va in vs + for vb in mono_to_vars[b] + push!(quad_terms, MOI.ScalarQuadraticTerm(1.0, va, vb)) + end + end + end + end + end + MOI.add_constraint(optimizer, MOI.ScalarQuadraticFunction(quad_terms, affine_terms, 0.0), MOI.GreaterThan(1.0)) + end + + if verbose + println("done with constraint creation...") + end + + affine_terms = MOI.ScalarAffineTerm{Float64}[] + for var in keys(var_to_mono) + push!(affine_terms, MOI.ScalarAffineTerm(1.0, var)) + end + obj = MOI.ScalarAffineFunction(affine_terms, 0.0) + MOI.set(optimizer, MOI.ObjectiveFunction{typeof(obj)}(), obj) + MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MIN_SENSE) + + MOI.optimize!(optimizer) + status = MOI.get(optimizer, MOI.TerminationStatus()) + if status != MOI.OPTIMAL + println("OptimalQuadratization: the binary program selecting the auxiliary variables terminated with status $status.") + end + + count_new_vars = 0 + active_vars = MOI.VariableIndex[] + for (var, mono) in var_to_mono + active = MOI.get(optimizer, MOI.VariablePrimal(), var) + if active ≈ 1 + if verbose + println("$mono is active") + end + if haskey(quad, mono) && verbose + # This should never happen + println("[WARNING] quad already has $mono, which means that $mono was activated twice...") + end + if var_divisor[var] isa M + quad[mono] = var_divisor[var] + else + quad[mono] = var_to_mono[var_divisor[var]] + end + count_new_vars += 1 + push!(active_vars, var) + end + end + + if verbose + println("should generate $count_new_vars new variables") + end + + mono_set = Set{M}(monos) + for (i, va) in enumerate(active_vars) + a = var_to_mono[va] + for vb in active_vars[i:end] + b = var_to_mono[vb] + mono = a*b + if mono in mono_set && !haskey(quad, mono) + quad[mono] = a + end + end + end + for mono in monos + if haskey(quad, mono) + continue + end + for (var, exp) in MP.powers(mono) + if exp > 0 && haskey(quad, MP.div_multiple(mono, var)) + quad[mono] = MP.div_multiple(mono, var) + break + end + end + end + + + _decompose_quadratic!(quad, vars, all_monos) + for mono in keys(quad) + if quad[mono] === nothing + error("OptimalQuadratization: no decomposition was selected for $mono.") + end + end + + MOI.empty!(optimizer) + + return quad +end + + +function decompose(monos::AbstractVector{M}, method::CompositeQuadratization; verbose::Bool=true) where {M<:MP.AbstractMonomial} + if verbose + println("Running composite quadratization scheme...") + if method.smart + print("Smart ") + end + print(lowercase(string(method.decomposition))) + if method.variable_reuse + print(" with variable reuse") + end + print("\n") + end + + vars = MP.variables(monos) + quad = DataStructures.OrderedDict{M,Union{Nothing,M}}() + + if method.smart + quad = smart_decompo_phase(monos, variable_reuse=method.variable_reuse) + monos = filter(x -> isnothing(quad[x]), collect(keys(quad))) + end + + if method.decomposition == Greedy + q = greedy_phase(monos, verbose=verbose) + elseif method.decomposition == Balanced # implies `method.smart` + q = cut_in_half_phase(monos, variable_reuse=method.variable_reuse, verbose=verbose) + else + q = fallback_phase(monos, verbose=verbose) + end + + for (k, v) in q + if isnothing(get(quad, k, nothing)) + quad[k] = v + end + end + + return _decompose_quadratic!(quad, vars, monos) +end + + include("MOI_wrapper.jl") end diff --git a/test/qcqp.jl b/test/qcqp.jl index 5d695a3..a464ca8 100644 --- a/test/qcqp.jl +++ b/test/qcqp.jl @@ -16,25 +16,66 @@ function test_solver_name(_, _, _) @test MOI.get(optimizer, MOI.SolverName()) == "PolyJuMP.QCQP with Mock" end -function _test_decompose(monos, exps) +function _test_decompose(monos, args...) 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 + expected = PolyJuMP.QCQP.DataStructures.OrderedDict{M,M}() + for x in vars + expected[x] = x + for y in vars + expected[x*y] = x + end + end + for (quad, exps) in args + for exp in exps + expected[exp[1]] = exp[2] + end + quad = PolyJuMP.QCQP.decompose(monos, quad) + @test quad == expected + end 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], + monos = [x^2, y^3, x^2 * y^3] + + _test_decompose( + monos, + PolyJuMP.QCQP.CompositeQuadratization(decompo=PolyJuMP.QCQP.Sequential, smart=false, variable_reuse=false) => + [y^3 => y, x^2 * y^3 => y, x^2 * y^2 => x^2], + ) + _test_decompose( + monos, + PolyJuMP.QCQP.CompositeQuadratization(decompo=PolyJuMP.QCQP.Sequential, smart=true, variable_reuse=false) => + [y^3 => y, x^2 * y^3 => y, x^2 * y^2 => x * y], + ) + _test_decompose( + monos, + PolyJuMP.QCQP.CompositeQuadratization(decompo=PolyJuMP.QCQP.Sequential, smart=false, variable_reuse=true) => + [y^3 => y, x^2 * y^3 => y, x^2 * y^2 => x^2], ) + _test_decompose( + monos, + PolyJuMP.QCQP.CompositeQuadratization(decompo=PolyJuMP.QCQP.Sequential, smart=true, variable_reuse=true) => + [y^3 => y, x^2 * y^3 => x^2], + ) + # Without the smart phase, the greedy phase has to create new variables + _test_decompose( + monos, + PolyJuMP.QCQP.CompositeQuadratization(decompo=PolyJuMP.QCQP.Greedy, smart=false, variable_reuse=false) => + [y^3 => y, x^2 * y^3 => x * y, x * y^2 => y], + ) + for decompo in (PolyJuMP.QCQP.Greedy, PolyJuMP.QCQP.Balanced) + _test_decompose( + monos, + PolyJuMP.QCQP.CompositeQuadratization(decompo=decompo, smart=true, variable_reuse=false) => + [y^3 => y, x^2 * y^3 => y, x^2 * y^2 => x * y], + ) + _test_decompose( + monos, + PolyJuMP.QCQP.CompositeQuadratization(decompo=decompo, smart=true, variable_reuse=true) => + [y^3 => y, x^2 * y^3 => x^2], + ) + end end MOI.Utilities.@model( @@ -137,7 +178,7 @@ function test_objective_and_constraint(x, y, T) PolyJuMP.optimize!(model) vis = MOI.get(inner, MOI.ListOfVariableIndices()) @test length(vis) == 7 - a, b, b2, a2, a3, b3, a6 = vis + a, b, b2, a2, b3, a3, 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) @@ -153,8 +194,8 @@ function test_objective_and_constraint(x, y, 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[3]) ≈ b3 - o * b2 * b + @test MOI.get(inner, MOI.ConstraintFunction(), cis[4]) ≈ a3 - o * a2 * a @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))