Fixes

.github/workflows/documentation.yml

@@ -22,6 +22,7 @@ jobs:
         run: |
           using Pkg
           Pkg.add([
+              PackageSpec(url="https://github.com/blegat/HomotopyContinuation.jl", rev="dynamicpolynomials_v0.6"),
               PackageSpec(name="MathOptInterface", rev="bl/set_map_type"),
               PackageSpec(name="StarAlgebras", rev="main"),
               PackageSpec(name="SymbolicWedderburn", rev="mk/StarAlgebras_0_3"),

docs/Project.toml

@@ -29,9 +29,11 @@ PolyJuMP = "ddf597a6-d67e-5340-b84c-e37d84115374"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 SCS = "c946c3f1-0d1f-5ce8-9dea-7daa1f7e2d13"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
+StarAlgebras = "0c0c59c1-dc5f-42e9-9a8b-b5dc384a6cd1"
 SumOfSquares = "4b9e565b-77fc-50a5-a571-1244f986bda1"
 SymbolicWedderburn = "858aa9a9-4c7c-4c62-b466-2421203962a2"
 TypedPolynomials = "afbbf031-7a57-5f58-a1b9-b774a0fad08d"
 
 [compat]
 Documenter = "1"
+DynamicPolynomials = "0.6"

docs/src/tutorials/Extension/certificate.jl

@@ -50,6 +50,7 @@ solution_summary(model)
 
 # We now define the Schmüdgen's certificate:
 
+import StarAlgebras as SA
 import MultivariateBases as MB
 const SOS = SumOfSquares
 const SOSC = SOS.Certificate
@@ -79,7 +80,7 @@ function SOSC.multiplier_basis(certificate::Schmüdgen, index::SOSC.PreorderInde
     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}}, ::Type{M}) where {IC,CT,BT,M}
-    return MB.similar_type(BT, M)
+    return MB.explicit_basis_type(BT)
 end
 
 function SOSC.generator(::Schmüdgen, index::SOSC.PreorderIndex, domain::SOSC.WithVariables)
@@ -97,8 +98,9 @@ SOS.matrix_cone_type(::Type{<:Schmüdgen{IC, CT}}) where {IC, CT} = SOS.matrix_c
 model = SOSModel(solver)
 @variable(model, α)
 @objective(model, Max, α)
-ideal_certificate = SOSC.Newton(SOSCone(), MB.MonomialBasis, tuple())
-certificate = Schmüdgen(ideal_certificate, SOSCone(), MB.MonomialBasis, maxdegree(p))
+basis = MB.FullBasis{MB.Monomial,typeof(x * y)}()
+ideal_certificate = SOSC.Newton(SOSCone(), basis, tuple())
+certificate = Schmüdgen(ideal_certificate, SOSCone(), basis, maxdegree(p))
 @constraint(model, c, p >= α, domain = S, certificate = certificate)
 optimize!(model)
 @test termination_status(model) == MOI.OPTIMAL #src

docs/src/tutorials/Extension/hypercube.jl

@@ -64,6 +64,7 @@ S.I.algo
 # However, we still need to divide polynomials by the Gröbner basis
 # which can be simplified in this case.
 
+import MutableArithmetics as MA
 const MP = MultivariatePolynomials
 const SS = SemialgebraicSets
 struct HypercubeIdeal{V} <: SS.AbstractPolynomialIdeal
@@ -76,6 +77,10 @@ MP.variables(set::HypercubeSet) = MP.variables(set.ideal)
 MP.variables(ideal::HypercubeIdeal) = ideal.variables
 Base.similar(set::HypercubeSet, ::Type) = set
 SS.ideal(set::HypercubeSet) = set.ideal
+function MA.promote_operation(::typeof(SS.ideal), ::Type{HypercubeSet{V}}) where {V}
+    return HypercubeIdeal{V}
+end
+SS.similar_type(S::Type{<:HypercubeSet}, ::Type) = S
 function Base.rem(p, set::HypercubeIdeal)
     return MP.polynomial(map(MP.terms(p)) do term
         mono = MP.monomial(term)

docs/src/tutorials/Extension/univariate_solver.jl

@@ -19,6 +19,7 @@ module MyUnivariateSolver
 import LinearAlgebra
 import MathOptInterface as MOI
 import MultivariatePolynomials as MP
+import MultivariateBases as MB
 import SumOfSquares as SOS
 
 function decompose(p::MP.AbstractPolynomial, tol=1e-6)
@@ -56,7 +57,7 @@ function decompose(p::MP.AbstractPolynomial, tol=1e-6)
     end
     q1 = MP.map_coefficients(real, q)
     q2 = MP.map_coefficients(imag, q)
-    return SOS.SOSDecomposition([q1, q2])
+    return SOS.SOSDecomposition([MB.algebra_element(q1), MB.algebra_element(q2)])
 end
 
 mutable struct Optimizer <: MOI.AbstractOptimizer
@@ -84,7 +85,7 @@ function MOI.add_constraint(optimizer::Optimizer, func::MOI.VectorAffineFunction
     if !isempty(func.terms)
         error("Only supports constant polynomials")
     end
-    optimizer.p = MP.polynomial(func.constants, set.monomials)
+    optimizer.p = MP.polynomial(MB.algebra_element(func.constants, set.basis))
     return MOI.ConstraintIndex{typeof(func),typeof(set)}(1) # There will be only ever one constraint so the index does not matter.
 end
 

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

@@ -40,7 +40,7 @@ Q = value_matrix(q) #src
 
 sosdec = SOSDecomposition(q)
 @test isapprox(sosdec, sos_decomposition(con_ref)) #src
-@test isapprox(sum(sosdec.ps.^2), p; rtol=1e-4, ztol=1e-6) #src
+@test isapprox(polynomial(sosdec), p; rtol=1e-4, ztol=1e-6) #src
 
 # We now seek for the SOS decomposition of the following polynomial:
 

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

@@ -76,7 +76,7 @@ Q12 = η1.Q * η.atoms[1].weight + η2.Q * η.atoms[2].weight
 
 model = SOSModel(CSDP.Optimizer)
 @variable(model, σ)
-@constraint(model, cheby_cref, p >= σ, basis = ChebyshevBasisFirstKind)
+@constraint(model, cheby_cref, p >= σ, basis = Chebyshev)
 @objective(model, Max, σ)
 optimize!(model)
 solution_summary(model)

docs/src/tutorials/Noncommutative and Hermitian/noncommutative_variables.jl

@@ -14,6 +14,7 @@
 # $(x * y + x^2)^2 = x^4 + x^3y + xyx^2 + xyxy$ is tested to be sum-of-squares.
 
 using Test #src
+import StarAlgebras as SA #src
 using DynamicPolynomials
 @ncpolyvar x y
 p = (x * y + x^2)^2
@@ -45,7 +46,7 @@ sos_decomposition(con_ref)
 
 dec = sos_decomposition(con_ref, 1e-6)                    #src
 @test length(dec.ps) == 1                                 #src
-@test sign(first(coefficients(dec.ps[1]))) * dec.ps[1] ≈ x * y + x^2 rtol=1e-5 atol=1e-5 #src
+@test sign(first(SA.coeffs(dec.ps[1]))) * dec.ps[1] ≈ x * y + x^2 rtol=1e-5 atol=1e-5 #src
 sos_decomposition(con_ref, 1e-6)
 
 # ## Example 2.2

docs/src/tutorials/Noncommutative and Hermitian/sums_of_hermitian_squares.jl

@@ -5,6 +5,7 @@
 # **Contributed by**: Benoît Legat
 
 using Test #src
+import StarAlgebras as SA #src
 
 using SumOfSquares
 
@@ -19,6 +20,6 @@ con_ref = @constraint(model, p in cone)
 optimize!(model)
 dec = sos_decomposition(con_ref, 1e-6) #src
 @test length(dec.ps) == 1 #src
-@test dec.ps[1]' * dec.ps[1] ≈ p atol=1e-6 rtol=1e-6 #src
-@test sign(real(first(coefficients(dec.ps[1])))) * dec.ps[1] ≈ y + im * x atol=1e-6 rtol=1e-6 #src
+@test polynomial(dec.ps[1])' * polynomial(dec.ps[1]) ≈ p atol=1e-6 rtol=1e-6 #src
+@test sign(real(first(SA.coeffs(dec.ps[1])))) * dec.ps[1] ≈ y + im * x atol=1e-6 rtol=1e-6 #src
 sos_decomposition(con_ref, 1e-6)

src/Bridges/Constraint/sos_polynomial.jl

@@ -175,15 +175,34 @@ end
 
 function _get(
     model::MOI.ModelLike,
-    attr,
+    attr::SOS.SOSDecompositionAttribute,
+    constraint::MOI.ConstraintIndex,
+    index::Int,
+)
+    return MOI.get(
+        model,
+        typeof(attr)(;
+            ranktol = attr.ranktol,
+            dec = attr.dec,
+            multiplier_index = index,
+            result_index = attr.result_index,
+        ),
+        constraint,
+    )
+end
+
+function _get(
+    model::MOI.ModelLike,
+    attr::Union{
+        SOS.GramMatrixAttribute,
+        SOS.MomentMatrixAttribute,
+    },
     constraint::MOI.ConstraintIndex,
     index::Int,
 )
     return MOI.get(
         model,
         typeof(attr)(;
-            attr.ranktol,
-            attr.dec,
             multiplier_index = index,
             result_index = attr.result_index,
         ),

src/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl

@@ -286,6 +286,17 @@ function MOI.get(
 )
     return MOI.get(model, attr, bridge.constraint)
 end
+
+function _gram(
+    f::Function,
+    Q::Vector{MOI.VariableIndex},
+    gram_basis,
+    T::Type,
+    MCT,
+)
+    return SOS.build_gram_matrix(convert(Vector{T}, f(Q)), gram_basis, MCT, T)
+end
+
 function _gram(
     f::Function,
     Qs::Vector{Vector{MOI.VariableIndex}},

src/Certificate/Certificate.jl

@@ -49,9 +49,13 @@ function within_total_bounds(mono::MP.AbstractMonomial, bounds::DegreeBounds)
     return bounds.mindegree <= MP.degree(mono) <= bounds.maxdegree
 end
 
+_vec(v::AbstractVector) = v
+# For `TypedPolynomials`
+_vec(v::Tuple) = MP.variable_union_type(first(v))[v...]
+
 function _divides(a, b)
     # `MP.divides(a, b)` is not implemented yet for noncommutative
-    vars = unique!(sort(MP.variables(a)))
+    vars = unique!(sort(_vec(MP.variables(a))))
     comm = is_commutative(vars)
     return all(vars) do v
         return _degree(a, v, comm) <= _degree(b, v, comm)

src/Certificate/ideal.jl

@@ -180,32 +180,8 @@ function Newton(cone, basis, variable_groups::Tuple)
     )
 end
 
-function _weight_type(::Type{T}, ::Type{BT}) where {T,BT}
-    return SA.AlgebraElement{
-        MA.promote_operation(
-            MB.algebra,
-            MA.promote_operation(MB.implicit_basis, BT),
-        ),
-        T,
-        MA.promote_operation(
-            MB.sparse_coefficients,
-            MP.polynomial_type(MP.monomial_type(BT), T),
-        ),
-    }
-end
-
-function _half_newton_polytope(a::SA.AlgebraElement, vars, filter)
-    return half_newton_polytope(
-        a,
-        _weight_type(Bool, typeof(SA.basis(a)))[],
-        vars,
-        _maxdegree(a, vars),
-        filter,
-    )[1]
-end
-
 function gram_basis(certificate::Newton, poly)
-    return _half_newton_polytope(
+    return half_newton_polytope(
         _algebra_element(poly),
         MP.variables(poly),
         certificate.newton,

src/Certificate/newton_polytope.jl

@@ -807,3 +807,46 @@ end
 function _sub(basis::SubBasis{B}, I) where {B}
     return SubBasis{B}(basis.monomials[I])
 end
+
+function _weight_type(::Type{T}, ::Type{BT}) where {T,BT}
+    return SA.AlgebraElement{
+        MA.promote_operation(
+            MB.algebra,
+            MA.promote_operation(MB.implicit_basis, BT),
+        ),
+        T,
+        MA.promote_operation(
+            MB.sparse_coefficients,
+            MP.polynomial_type(MP.monomial_type(BT), T),
+        ),
+    }
+end
+
+function half_newton_polytope(a::SA.AlgebraElement, vars, filter)
+    return half_newton_polytope(
+        a,
+        _weight_type(Bool, typeof(SA.basis(a)))[],
+        vars,
+        _maxdegree(a, vars),
+        filter,
+    )[1]
+end
+
+function half_newton_polytope(a::SA.AlgebraElement, filter)
+    return half_newton_polytope(a, MP.variables(a), filter)
+end
+
+function half_newton_polytope(basis::MB.SubBasis, args...)
+    a = MB.algebra_element(
+        SA.SparseCoefficients(basis.monomials, ones(length(basis))),
+        MB.implicit_basis(basis),
+    )
+    return half_newton_polytope(a, args...)
+end
+
+function monomials_half_newton_polytope(monos, args...)
+    return half_newton_polytope(
+        MB.SubBasis{MB.Monomial}(monos),
+        args...,
+    ).monomials
+end

src/rand.jl

@@ -8,27 +8,14 @@ function randpsd(n; r = n, eps = 0.1)
     return Q' * Diagonal(d) * Q
 end
 
-function _monomials_half_newton_polytope(_monos, filter)
-    monos = MP.monomial_vector(_monos)
-    basis = MB.FullBasis{MB.Monomial,eltype(monos)}()
-    return SumOfSquares.Certificate._half_newton_polytope(
-        MB.algebra_element(
-            SA.SparseCoefficients(monos, ones(length(monos))),
-            basis,
-        ),
-        MP.variables(monos),
-        filter,
-    ).monomials
-end
-
 function _randsos(
     X::AbstractVector{<:MP.AbstractMonomial};
     r = -1,
     monotype = :Classic,
     eps = 0.1,
 )
     if monotype == :Classic
-        x = _monomials_half_newton_polytope(
+        x = Certificate.monomials_half_newton_polytope(
             X,
             Certificate.NewtonDegreeBounds(tuple()),
         )

src/sosdec.jl

@@ -121,10 +121,12 @@ function MP.polynomial(d::SOSDecomposition)
     return MP.polynomial(MB.algebra_element(d))
 end
 
-function MB.algebra_element(decomp::SOSDecomposition)
-    res = zero(first(decomp.ps))
+function MB.algebra_element(decomp::SOSDecomposition{A,T,V,U}) where {A,T,V,U}
+    basis = MB.implicit_basis(SA.basis(first(decomp.ps)))
+    res = zero(U, MB.algebra(basis))
     for p in decomp.ps
-        MA.operate!(SA.UnsafeAddMul(*), res, SA.star(p), p)
+        implicit = MB.algebra_element(SA.coeffs(p, basis), basis)
+        MA.operate!(SA.UnsafeAddMul(*), res, SA.star(implicit), implicit)
     end
     MA.operate!(SA.canonical, SA.coeffs(res))
     return res

test/Project.toml

@@ -16,4 +16,4 @@ SymbolicWedderburn = "858aa9a9-4c7c-4c62-b466-2421203962a2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
-DynamicPolynomials = "0.5"
+DynamicPolynomials = "0.5,0.6"