Fix constraint with less variables than variety (#239)

* Fix constraint with less variables than variety

* Fixes

* Fix with_variables

* Fixes

docs/src/tutorials/Extension/certificate.jl

@@ -54,30 +54,28 @@ end
 SOSC.cone(certificate::Schmüdgen) = certificate.cone
 
 function SOSC.preprocessed_domain(::Schmüdgen, domain::BasicSemialgebraicSet, p)
-    return SOSC.DomainWithVariables(domain, variables(p))
+    return SOSC.with_variables(domain, p)
 end
 
-function SOSC.preorder_indices(::Schmüdgen, domain::SOSC.DomainWithVariables)
-    n = length(domain.domain.p)
+function SOSC.preorder_indices(::Schmüdgen, domain::SOSC.WithVariables)
+    n = length(domain.inner.p)
     if n >= Sys.WORD_SIZE
         error("There are $(2^n - 1) products in Schmüdgen's certificate, they cannot even be indexed with `$Int`.")
     end
     return map(SOSC.PreorderIndex, 1:(2^n-1))
 end
 
-function SOSC.multiplier_basis(certificate::Schmüdgen, index::SOSC.PreorderIndex, domain::SOSC.DomainWithVariables)
+function SOSC.multiplier_basis(certificate::Schmüdgen, index::SOSC.PreorderIndex, domain::SOSC.WithVariables)
     q = SOSC.generator(certificate, index, domain)
-    vars = sort!([domain.variables..., variables(q)...], rev = true)
-    unique!(vars)
-    return SOSC.maxdegree_gram_basis(certificate.basis, vars, SOSC.multiplier_maxdegree(certificate.maxdegree, q))
+    return SOSC.maxdegree_gram_basis(certificate.basis, variables(domain), SOSC.multiplier_maxdegree(certificate.maxdegree, q))
 end
 function SOSC.multiplier_basis_type(::Type{Schmüdgen{IC, CT, BT}}) where {IC, CT, BT}
     return BT
 end
 
-function SOSC.generator(::Schmüdgen, index::SOSC.PreorderIndex, domain::SOSC.DomainWithVariables)
-    I = [i for i in eachindex(domain.domain.p) if !iszero(index.value & (1 << (i - 1)))]
-    return prod([domain.domain.p[i] for i in eachindex(domain.domain.p) if !iszero(index.value & (1 << (i - 1)))])
+function SOSC.generator(::Schmüdgen, index::SOSC.PreorderIndex, domain::SOSC.WithVariables)
+    I = [i for i in eachindex(domain.inner.p) if !iszero(index.value & (1 << (i - 1)))]
+    return prod([domain.inner.p[i] for i in eachindex(domain.inner.p) if !iszero(index.value & (1 << (i - 1)))])
 end
 
 SOSC.ideal_certificate(certificate::Schmüdgen) = certificate.ideal_certificate

docs/src/tutorials/Extension/hypercube.jl

@@ -72,6 +72,8 @@ end
 struct HypercubeSet{V} <: SS.AbstractAlgebraicSet
     ideal::HypercubeIdeal{V}
 end
+MP.variables(set::HypercubeSet) = MP.variables(set.ideal)
+MP.variables(ideal::HypercubeIdeal) = ideal.variables
 MP.changecoefficienttype(set::HypercubeSet, ::Type) = set
 SS.ideal(set::HypercubeSet) = set.ideal
 function Base.rem(p, set::HypercubeIdeal)

docs/src/tutorials/Getting started/circle.jl

@@ -0,0 +1,39 @@
+# # Nonnegative over a variety
+
+#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/generated/Getting started/circle.ipynb)
+#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/Getting started/circle.ipynb)
+
+# The polynomial ``1 - y^2`` is nonnegative for all ``y`` in the unit circle.
+# This can be verified using Sum-of-Squares.
+
+using Test #src
+using DynamicPolynomials
+using SumOfSquares
+@polyvar x y
+S = @set x^2 + y^2 == 1
+
+# We need to pick an SDP solver, see [here](https://jump.dev/JuMP.jl/v0.23.0/installation/#Supported-solvers) for a list of the available choices.
+# The domain over which the nonnegativity of ``1 - y^2`` should be certified
+# is specified through the `domain` keyword argument.
+
+import CSDP
+model = SOSModel(CSDP.Optimizer)
+set_silent(model)
+con_ref = @constraint(model, 1 - y^2 >= 0, domain = S)
+optimize!(model)
+
+# We can see that the model was feasible:
+
+@test JuMP.termination_status(model) == MOI.OPTIMAL #src
+@test JuMP.primal_status(model) == MOI.FEASIBLE_POINT #src
+solution_summary(model)
+
+# The certificate can be obtained as follows:
+
+dec = sos_decomposition(con_ref, 1e-6) #src
+@test length(dec) == 1 #src
+@test first(dec) ≈ x rtol = 1e-6 #src
+sos_decomposition(con_ref, 1e-6)
+
+# It returns ``x^2`` which is a valid certificate as:
+# $$ 1 - y^2 \equiv x^2 \pmod{x^2 + y^2 - 1} $$

src/Bridges/Constraint/sos_polynomial.jl

@@ -40,7 +40,10 @@ function MOI.Bridges.Constraint.bridge_constraint(
     # `set.domain.V` is `FullSpace` or `FixedPolynomialsSet`.
     # FIXME convert needed because the coefficient type of `r` is `Any` otherwise if `domain` is `AlgebraicSet`
     r = SOS.Certificate.reduced_polynomial(s.certificate, p, MP.changecoefficienttype(s.domain, T))
-    gram_basis = SOS.Certificate.gram_basis(s.certificate, r)
+    gram_basis = SOS.Certificate.gram_basis(
+        s.certificate,
+        SOS.Certificate.with_variables(r, s.domain),
+    )
     g, Q, cQ = SOS.add_gram_matrix(model, MCT, gram_basis, T)
     # MOI does not modify the coefficients of the functions so we can modify `r`.
     # without altering `f`.

src/Certificate/Certificate.jl

@@ -78,24 +78,47 @@ end
 
 cone(certificate::Putinar) = certificate.cone
 
-struct DomainWithVariables{S, V}
-    domain::S
+struct WithVariables{S, V}
+    inner::S
     variables::V
 end
+
+function MP.variables(v::WithVariables)
+    return v.variables
+end
+function MP.monomials(v::WithVariables)
+    return MP.monomials(v.inner)
+end
+
+
+_cat(a::Vector, b::Vector) = vcat(a, b)
+_cat(a::Vector, ::Tuple{}) = copy(a)
+_cat(a::Tuple, b::Tuple) = [a..., b...]
+_cat(a::Tuple, ::Tuple{}) = [a...]
+
+_vars(::SemialgebraicSets.FullSpace) = tuple()
+_vars(x) = MP.variables(x)
+
+function with_variables(inner, outer)
+    vars = sort!(_cat(_vars(inner), _vars(outer)), rev = true)
+    unique!(vars)
+    return WithVariables(inner, vars)
+end
+
 function preprocessed_domain(::Putinar, domain::BasicSemialgebraicSet, p)
-    return DomainWithVariables(domain, MP.variables(p))
+    return with_variables(domain, p)
 end
 
-function preorder_indices(::Putinar, domain::DomainWithVariables)
-    return map(PreorderIndex, eachindex(domain.domain.p))
+function preorder_indices(::Putinar, domain::WithVariables)
+    return map(PreorderIndex, eachindex(domain.inner.p))
 end
 
 function maxdegree_gram_basis(B::Type, variables, maxdegree::Int)
     return MB.maxdegree_basis(B, variables, div(maxdegree, 2))
 end
 multiplier_maxdegree(maxdegree, q) = maxdegree - MP.maxdegree(q)
-function multiplier_basis(certificate::Putinar, index::PreorderIndex, domain::DomainWithVariables)
-    q = domain.domain.p[index.value]
+function multiplier_basis(certificate::Putinar, index::PreorderIndex, domain::WithVariables)
+    q = domain.inner.p[index.value]
     vars = sort!([domain.variables..., MP.variables(q)...], rev = true)
     unique!(vars)
     return maxdegree_gram_basis(certificate.basis, vars, multiplier_maxdegree(certificate.maxdegree, q))
@@ -104,8 +127,8 @@ function multiplier_basis_type(::Type{Putinar{IC, CT, BT}}) where {IC, CT, BT}
     return BT
 end
 
-function generator(::Putinar, index::PreorderIndex, domain::DomainWithVariables)
-    return domain.domain.p[index.value]
+function generator(::Putinar, index::PreorderIndex, domain::WithVariables)
+    return domain.inner.p[index.value]
 end
 
 ideal_certificate(certificate::Putinar) = certificate.ideal_certificate

src/Certificate/Sparsity/ideal.jl

@@ -35,6 +35,9 @@ end
 function sparsity(poly::MP.AbstractPolynomial, sp::Union{SignSymmetry, Monomial}, certificate::SumOfSquares.Certificate.AbstractIdealCertificate)
     return sparsity(MP.monomials(poly), sp, SumOfSquares.Certificate.gram_basis(certificate, poly))
 end
+function sparsity(v::SumOfSquares.Certificate.WithVariables, sp, certificate)
+    return sparsity(v.inner, sp, certificate)
+end
 function SumOfSquares.Certificate.gram_basis(certificate::Ideal, poly)
     return sparsity(poly, certificate.sparsity, certificate.certificate)
 end

test/certificate.jl

@@ -4,6 +4,22 @@ const MP = MultivariatePolynomials
 import MultivariateBases
 const MB = MultivariateBases
 
+@testset "with_variables" begin
+    @polyvar x y z
+    p = x + z
+    v = SumOfSquares.Certificate.with_variables(p, y)
+    @test v.inner === p
+    @test MP.variables(v.inner) == [x, z]
+    @test MP.variables(v) == [x, y, z]
+
+    @ncpolyvar a b
+    q = a * b + b * a
+    v = SumOfSquares.Certificate.with_variables(q, FullSpace())
+    @test v.inner === q
+    @test MP.variables(v.inner) == [a, b, a]
+    @test MP.variables(v) == [a, b]
+end
+
 @testset "Monomial selection for certificate" begin
     @polyvar x y z
     @ncpolyvar a b