From 123f7e0c0d7e0500cf37fece625567d4dbb0a69f Mon Sep 17 00:00:00 2001 From: votroto Date: Sat, 9 Mar 2024 15:27:10 +0100 Subject: [PATCH 01/16] fix variable order (first draft) --- src/QCQP/MOI_wrapper.jl | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/src/QCQP/MOI_wrapper.jl b/src/QCQP/MOI_wrapper.jl index 2b14c7d..d712a4a 100644 --- a/src/QCQP/MOI_wrapper.jl +++ b/src/QCQP/MOI_wrapper.jl @@ -237,6 +237,7 @@ end function _add_variables!( p::PolyJuMP.ScalarPolynomialFunction{T,P}, d, + va, ) where {T,P} if isnothing(d) d = Dict{MP.monomial_type(P),MOI.VariableIndex}() @@ -246,8 +247,8 @@ function _add_variables!( d = convert(Dict{M,MOI.VariableIndex}, d) end end - for (v, vi) in zip(MP.variables(p.polynomial), p.variables) - d[v] = vi + for (v, vi) in va + d[vi] = v end return d end @@ -333,7 +334,7 @@ function MOI.Utilities.final_touch(model::Optimizer{T}, _) where {T} monos = nothing if !isnothing(model.objective) func, index_to_var = _subs!(model.objective, index_to_var) - vars = _add_variables!(func, vars) + vars = _add_variables!(func, vars, index_to_var) monos = _add_monomials!(func, monos) end if !isempty(model.constraints) @@ -350,7 +351,7 @@ function MOI.Utilities.final_touch(model::Optimizer{T}, _) where {T} ) func = MOI.get(model, MOI.ConstraintFunction(), ci) func, index_to_var = _subs!(func, index_to_var) - vars = _add_variables!(func, vars) + vars = _add_variables!(func, vars, index_to_var) monos = _add_monomials!(func, monos) end end From 71bbe0046b5de22b29060a3765ebfc96d13c4fc8 Mon Sep 17 00:00:00 2001 From: votroto Date: Mon, 11 Mar 2024 15:26:25 +0100 Subject: [PATCH 02/16] fix qcqp tests --- test/qcqp.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/qcqp.jl b/test/qcqp.jl index 6f48007..fce6d3d 100644 --- a/test/qcqp.jl +++ b/test/qcqp.jl @@ -149,8 +149,8 @@ function test_no_monomials(x, y, T) model = PolyJuMP.JuMP.GenericModel{T}() do return PolyJuMP.QCQP.Optimizer{T}(MOI.Utilities.MockOptimizer(inner)) end - PolyJuMP.@variable(model, 0 <= x[1:2] <= 1) - PolyJuMP.@constraint(model, x[1] * x[2] == 0.5) + PolyJuMP.@variable(model, 0 <= x[1:2] <= 2) + PolyJuMP.@constraint(model, x[1] * x[2] == 1) PolyJuMP.@objective(model, Min, sum(x)) PolyJuMP.optimize!(model) @test MOI.get(inner, MOI.NumberOfVariables()) == 2 @@ -209,7 +209,7 @@ function test_scalar_nonlinear_function(x, y, T) PolyJuMP.@expression(model, g, x^2) PolyJuMP.@constraint(model, f * g == 0) PolyJuMP.optimize!(model) - F, S = ScalarQuadraticFunction{T}, EqualTo{T} + F, S = MOI.ScalarQuadraticFunction{T}, MOI.EqualTo{T} @test MOI.get(inner, MOI.NumberOfConstraints{F,S}()) == 2 @test MOI.get(inner, MOI.NumberOfVariables()) == 2 return From 7a264f0ba0c369d872f186d542c54dfde63958a5 Mon Sep 17 00:00:00 2001 From: votroto Date: Mon, 11 Mar 2024 16:32:44 +0100 Subject: [PATCH 03/16] remove redundant function --- src/functions.jl | 23 ----------------------- 1 file changed, 23 deletions(-) diff --git a/src/functions.jl b/src/functions.jl index f0d2a2a..e52f148 100644 --- a/src/functions.jl +++ b/src/functions.jl @@ -64,29 +64,6 @@ 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}, From 6cfd7ea013febe17ac230f262c2ed8d00eb029ed Mon Sep 17 00:00:00 2001 From: votroto Date: Tue, 12 Mar 2024 09:54:55 +0100 Subject: [PATCH 04/16] avoid rational poly if possible --- src/nl_to_polynomial.jl | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/nl_to_polynomial.jl b/src/nl_to_polynomial.jl index 71f0a21..14776ae 100644 --- a/src/nl_to_polynomial.jl +++ b/src/nl_to_polynomial.jl @@ -94,7 +94,11 @@ function _to_polynomial!( return a^b elseif _is_operator(expr, :/) && length(operands) == 2 a, b = _to_polynomial!.(Ref(d), T, operands) - return a / b + if iszero(MP.rem(a, b)) + return MP.div(a, b) + else + return a / b + end elseif _is_variable(expr) return _to_polynomial!(d, T, operands[1]) else From 06602c2ee3ce56689ebe5e664eb05c86864c4c0b Mon Sep 17 00:00:00 2001 From: votroto Date: Tue, 12 Mar 2024 09:59:12 +0100 Subject: [PATCH 05/16] resort moi by updated vars --- src/QCQP/MOI_wrapper.jl | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) diff --git a/src/QCQP/MOI_wrapper.jl b/src/QCQP/MOI_wrapper.jl index d712a4a..f6ab698 100644 --- a/src/QCQP/MOI_wrapper.jl +++ b/src/QCQP/MOI_wrapper.jl @@ -228,8 +228,16 @@ function _subs!( end end if !isempty(old_var) + old_var_map = Dict(zip(old_var, new_var)) + all_old_vars = MP.variables(p.polynomial) + + updated_vars = [get(old_var_map, v, v) for v in all_old_vars] + variable_map = collect(zip(p.variables, updated_vars)) + sort!(variable_map, by = x -> x[2], rev = true) + moi_variables = [x[1] for x in variable_map] + poly = MP.subs(p.polynomial, old_var => new_var) - p = PolyJuMP.ScalarPolynomialFunction(poly, p.variables) + p = PolyJuMP.ScalarPolynomialFunction(poly, moi_variables) end return p, index_to_var end @@ -247,8 +255,8 @@ function _add_variables!( d = convert(Dict{M,MOI.VariableIndex}, d) end end - for (v, vi) in va - d[vi] = v + for (v, vi) in zip(MP.variables(p.polynomial), p.variables) + d[v] = vi end return d end From c0886b373cd2b29c6a263ff21c2b3dc84e226685 Mon Sep 17 00:00:00 2001 From: votroto Date: Tue, 12 Mar 2024 10:13:50 +0100 Subject: [PATCH 06/16] remove redundant arg --- src/QCQP/MOI_wrapper.jl | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/QCQP/MOI_wrapper.jl b/src/QCQP/MOI_wrapper.jl index f6ab698..0910510 100644 --- a/src/QCQP/MOI_wrapper.jl +++ b/src/QCQP/MOI_wrapper.jl @@ -245,7 +245,6 @@ end function _add_variables!( p::PolyJuMP.ScalarPolynomialFunction{T,P}, d, - va, ) where {T,P} if isnothing(d) d = Dict{MP.monomial_type(P),MOI.VariableIndex}() @@ -342,7 +341,7 @@ function MOI.Utilities.final_touch(model::Optimizer{T}, _) where {T} monos = nothing if !isnothing(model.objective) func, index_to_var = _subs!(model.objective, index_to_var) - vars = _add_variables!(func, vars, index_to_var) + vars = _add_variables!(func, vars) monos = _add_monomials!(func, monos) end if !isempty(model.constraints) @@ -359,7 +358,7 @@ function MOI.Utilities.final_touch(model::Optimizer{T}, _) where {T} ) func = MOI.get(model, MOI.ConstraintFunction(), ci) func, index_to_var = _subs!(func, index_to_var) - vars = _add_variables!(func, vars, index_to_var) + vars = _add_variables!(func, vars) monos = _add_monomials!(func, monos) end end From 687fdbc249647b521bf40ffac2ccad3c22311d8f Mon Sep 17 00:00:00 2001 From: votroto Date: Tue, 12 Mar 2024 14:56:24 +0100 Subject: [PATCH 07/16] correct moi order with dicts --- src/QCQP/MOI_wrapper.jl | 14 ++++++-------- src/nl_to_polynomial.jl | 7 +++---- 2 files changed, 9 insertions(+), 12 deletions(-) diff --git a/src/QCQP/MOI_wrapper.jl b/src/QCQP/MOI_wrapper.jl index 0910510..42fadd5 100644 --- a/src/QCQP/MOI_wrapper.jl +++ b/src/QCQP/MOI_wrapper.jl @@ -228,16 +228,14 @@ function _subs!( end end if !isempty(old_var) - old_var_map = Dict(zip(old_var, new_var)) - all_old_vars = MP.variables(p.polynomial) - - updated_vars = [get(old_var_map, v, v) for v in all_old_vars] - variable_map = collect(zip(p.variables, updated_vars)) - sort!(variable_map, by = x -> x[2], rev = true) - moi_variables = [x[1] for x in variable_map] + to_old_map = Dict(zip(new_var, old_var)) + to_moi_map = Dict(zip(MP.variables(p.polynomial), p.variables)) poly = MP.subs(p.polynomial, old_var => new_var) - p = PolyJuMP.ScalarPolynomialFunction(poly, moi_variables) + all_new_vars = MP.variables(poly) + moi_vars = [to_moi_map[get(to_old_map, v, v)] for v in all_new_vars] + + p = PolyJuMP.ScalarPolynomialFunction(poly, moi_vars) end return p, index_to_var end diff --git a/src/nl_to_polynomial.jl b/src/nl_to_polynomial.jl index 14776ae..b67cbad 100644 --- a/src/nl_to_polynomial.jl +++ b/src/nl_to_polynomial.jl @@ -98,7 +98,7 @@ function _to_polynomial!( return MP.div(a, b) else return a / b - end + end elseif _is_variable(expr) return _to_polynomial!(d, T, operands[1]) else @@ -113,9 +113,8 @@ function _to_polynomial(expr, ::Type{T}) where {T} 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], rev = true) - variables = [x[1] for x in variable_map] + inv = Dict(v => k for (k, v) in d) + variables = [inv[v] for v in MP.variables(poly)] P = MP.polynomial_type(V, T) return ScalarPolynomialFunction{T,P}(poly, variables) end From e141a09f87e79496b137763abb5ef1a5facca7e9 Mon Sep 17 00:00:00 2001 From: votroto Date: Wed, 13 Mar 2024 09:41:53 +0100 Subject: [PATCH 08/16] replace rem+div by divrem --- src/nl_to_polynomial.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/nl_to_polynomial.jl b/src/nl_to_polynomial.jl index b67cbad..c9ed472 100644 --- a/src/nl_to_polynomial.jl +++ b/src/nl_to_polynomial.jl @@ -94,8 +94,9 @@ function _to_polynomial!( return a^b elseif _is_operator(expr, :/) && length(operands) == 2 a, b = _to_polynomial!.(Ref(d), T, operands) - if iszero(MP.rem(a, b)) - return MP.div(a, b) + divisor, remainder = Base.divrem(a, b) + if iszero(remainder) + return divisor else return a / b end From 8af724ef7ff71a85245f6dc25943283515ab1caf Mon Sep 17 00:00:00 2001 From: votroto Date: Wed, 13 Mar 2024 17:09:04 +0100 Subject: [PATCH 09/16] add tests + refactor --- src/QCQP/MOI_wrapper.jl | 57 ++++++------- src/nl_to_polynomial.jl | 6 +- test/qcqp_extra.jl | 174 ++++++++++++++++++++++++++++++++++++++++ 3 files changed, 208 insertions(+), 29 deletions(-) create mode 100644 test/qcqp_extra.jl diff --git a/src/QCQP/MOI_wrapper.jl b/src/QCQP/MOI_wrapper.jl index 42fadd5..6d73953 100644 --- a/src/QCQP/MOI_wrapper.jl +++ b/src/QCQP/MOI_wrapper.jl @@ -211,6 +211,28 @@ function _subs!( ) end +""" + _subs_ensure_moi_order(p::PolyJuMP.ScalarPolynomialFunction, old, new) + +Substitutes old `MP.variables(p.polynomial)` with new vars, while re-sorting the +MOI `p.variables` to get them in the correct order after substitution. +""" +function _subs_ensure_moi_order(p::PolyJuMP.ScalarPolynomialFunction, old, new) + if isempty(old) + return p + end + + poly = MP.subs(p.polynomial, old => new) + + all_new_vars = MP.variables(poly) + to_old_map = Dict(zip(new, old)) + to_moi_map = Dict(zip(MP.variables(p.polynomial), p.variables)) + moi_vars = [to_moi_map[get(to_old_map, v, v)] for v in all_new_vars] + + return PolyJuMP.ScalarPolynomialFunction(poly, moi_vars) +end + + function _subs!( p::PolyJuMP.ScalarPolynomialFunction, index_to_var::Dict{K,V}, @@ -227,16 +249,7 @@ function _subs!( index_to_var[vi] = var end end - if !isempty(old_var) - to_old_map = Dict(zip(new_var, old_var)) - to_moi_map = Dict(zip(MP.variables(p.polynomial), p.variables)) - - poly = MP.subs(p.polynomial, old_var => new_var) - all_new_vars = MP.variables(poly) - moi_vars = [to_moi_map[get(to_old_map, v, v)] for v in all_new_vars] - - p = PolyJuMP.ScalarPolynomialFunction(poly, moi_vars) - end + p = _subs_ensure_moi_order(p, old_var, new_var) return p, index_to_var end @@ -342,23 +355,13 @@ function MOI.Utilities.final_touch(model::Optimizer{T}, _) where {T} 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 + for S in keys(model.constraints) + F = PolyJuMP.ScalarPolynomialFunction{T,model.constraints[S][1]} + for ci in MOI.get(model, MOI.ListOfConstraintIndices{F,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 if !isnothing(monos) diff --git a/src/nl_to_polynomial.jl b/src/nl_to_polynomial.jl index c9ed472..846335d 100644 --- a/src/nl_to_polynomial.jl +++ b/src/nl_to_polynomial.jl @@ -114,8 +114,10 @@ function _to_polynomial(expr, ::Type{T}) where {T} end function _scalar_polynomial(d::Dict{K,V}, ::Type{T}, poly) where {T,K,V} - inv = Dict(v => k for (k, v) in d) - variables = [inv[v] for v in MP.variables(poly)] + var_set = Set(MP.variables(poly)) + variable_map = Tuple{K, V}[(k, v) for (k, v) in d if v in var_set] + 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) end diff --git a/test/qcqp_extra.jl b/test/qcqp_extra.jl new file mode 100644 index 0000000..d749c97 --- /dev/null +++ b/test/qcqp_extra.jl @@ -0,0 +1,174 @@ +module TestQCQPExtra + +using Test + +import MathOptInterface as MOI +import MultivariatePolynomials as MP +import PolyJuMP +import Random + +# Random.seed! set below! + +function test_unconstrained_before_projection(T) + inner = Model{T}() + optimizer = MOI.Utilities.MockOptimizer(inner) + model = PolyJuMP.JuMP.GenericModel{T}( + () -> PolyJuMP.QCQP.Optimizer{T}(optimizer), + ) + PolyJuMP.@variable(model, -1 <= a[1:2] <= 1) + PolyJuMP.@objective(model, Min, a[1]^2*a[2]^2) + 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 test_unconstrained_after_projection(T) + inner = Model{T}() + optimizer = MOI.Utilities.MockOptimizer(inner) + model = PolyJuMP.JuMP.GenericModel{T}( + () -> PolyJuMP.QCQP.Optimizer{T}(optimizer), + ) + PolyJuMP.@variable(model, -1 <= a <= 1) + PolyJuMP.@objective(model, Min, a^2) + PolyJuMP.optimize!(model) + + vis = MOI.get(inner, MOI.ListOfVariableIndices()) + @test length(vis) == 1 + + F = MOI.ScalarQuadraticFunction{T} + S = MOI.EqualTo{T} + cis = MOI.get(inner, MOI.ListOfConstraintIndices{F,S}()) + @test length(cis) == 0 +end + +function _random_polynomial(vars, T) + ms = Random.shuffle(MP.monomials(vars, 1:length(vars))) + return sum(ms[i] * T(randn()) for i in eachindex(ms) if rand(Bool)) +end + +function test_subs!_preserves_moi_sync(xs, ys, T) + p = _random_polynomial(xs, T) + + mois = MOI.VariableIndex.(eachindex(xs)) + vals = T.(randn(length(mois))) + + mask = rand(Bool, length(xs)) + is = Random.shuffle(eachindex(xs)[mask]) + index = Dict{eltype(mois), eltype(xs)}(zip(mois[is], ys[is])) + + moi_map = Dict(zip(xs, mois)) + moivars = [moi_map[v] for v in MP.variables(p)] + + before = PolyJuMP.ScalarPolynomialFunction(p, moivars) + after, _ = PolyJuMP.QCQP._subs!(before, index) + + bmap = [vals[v.value] for v in before.variables] + amap = [vals[v.value] for v in after.variables] + + bvalue = before.polynomial(MP.variables(before.polynomial) => bmap) + avalue = after.polynomial(MP.variables(after.polynomial) => amap) + + # avoid verbose fails + @test isapprox(Float64(bvalue), Float64(avalue)) +end + +function test_scalar_polynomial_function(xs, T) + pick = rand(eachindex(xs)) + ids = Random.shuffle(eachindex(xs)) + poly = sum(T(randn()) * xs[i] for i in ids if i != pick) + mois = MOI.VariableIndex.(eachindex(xs)) + moi_to_vars = Dict(zip(mois, xs)) + + spf = PolyJuMP._scalar_polynomial(moi_to_vars, Any, poly) + + expected = MP.variables(poly) + actual = [moi_to_vars[m] for m in spf.variables] + + @test length(MP.variables(spf.polynomial)) == length(spf.variables) + @test expected == actual +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 run_tests_e2e() + for name in names(@__MODULE__; all = true) + if startswith("$name", "test_e2e") + @testset "$(name) $T" for T in [Int, Float64] + getfield(@__MODULE__, name)(T) + end + end + end +end + +function run_test_scalar_polynomial_function(xs, samples) + @testset "qcqp extra $T" for T in [Float64, BigFloat] + for i in 1:samples + test_scalar_polynomial_function(xs, T) + end + end +end + +function run_tests_subs(xs, ys, samples, TS) + Random.seed!(2024) + for name in names(@__MODULE__; all = true) + if startswith("$name", "test_subs") + @testset "$(name) $T" for T in TS + for i in 1:samples + getfield(@__MODULE__, name)(xs, ys, T) + end + end + end + end +end + +end # module + +using Test + +@testset "TestQCQPFinalTouch" begin + TestQCQPExtra.run_tests_e2e() +end + +import DynamicPolynomials +@testset "DynamicPolynomials" begin + ids = 1:4 + DynamicPolynomials.@polyvar(x[ids]) + DynamicPolynomials.@polyvar(y[ids]) + samples = 10 + types = [Float64, BigFloat] # Rational fails with DynamicPolynomials + TestQCQPExtra.run_tests_subs(x, y, samples, types) + + TestQCQPExtra.run_test_scalar_polynomial_function(x, samples) +end + +import TypedPolynomials +@testset "TypedPolynomials" begin + TypedPolynomials.@polyvar(z[1:4]) + TypedPolynomials.@polyvar(w[1:4]) + types = [Float64, BigFloat, Rational{BigInt}] + samples = 10 + TestQCQPExtra.run_tests_subs(z, w, samples, types) +end \ No newline at end of file From f57eadeda819c327e86a4e557b8afdf263d34e98 Mon Sep 17 00:00:00 2001 From: Oscar Dowson Date: Fri, 15 Mar 2024 13:54:37 +1300 Subject: [PATCH 10/16] Apply suggestions from code review --- src/QCQP/MOI_wrapper.jl | 1 - src/nl_to_polynomial.jl | 2 +- test/qcqp_extra.jl | 10 ++++------ 3 files changed, 5 insertions(+), 8 deletions(-) diff --git a/src/QCQP/MOI_wrapper.jl b/src/QCQP/MOI_wrapper.jl index 6d73953..85f314d 100644 --- a/src/QCQP/MOI_wrapper.jl +++ b/src/QCQP/MOI_wrapper.jl @@ -232,7 +232,6 @@ function _subs_ensure_moi_order(p::PolyJuMP.ScalarPolynomialFunction, old, new) return PolyJuMP.ScalarPolynomialFunction(poly, moi_vars) end - function _subs!( p::PolyJuMP.ScalarPolynomialFunction, index_to_var::Dict{K,V}, diff --git a/src/nl_to_polynomial.jl b/src/nl_to_polynomial.jl index 846335d..89749ac 100644 --- a/src/nl_to_polynomial.jl +++ b/src/nl_to_polynomial.jl @@ -115,7 +115,7 @@ end function _scalar_polynomial(d::Dict{K,V}, ::Type{T}, poly) where {T,K,V} var_set = Set(MP.variables(poly)) - variable_map = Tuple{K, V}[(k, v) for (k, v) in d if v in var_set] + variable_map = Tuple{K,V}[(k, v) for (k, v) in d if v in var_set] sort!(variable_map, by = x -> x[2], rev = true) variables = [x[1] for x in variable_map] P = MP.polynomial_type(V, T) diff --git a/test/qcqp_extra.jl b/test/qcqp_extra.jl index d749c97..b3d6964 100644 --- a/test/qcqp_extra.jl +++ b/test/qcqp_extra.jl @@ -16,12 +16,11 @@ function test_unconstrained_before_projection(T) () -> PolyJuMP.QCQP.Optimizer{T}(optimizer), ) PolyJuMP.@variable(model, -1 <= a[1:2] <= 1) - PolyJuMP.@objective(model, Min, a[1]^2*a[2]^2) + PolyJuMP.@objective(model, Min, a[1]^2 * a[2]^2) 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}()) @@ -35,12 +34,11 @@ function test_unconstrained_after_projection(T) () -> PolyJuMP.QCQP.Optimizer{T}(optimizer), ) PolyJuMP.@variable(model, -1 <= a <= 1) - PolyJuMP.@objective(model, Min, a^2) + PolyJuMP.@objective(model, Min, a^2) PolyJuMP.optimize!(model) vis = MOI.get(inner, MOI.ListOfVariableIndices()) @test length(vis) == 1 - F = MOI.ScalarQuadraticFunction{T} S = MOI.EqualTo{T} cis = MOI.get(inner, MOI.ListOfConstraintIndices{F,S}()) @@ -60,7 +58,7 @@ function test_subs!_preserves_moi_sync(xs, ys, T) mask = rand(Bool, length(xs)) is = Random.shuffle(eachindex(xs)[mask]) - index = Dict{eltype(mois), eltype(xs)}(zip(mois[is], ys[is])) + index = Dict{eltype(mois),eltype(xs)}(zip(mois[is], ys[is])) moi_map = Dict(zip(xs, mois)) moivars = [moi_map[v] for v in MP.variables(p)] @@ -171,4 +169,4 @@ import TypedPolynomials types = [Float64, BigFloat, Rational{BigInt}] samples = 10 TestQCQPExtra.run_tests_subs(z, w, samples, types) -end \ No newline at end of file +end From 481fbf3295e9d542df329308cb0adaf147d3e1d7 Mon Sep 17 00:00:00 2001 From: Oscar Dowson Date: Fri, 15 Mar 2024 20:25:46 +1300 Subject: [PATCH 11/16] Update runtests.jl --- test/runtests.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/test/runtests.jl b/test/runtests.jl index 99b61bb..a5da875 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -8,3 +8,4 @@ include("Mock/mock_tests.jl") include("kkt.jl") include("sage.jl") include("qcqp.jl") +include("qcqp_extra.jl") From 23a69fa6ae4773b1f8f2701f800745585427cb79 Mon Sep 17 00:00:00 2001 From: Oscar Dowson Date: Fri, 15 Mar 2024 20:41:42 +1300 Subject: [PATCH 12/16] Update Project.toml --- test/Project.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/test/Project.toml b/test/Project.toml index 125b4cf..995e4fd 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -10,6 +10,7 @@ MultivariateMoments = "f4abf1af-0426-5881-a0da-e2f168889b5e" MultivariatePolynomials = "102ac46a-7ee4-5c85-9060-abc95bfdeaa3" MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0" PolyJuMP = "ddf597a6-d67e-5340-b84c-e37d84115374" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" SemialgebraicSets = "8e049039-38e8-557d-ae3a-bc521ccf6204" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" TypedPolynomials = "afbbf031-7a57-5f58-a1b9-b774a0fad08d" From f36419d626767ed00ffa8139c548eb36bafd2318 Mon Sep 17 00:00:00 2001 From: Oscar Dowson Date: Mon, 18 Mar 2024 09:59:56 +1300 Subject: [PATCH 13/16] Update qcqp_extra.jl --- test/qcqp_extra.jl | 121 +++++++++++++++++++++------------------------ 1 file changed, 57 insertions(+), 64 deletions(-) diff --git a/test/qcqp_extra.jl b/test/qcqp_extra.jl index b3d6964..c345886 100644 --- a/test/qcqp_extra.jl +++ b/test/qcqp_extra.jl @@ -7,7 +7,59 @@ import MultivariatePolynomials as MP import PolyJuMP import Random -# Random.seed! set below! + +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 run_tests_e2e() + for name in names(@__MODULE__; all = true) + if startswith("$name", "test_e2e") + @testset "$(name) $T" for T in [Int, Float64] + getfield(@__MODULE__, name)(T) + end + end + end + return +end + +function run_test_scalar_polynomial_function(xs, samples) + @testset "qcqp extra $T" for T in [Float64, BigFloat] + for i in 1:samples + test_scalar_polynomial_function(xs, T) + end + end + return +end + +function run_tests_subs(xs, ys, samples, TS) + Random.seed!(2024) + for name in names(@__MODULE__; all = true) + if startswith("$name", "test_subs") + @testset "$(name) $T" for T in TS + for i in 1:samples + getfield(@__MODULE__, name)(xs, ys, T) + end + end + end + end + return +end function test_unconstrained_before_projection(T) inner = Model{T}() @@ -18,13 +70,13 @@ function test_unconstrained_before_projection(T) PolyJuMP.@variable(model, -1 <= a[1:2] <= 1) PolyJuMP.@objective(model, Min, a[1]^2 * a[2]^2) 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 + return end function test_unconstrained_after_projection(T) @@ -36,13 +88,13 @@ function test_unconstrained_after_projection(T) PolyJuMP.@variable(model, -1 <= a <= 1) PolyJuMP.@objective(model, Min, a^2) PolyJuMP.optimize!(model) - vis = MOI.get(inner, MOI.ListOfVariableIndices()) @test length(vis) == 1 F = MOI.ScalarQuadraticFunction{T} S = MOI.EqualTo{T} cis = MOI.get(inner, MOI.ListOfConstraintIndices{F,S}()) @test length(cis) == 0 + return end function _random_polynomial(vars, T) @@ -52,28 +104,22 @@ end function test_subs!_preserves_moi_sync(xs, ys, T) p = _random_polynomial(xs, T) - mois = MOI.VariableIndex.(eachindex(xs)) vals = T.(randn(length(mois))) - mask = rand(Bool, length(xs)) is = Random.shuffle(eachindex(xs)[mask]) index = Dict{eltype(mois),eltype(xs)}(zip(mois[is], ys[is])) - moi_map = Dict(zip(xs, mois)) moivars = [moi_map[v] for v in MP.variables(p)] - before = PolyJuMP.ScalarPolynomialFunction(p, moivars) after, _ = PolyJuMP.QCQP._subs!(before, index) - bmap = [vals[v.value] for v in before.variables] amap = [vals[v.value] for v in after.variables] - bvalue = before.polynomial(MP.variables(before.polynomial) => bmap) avalue = after.polynomial(MP.variables(after.polynomial) => amap) - # avoid verbose fails @test isapprox(Float64(bvalue), Float64(avalue)) + return end function test_scalar_polynomial_function(xs, T) @@ -82,64 +128,12 @@ function test_scalar_polynomial_function(xs, T) poly = sum(T(randn()) * xs[i] for i in ids if i != pick) mois = MOI.VariableIndex.(eachindex(xs)) moi_to_vars = Dict(zip(mois, xs)) - spf = PolyJuMP._scalar_polynomial(moi_to_vars, Any, poly) - expected = MP.variables(poly) actual = [moi_to_vars[m] for m in spf.variables] - @test length(MP.variables(spf.polynomial)) == length(spf.variables) @test expected == actual -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 run_tests_e2e() - for name in names(@__MODULE__; all = true) - if startswith("$name", "test_e2e") - @testset "$(name) $T" for T in [Int, Float64] - getfield(@__MODULE__, name)(T) - end - end - end -end - -function run_test_scalar_polynomial_function(xs, samples) - @testset "qcqp extra $T" for T in [Float64, BigFloat] - for i in 1:samples - test_scalar_polynomial_function(xs, T) - end - end -end - -function run_tests_subs(xs, ys, samples, TS) - Random.seed!(2024) - for name in names(@__MODULE__; all = true) - if startswith("$name", "test_subs") - @testset "$(name) $T" for T in TS - for i in 1:samples - getfield(@__MODULE__, name)(xs, ys, T) - end - end - end - end + return end end # module @@ -158,7 +152,6 @@ import DynamicPolynomials samples = 10 types = [Float64, BigFloat] # Rational fails with DynamicPolynomials TestQCQPExtra.run_tests_subs(x, y, samples, types) - TestQCQPExtra.run_test_scalar_polynomial_function(x, samples) end From d125564902511d6dee5daf3604d07d7055d051c0 Mon Sep 17 00:00:00 2001 From: Oscar Dowson Date: Mon, 18 Mar 2024 10:01:04 +1300 Subject: [PATCH 14/16] Update MOI_wrapper.jl --- src/QCQP/MOI_wrapper.jl | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/src/QCQP/MOI_wrapper.jl b/src/QCQP/MOI_wrapper.jl index 85f314d..82b455f 100644 --- a/src/QCQP/MOI_wrapper.jl +++ b/src/QCQP/MOI_wrapper.jl @@ -221,14 +221,11 @@ function _subs_ensure_moi_order(p::PolyJuMP.ScalarPolynomialFunction, old, new) if isempty(old) return p end - poly = MP.subs(p.polynomial, old => new) - all_new_vars = MP.variables(poly) to_old_map = Dict(zip(new, old)) to_moi_map = Dict(zip(MP.variables(p.polynomial), p.variables)) moi_vars = [to_moi_map[get(to_old_map, v, v)] for v in all_new_vars] - return PolyJuMP.ScalarPolynomialFunction(poly, moi_vars) end @@ -354,8 +351,8 @@ function MOI.Utilities.final_touch(model::Optimizer{T}, _) where {T} vars = _add_variables!(func, vars) monos = _add_monomials!(func, monos) end - for S in keys(model.constraints) - F = PolyJuMP.ScalarPolynomialFunction{T,model.constraints[S][1]} + for (S, constraints) in model.constraints + F = PolyJuMP.ScalarPolynomialFunction{T,constraints[1]} for ci in MOI.get(model, MOI.ListOfConstraintIndices{F,S}()) func = MOI.get(model, MOI.ConstraintFunction(), ci) func, index_to_var = _subs!(func, index_to_var) From f0c5940c72917184f75f3cc08c21f62645885fe6 Mon Sep 17 00:00:00 2001 From: Oscar Dowson Date: Mon, 18 Mar 2024 10:17:12 +1300 Subject: [PATCH 15/16] Update test/qcqp_extra.jl --- test/qcqp_extra.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/test/qcqp_extra.jl b/test/qcqp_extra.jl index c345886..22b0c1e 100644 --- a/test/qcqp_extra.jl +++ b/test/qcqp_extra.jl @@ -7,7 +7,6 @@ import MultivariatePolynomials as MP import PolyJuMP import Random - MOI.Utilities.@model( Model, (), From 2a194a11f30b92b0d95345379a13dff091efe510 Mon Sep 17 00:00:00 2001 From: Oscar Dowson Date: Mon, 18 Mar 2024 10:21:07 +1300 Subject: [PATCH 16/16] Update nl_to_polynomial.jl --- src/nl_to_polynomial.jl | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/src/nl_to_polynomial.jl b/src/nl_to_polynomial.jl index 89749ac..d19799e 100644 --- a/src/nl_to_polynomial.jl +++ b/src/nl_to_polynomial.jl @@ -94,12 +94,14 @@ function _to_polynomial!( return a^b elseif _is_operator(expr, :/) && length(operands) == 2 a, b = _to_polynomial!.(Ref(d), T, operands) - divisor, remainder = Base.divrem(a, b) - if iszero(remainder) - return divisor - else - return a / b - end + return a / b + # TODO(odow): see PR#113 + # divisor, remainder = Base.divrem(a, b) + # if iszero(remainder) + # return divisor + # else + # return a / b + # end elseif _is_variable(expr) return _to_polynomial!(d, T, operands[1]) else