Change interpretation and default maxdegree (#227)

* Change interpretation and default maxdegree

* Fix

* Add default_maxdegree to docs

* Add tolerance to test in doc

* Fix for Julia v1.0

Author: blegat

docs/src/constraints.md

@@ -248,6 +248,11 @@ ArXiv e-prints, **2017**.
 
 ## Reference
 
+Default choice for the `maxdegree` keyword:
+```@docs
+SumOfSquares.default_maxdegree
+```
+
 Special case that is second-order cone representable:
 ```@docs
 SumOfSquares.PositiveSemidefinite2x2ConeTriangle

docs/src/tutorials/Symmetry/permutation_symmetry.jl

@@ -43,7 +43,7 @@ model = Model(solver)
 pattern = Symmetry.Pattern(G, Symmetry.VariablePermutation())
 con_ref = @constraint(model, poly - t in SOSCone(), symmetry = pattern)
 optimize!(model)
-@test value(t) ≈ -1 #src
+@test value(t) ≈ -1 rtol=1e-6 #src
 value(t)
 
 # We indeed find `-1`, let's verify that symmetry was exploited:

src/Certificate/Certificate.jl

@@ -97,12 +97,7 @@ end
 function maxdegree_gram_basis(B::Type, variables, maxdegree::Int)
     return MB.maxdegree_basis(B, variables, div(maxdegree, 2))
 end
-function multiplier_maxdegree(maxdegree, q)
-    maxdegree_s2 = maxdegree - MP.maxdegree(q)
-    # If maxdegree_s2 is odd, `div(maxdegree_s2, 2)` would make s^2 have degree up to maxdegree_s2-1
-    # for this reason, we take `div(maxdegree_s2 + 1, 2)` so that s^2 have degree up to maxdegree_s2+1
-    return maxdegree_s2 + 1
-end
+multiplier_maxdegree(maxdegree, q) = maxdegree - MP.maxdegree(q)
 function get(certificate::Putinar, ::MultiplierBasis, index::PreorderIndex, domain::DomainWithVariables)
     q = domain.domain.p[index.value]
     vars = sort!([domain.variables..., MP.variables(q)...], rev = true)

src/constraint.jl

@@ -129,13 +129,41 @@ function default_certificate(::BasicSemialgebraicSet, ::Certificate.NoSparsity,
     return Certificate.Putinar(
         ideal_certificate, cone, basis, maxdegree)
 end
+
+# Julia v1.0 does not support `init` keyword
+#_max_maxdegree(p) = maximum(MP.maxdegree, p, init=0)
+_max_maxdegree(p) = mapreduce(MP.maxdegree, max, p, init=0)
+
+_maxdegree(domain) = 0
+
+function _maxdegree(domain::AlgebraicSet)
+    return _max_maxdegree(domain.I.p)
+end
+
+function _maxdegree(domain::BasicSemialgebraicSet)
+    return max(_max_maxdegree(domain.p), _maxdegree(domain.V))
+end
+
+"""
+    default_maxdegree(monos, domain)
+
+Return the default `maxdegree` to use for certifying a polynomial with
+monomials `monos` to be Sum-of-Squares over the domain `domain`.
+By default, the maximum of the maxdegree of `monos` and of all multipliers
+of the domain are used so that at least constant multipliers can be used
+with a Putinar certificate.
+"""
+function default_maxdegree(monos, domain)
+    return max(MP.maxdegree(monos), _maxdegree(domain))
+end
+
 function JuMP.moi_set(
     cone::SOSLikeCone,
     monos::AbstractVector{<:MP.AbstractMonomial};
     domain::AbstractSemialgebraicSet=FullSpace(),
     basis=MonomialBasis,
     newton_polytope::Tuple=tuple(),
-    maxdegree::Union{Nothing, Int}=MP.maxdegree(monos),
+    maxdegree::Union{Nothing, Int}=default_maxdegree(monos, domain),
     sparsity::Certificate.Sparsity=Certificate.NoSparsity(),
     symmetry::Union{Nothing,Certificate.Symmetry.Pattern}=nothing,
     newton_of_remainder::Bool=false,

test/Tests/chebyshev.jl

@@ -21,8 +21,8 @@ function chebyshev_test(optimizer, config::MOIT.Config)
     p = p1 + γ * x^ndeg # the leading coeff of p is γ
 
     dom = @set x >= -1 && x <= 1
-    @constraint(model, 1 - p in SOSCone(), domain = dom)
-    @constraint(model, p + 1 in SOSCone(), domain = dom)
+    @constraint(model, 1 - p in SOSCone(), domain = dom, maxdegree = 9)
+    @constraint(model, p + 1 in SOSCone(), domain = dom, maxdegree = 9)
 
     @objective(model, Max, γ)
 

test/Tests/rearrangement.jl

@@ -18,7 +18,7 @@ function rearrangement_test(optimizer, config::MOIT.Config)
     S = @set y ≤ x && y ≤ z && x ≥ 0 && z ≥ 0
     # In fact, it is nonnegative everywhere so let's add terms `x + z`:
     cref = @constraint(model, x^2 + z^2 + 2y^2 - 2x * y - 2y * z + x + z in SOSCone(),
-                       domain = S, sparsity = VariableSparsity())
+                       domain = S, sparsity = VariableSparsity(), maxdegree = 3)
 
     optimize!(model)
 

test/equalitypolyconstr.jl

@@ -3,7 +3,7 @@
 
     m = SOSModel(factory)
 
-    @constraint(m, x^3 + x*y^2 == x, domain=(@set x^2 + y^2 >= 1 && x^2 + y^2 <= 1))
+    @constraint(m, x^3 + x*y^2 == x, domain=(@set x^2 + y^2 >= 1 && x^2 + y^2 <= 1), maxdegree=4)
 
     JuMP.optimize!(m)