Use WeightedSOSCone by default

.github/workflows/ci.yml

@@ -27,7 +27,7 @@ jobs:
 #            arch: x64
     steps:
       - uses: actions/checkout@v4
-      - uses: julia-actions/setup-julia@v1
+      - uses: julia-actions/setup-julia@v2
         with:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
@@ -37,8 +37,12 @@ jobs:
         run: |
           using Pkg
           Pkg.add([
+              PackageSpec(name="StarAlgebras", rev="main"),
+              PackageSpec(name="SymbolicWedderburn", rev="mk/StarAlgebras_0_3"),
               PackageSpec(name="MultivariateBases", rev="master"),
-              PackageSpec(name="MathOptInterface", rev="master"),
+              PackageSpec(name="SemialgebraicSets", rev="master"),
+              PackageSpec(name="MultivariateMoments", rev="master"),
+              PackageSpec(name="PolyJuMP", rev="master"),
           ])
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1

.github/workflows/documentation.yml

@@ -22,6 +22,13 @@ jobs:
         run: |
           using Pkg
           Pkg.add([
+              PackageSpec(url="https://github.com/blegat/HomotopyContinuation.jl", rev="dynamicpolynomials_v0.6"),
+              PackageSpec(name="StarAlgebras", rev="main"),
+              PackageSpec(name="SymbolicWedderburn", rev="mk/StarAlgebras_0_3"),
+              PackageSpec(name="MultivariateBases", rev="master"),
+              PackageSpec(name="SemialgebraicSets", rev="master"),
+              PackageSpec(name="MultivariateMoments", rev="master"),
+              PackageSpec(name="PolyJuMP", rev="master"),
               PackageSpec(path=pwd()),
           ])
           Pkg.instantiate()

.github/workflows/examples.yml

@@ -17,6 +17,13 @@ jobs:
         run: |
           using Pkg
           Pkg.add([
+              PackageSpec(url="https://github.com/blegat/HomotopyContinuation.jl", rev="dynamicpolynomials_v0.6"),
+              PackageSpec(name="StarAlgebras", rev="main"),
+              PackageSpec(name="SymbolicWedderburn", rev="mk/StarAlgebras_0_3"),
+              PackageSpec(name="MultivariateBases", rev="master"),
+              PackageSpec(name="SemialgebraicSets", rev="master"),
+              PackageSpec(name="MultivariateMoments", rev="master"),
+              PackageSpec(name="PolyJuMP", rev="master"),
               PackageSpec(path=pwd()),
           ])
           Pkg.instantiate()

Project.toml

@@ -16,6 +16,7 @@ PolyJuMP = "ddf597a6-d67e-5340-b84c-e37d84115374"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SemialgebraicSets = "8e049039-38e8-557d-ae3a-bc521ccf6204"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+StarAlgebras = "0c0c59c1-dc5f-42e9-9a8b-b5dc384a6cd1"
 SymbolicWedderburn = "858aa9a9-4c7c-4c62-b466-2421203962a2"
 
 [compat]
@@ -29,5 +30,5 @@ MutableArithmetics = "1"
 PolyJuMP = "0.7"
 Reexport = "0.2, 1.0"
 SemialgebraicSets = "0.3"
-SymbolicWedderburn = "0.3"
+SymbolicWedderburn = "0.4"
 julia = "1.6"

bench/sos_polynomial.jl

@@ -5,7 +5,7 @@ function sos_polynomial(n)
     @polyvar x
     monos = monomials(x, 0:2n)
     set = SumOfSquares.SOSPolynomialSet(
-        SumOfSquares.FullSpace(), monos,
+        SumOfSquares.FullSpace(), MB.SubBasis{MB.Monomial}(monos),
         SumOfSquares.Certificate.Remainder(
             SumOfSquares.SOSCone(),
             SumOfSquares.MonomialBasis,

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/make.jl

@@ -13,7 +13,7 @@ const _TUTORIAL_SUBDIR = [
     "Other Applications",
     "Noncommutative and Hermitian",
     "Sparsity",
-    "Symmetry",
+    #"Symmetry",
     "Extension",
 ]
 

docs/src/constraints.md

@@ -41,13 +41,13 @@ CachingOptimizer state: NO_OPTIMIZER
 Solver name: No optimizer attached.
 
 julia> @variable(model, p, Poly(X))
-(_[1]) + (_[2])x₃ + (_[3])x₂ + (_[4])x₁ + (_[5])x₃² + (_[6])x₂x₃ + (_[7])x₂² + (_[8])x₁x₃ + (_[9])x₁x₂ + (_[10])x₁²
+(_[1])·1 + (_[2])·x₃ + (_[3])·x₂ + (_[4])·x₁ + (_[5])·x₃² + (_[6])·x₂x₃ + (_[7])·x₂² + (_[8])·x₁x₃ + (_[9])·x₁x₂ + (_[10])·x₁²
 
 julia> @variable(model, q, Poly(X))
-(_[11]) + (_[12])x₃ + (_[13])x₂ + (_[14])x₁ + (_[15])x₃² + (_[16])x₂x₃ + (_[17])x₂² + (_[18])x₁x₃ + (_[19])x₁x₂ + (_[20])x₁²
+(_[11])·1 + (_[12])·x₃ + (_[13])·x₂ + (_[14])·x₁ + (_[15])·x₃² + (_[16])·x₂x₃ + (_[17])·x₂² + (_[18])·x₁x₃ + (_[19])·x₁x₂ + (_[20])·x₁²
 
 julia> @constraint(model, p + q == 1)
-(_[1] + _[11] - 1) + (_[2] + _[12])x₃ + (_[3] + _[13])x₂ + (_[4] + _[14])x₁ + (_[5] + _[15])x₃² + (_[6] + _[16])x₂x₃ + (_[7] + _[17])x₂² + (_[8] + _[18])x₁x₃ + (_[9] + _[19])x₁x₂ + (_[10] + _[20])x₁² ∈ PolyJuMP.ZeroPoly()
+(_[1] + _[11] - 1)·1 + (_[2] + _[12])·x₃ + (_[3] + _[13])·x₂ + (_[4] + _[14])·x₁ + (_[5] + _[15])·x₃² + (_[6] + _[16])·x₂x₃ + (_[7] + _[17])·x₂² + (_[8] + _[18])·x₁x₃ + (_[9] + _[19])·x₁x₂ + (_[10] + _[20])·x₁² ∈ PolyJuMP.ZeroPoly()
 ```
 
 Vectorized constraints can also be used as well as vector of constraints,
@@ -66,7 +66,7 @@ Polynomials can be constrained to be sum-of-squares with the `in` syntax.
 For instance, to constrain a polynomial `p` to be sum-of-squares, do
 ```jldoctest constraint-pq
 julia> @constraint(model, p in SOSCone())
-(_[1]) + (_[2])x₃ + (_[3])x₂ + (_[4])x₁ + (_[5])x₃² + (_[6])x₂x₃ + (_[7])x₂² + (_[8])x₁x₃ + (_[9])x₁x₂ + (_[10])x₁² is SOS
+(_[1])·1 + (_[2])·x₃ + (_[3])·x₂ + (_[4])·x₁ + (_[5])·x₃² + (_[6])·x₂x₃ + (_[7])·x₂² + (_[8])·x₁x₃ + (_[9])·x₁x₂ + (_[10])·x₁² is SOS
 ```
 
 ### Automatically interpreting polynomial nonnegativity as a sum-of-squares constraint
@@ -147,13 +147,13 @@ julia> @variable(model, β)
 β
 
 julia> @constraint(model, α * x^2 + β * y^2 ≥ (α - β) * x * y)
-(β)y² + (-α + β)xy + (α)x² is SOS
+(β)·y² + (-α + β)·xy + (α)·x² is SOS
 ```
 where `α` and `β` are JuMP decision variables and `x` and `y` are polynomial
 variables. Since the polynomial is a quadratic form, the sum-of-squares
 certificate is also a quadratic form (see [Blekherman2012; Section~3.3.4](@cite)). Hence the
 default polynomial basis used for the [Nonnegative polynomial variables]
-certificate is `MonomialBasis([x, y])`, that is, we search for a positive
+certificate is `MultivariateBases.SubBasis{MultivariateBases.Monomial}([x, y])`, that is, we search for a positive
 semidefinite matrix `Q` such that
 ```math
 \alpha x^2 + \beta y^2 - (\alpha - \beta) x y = X^\top Q X
@@ -164,16 +164,16 @@ As the polynomial space is determined by the polynomial being constrained,
 only the basis *type* needs to be given. For instance, to use the scaled monomial
 basis in the example above, use
 ```jldoctest constraint-xy
-julia> @constraint(model, α * x^2 + β * y^2 ≥ (α - β) * x * y, basis = ScaledMonomialBasis)
-(β)y² + (-α + β)xy + (α)x² is SOS
+julia> @constraint(model, α * x^2 + β * y^2 ≥ (α - β) * x * y, basis = ScaledMonomial)
+(β)·ScaledMonomial(y²) + (-0.7071067811865475 α + 0.7071067811865475 β)·ScaledMonomial(xy) + (α)·ScaledMonomial(x²) is SOS
 ```
 
 ## Polynomial nonnegativity on a subset of the space
 
 By default, the constraint
 ```jldoctest constraint-xy
 julia> @constraint(model, x^3 - x^2 + 2x*y -y^2 + y^3 >= α)
-(-α) + (-1)y² + (2)xy + (-1)x² + (1)y³ + (1)x³ is SOS
+(-α)·1 + (-1)·y² + (2)·xy + (-1)·x² + (1)·y³ + (1)·x³ is SOS
 ```
 constrains the polynomial to be nonnegative for every real numbers `x` and `y`.
 However, the set of points `(x, y)` for which the polynomial is constrained
@@ -182,7 +182,7 @@ to be nonnegative can be specified by the `domain` keyword:
 julia> S = @set x >= 0 && y >= 0 && x + y >= 1;
 
 julia> @constraint(model, x^3 - x^2 + 2x*y -y^2 + y^3 >= α, domain = S)
-(-α) + (-1)y² + (2)xy + (-1)x² + (1)y³ + (1)x³ is SOS
+(-α)·1 + (-1)·y² + (2)·xy + (-1)·x² + (1)·y³ + (1)·x³ is SOS
 ```
 See [this notebook](https://github.com/jump-dev/SumOfSquares.jl/blob/master/examples/Polynomial_Optimization.ipynb)
 for a detailed example.

docs/src/reference/constraints.md

@@ -20,6 +20,7 @@ SOSDecomposition
 SOSDecompositionWithDomain
 SumOfSquares.SOSDecompositionAttribute
 sos_decomposition
+SumOfSquares.MultiplierIndexBoundsError
 ```
 
 SAGE decomposition attribute:

docs/src/tutorials/Extension/certificate.jl

@@ -50,13 +50,14 @@ 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
-struct Schmüdgen{IC <: SOSC.AbstractIdealCertificate, CT <: SOS.SOSLikeCone, BT <: SOS.AbstractPolynomialBasis} <: SOSC.AbstractPreorderCertificate
+struct Schmüdgen{IC<:SOSC.AbstractIdealCertificate,CT<:SOS.SOSLikeCone,B<:SA.AbstractBasis} <: SOSC.AbstractPreorderCertificate
     ideal_certificate::IC
     cone::CT
-    basis::Type{BT}
+    basis::B
     maxdegree::Int
 end
 
@@ -78,8 +79,8 @@ function SOSC.multiplier_basis(certificate::Schmüdgen, index::SOSC.PreorderInde
     q = SOSC.generator(certificate, index, domain)
     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
+function SOSC.multiplier_basis_type(::Type{Schmüdgen{IC, CT, BT}}, ::Type{M}) where {IC,CT,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/motzkin.jl

@@ -90,7 +90,7 @@ value(deno)
 # We can easily extend `Plots` by adding a recipe to plot bivariate polynomials.
 
 using RecipesBase
-@recipe function f(x::AbstractVector, y::AbstractVector, p::Polynomial)
+@recipe function f(x::AbstractVector, y::AbstractVector, p::AbstractPolynomial)
     x, y, (x, y) -> p(variables(p) => [x, y])
 end
 import Plots

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)

docs/src/tutorials/Systems and Control/lyapunov_function_search.jl

@@ -63,4 +63,4 @@ JuMP.primal_status(model)
 # We can now obtain this feasible solution with:
 
 value(V)
-@test iszero(remove_monomials(value(V), monos)) #src
+@test iszero(remove_monomials(polynomial(value(V)), monos)) #src