Updates to the new StarAlgebra

.github/workflows/ci.yml

@@ -32,21 +32,11 @@ jobs:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
       - uses: julia-actions/cache@v2
-      - name: MOI
-        shell: julia --project=@. {0}
-        run: |
-          using Pkg
-          Pkg.add([
-              PackageSpec(name="StarAlgebras", rev="main"),
-              PackageSpec(name="SymbolicWedderburn", rev="master"),
-              PackageSpec(name="MultivariateBases", rev="master"),
-              PackageSpec(name="MultivariateMoments", rev="master"),
-              PackageSpec(name="PolyJuMP", rev="master"),
-          ])
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1
-        with:
-          depwarn: error
+# # See https://github.com/oxfordcontrol/Clarabel.jl/pull/230
+#        with:
+#          depwarn: error
       - uses: julia-actions/julia-processcoverage@v1
       - uses: codecov/codecov-action@v4
         with:

.github/workflows/documentation.yml

@@ -22,11 +22,7 @@ jobs:
         run: |
           using Pkg
           Pkg.add([
-              PackageSpec(name="StarAlgebras", rev="main"),
-              PackageSpec(name="SymbolicWedderburn", rev="master"),
-              PackageSpec(name="MultivariateBases", rev="master"),
-              PackageSpec(name="MultivariateMoments", rev="master"),
-              PackageSpec(name="PolyJuMP", rev="master"),
+              PackageSpec(name="MultivariateBases", rev="bl/comparable"),
               PackageSpec(path=pwd()),
           ])
           Pkg.instantiate()

.github/workflows/examples.yml

@@ -17,11 +17,7 @@ jobs:
         run: |
           using Pkg
           Pkg.add([
-              PackageSpec(name="StarAlgebras", rev="main"),
-              PackageSpec(name="SymbolicWedderburn", rev="master"),
-              PackageSpec(name="MultivariateBases", rev="master"),
-              PackageSpec(name="MultivariateMoments", rev="master"),
-              PackageSpec(name="PolyJuMP", rev="master"),
+              PackageSpec(name="MultivariateBases", rev="bl/comparable"),
               PackageSpec(path=pwd()),
           ])
           Pkg.instantiate()

Project.toml

@@ -1,7 +1,7 @@
 name = "SumOfSquares"
 uuid = "4b9e565b-77fc-50a5-a571-1244f986bda1"
-repo = "https://github.com/jump-dev/SumOfSquares.jl.git"
 version = "0.7.3"
+repo = "https://github.com/jump-dev/SumOfSquares.jl.git"
 
 [deps]
 CliqueTrees = "60701a23-6482-424a-84db-faee86b9b1f8"
@@ -22,16 +22,16 @@ SymbolicWedderburn = "858aa9a9-4c7c-4c62-b466-2421203962a2"
 
 [compat]
 CliqueTrees = "1"
-DataStructures = "0.18"
+DataStructures = "0.19"
 JuMP = "1.10"
 MathOptInterface = "1.13"
-MultivariateBases = "0.2"
-MultivariateMoments = "0.4"
+MultivariateBases = "0.3"
+MultivariateMoments = "0.5"
 MultivariatePolynomials = "0.5"
 MutableArithmetics = "1"
-PolyJuMP = "0.7"
-Reexport = "0.2, 1.0"
+PolyJuMP = "0.8"
+Reexport = "1"
 SemialgebraicSets = "0.3"
 StarAlgebras = "0.3"
-SymbolicWedderburn = "0.4"
+SymbolicWedderburn = "0.5"
 julia = "1.10"

docs/Project.toml

@@ -5,7 +5,7 @@ Clarabel = "61c947e1-3e6d-4ee4-985a-eec8c727bd6e"
 ColorSchemes = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
 Cyclotomics = "da8f5974-afbb-4dc8-91d8-516d5257c83b"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
-DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244"
 Dualization = "191a621a-6537-11e9-281d-650236a99e60"

docs/src/tutorials/Extension/certificate.jl

@@ -20,18 +20,18 @@ S = @set x >= 0 && y >= 0 && x^2 + y^2 >= 2
 # We will now see how to find the optimal solution using Sum of Squares Programming.
 # We first need to pick an SDP solver, see [here](https://jump.dev/JuMP.jl/v1.12/installation/#Supported-solvers) for a list of the available choices.
 # Note that SumOfSquares generates a *standard form* SDP (i.e., SDP variables
-# and equality constraints) while SCS expects a *geometric form* SDP (i.e.,
+# and equality constraints) while Clarabel expects a *geometric form* SDP (i.e.,
 # free variables and symmetric matrices depending affinely on these variables
 # constrained to belong to the PSD cone).
 # JuMP will transform the standard from to the geometric form will create the PSD
 # variables as free variables and then constrain then to be PSD.
 # While this will work, since the dual of a standard from is in in geometric form,
 # dualizing the problem will generate a smaller SDP.
-# We use therefore `Dualization.dual_optimizer` so that SCS solves the dual problem.
+# We use therefore `Dualization.dual_optimizer` so that Clarabel solves the dual problem.
 
-import SCS
+import Clarabel
 using Dualization
-solver = dual_optimizer(SCS.Optimizer)
+solver = dual_optimizer(Clarabel.Optimizer)
 
 # A Sum-of-Squares certificate that $p \ge \alpha$ over the domain `S`, ensures that $\alpha$ is a lower bound to the polynomial optimization problem.
 # The following program searches for the largest lower bound.
@@ -98,7 +98,7 @@ SOS.matrix_cone_type(::Type{<:Schmüdgen{IC, CT}}) where {IC, CT} = SOS.matrix_c
 model = SOSModel(solver)
 @variable(model, α)
 @objective(model, Max, α)
-basis = MB.FullBasis{MB.Monomial,typeof(x * y)}()
+basis = MB.FullBasis{MB.Monomial}(x * y)
 ideal_certificate = SOSC.Newton(SOSCone(), basis, basis, tuple())
 certificate = Schmüdgen(ideal_certificate, SOSCone(), basis, maxdegree(p))
 @constraint(model, c, p >= α, domain = S, certificate = certificate)

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

@@ -51,7 +51,7 @@ gram = gram_matrix(cref)
 
 #-
 
-gram.basis.monomials' * gram.Q * gram.basis.monomials
+keys_as_monomials(gram.basis)' * gram.Q * keys_as_monomials(gram.basis)
 
 # where the matrix `gram.Q` is positive semidefinite, because `p` is SOS. If we
 # could only get the decomposition `gram.Q = V' * V`, the SOS decomposition would

docs/src/tutorials/Getting started/sum-of-squares_matrices.jl

@@ -65,13 +65,13 @@ p = vec(y)' * P * vec(y)
 
 X = monomials(p)
 unipartite = Certificate.NewtonDegreeBounds(tuple())
-@test Certificate.monomials_half_newton_polytope(X, unipartite) == [x * y[1], x * y[2], y[1] * y[2], x, y[1], y[2]] #src
+@test Certificate.monomials_half_newton_polytope(X, unipartite) == [y[2], y[1], x, y[1] * y[2], x * y[2], x * y[1]] #src
 Certificate.monomials_half_newton_polytope(X, unipartite) #!jl
 
 # Exploiting the multipartite structure gives 4 monomials.
 
 multipartite = Certificate.NewtonDegreeBounds(([x], y))
-@test Certificate.monomials_half_newton_polytope(X, multipartite) == [x * y[1], x * y[2], y[1], y[2]] #src
+@test Certificate.monomials_half_newton_polytope(X, multipartite) == [y[2], y[1], x * y[2], x * y[1]] #src
 Certificate.monomials_half_newton_polytope(X, multipartite) #!jl
 
 # In the example above, there were only 3 monomials, where does the difference come from ?
@@ -81,12 +81,12 @@ Certificate.monomials_half_newton_polytope(X, multipartite) #!jl
 # hence the whole column and row will be zero as well.
 # Therefore, we can remove this monomial.
 
-@test Certificate.monomials_half_newton_polytope(X, Certificate.NewtonFilter(multipartite)) == [x * y[1], x * y[2], y[1]] #src
+@test Certificate.monomials_half_newton_polytope(X, Certificate.NewtonFilter(multipartite)) == [y[1], x * y[2], x * y[1]] #src
 Certificate.monomials_half_newton_polytope(X, Certificate.NewtonFilter(multipartite)) #!jl
 
 # The same reasoning can be used for monomials `y[1]y[2]` and `x` therefore whether
 # we exploit the multipartite structure or not, we get only 3 monomials thanks
 # to this post filter.
 
-@test Certificate.monomials_half_newton_polytope(X, Certificate.NewtonFilter(unipartite)) == [x * y[1], x * y[2], y[1]] #src
+@test Certificate.monomials_half_newton_polytope(X, Certificate.NewtonFilter(unipartite)) == [y[1], x * y[2], x * y[1]] #src
 Certificate.monomials_half_newton_polytope(X, Certificate.NewtonFilter(unipartite)) #!jl

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

@@ -53,12 +53,12 @@ minimizers = [η.atoms[1].center; η.atoms[2].center]
 # Below are more details on what we mean by convex combination.
 # The moment matrix of the atomic measure at the first minimizer is:
 
-η1 = moment_matrix(dirac(monomials(x, 0:4), x => round(minimizers[1])), ν.basis.monomials)
+η1 = moment_matrix(dirac(monomials(x, 0:4), x => round(minimizers[1])), keys_as_monomials(ν.basis))
 η1.Q
 
 # The moment matrix of the atomic measure at the second minimizer is:
 
-η2 = moment_matrix(dirac(monomials(x, 0:4), x => round(minimizers[2])), ν.basis.monomials)
+η2 = moment_matrix(dirac(monomials(x, 0:4), x => round(minimizers[2])), keys_as_monomials(ν.basis))
 η2.Q
 
 # And the moment matrix is the convex combination of both:

docs/src/tutorials/Polynomial Optimization/polynomial_optimization.jl

@@ -155,7 +155,7 @@ optimize!(model)
 
 # We can see that the basis of the moment matrix didn't increase:
 
-@test length(moment_matrix(c4).basis.monomials) == 3 #src
+@test length(moment_matrix(c4).basis) == 3 #src
 moment_matrix(c4)
 
 # This is because of the Newton polytope reduction that determined that gram matrix will

docs/src/tutorials/Sparsity/sign_symmetry.jl

@@ -31,14 +31,14 @@ function sos_check(sparsity)
 end
 
 g = sos_check(Sparsity.NoPattern())
-@test g.basis.monomials == [x[1]^2, x[1] * x[2], x[2]^2, x[1], x[2], x[3], 1] #src
-g.basis.monomials
+@test keys_as_monomials(g.basis) == [x[1]^2, x[1] * x[2], x[2]^2, x[1], x[2], x[3], 1] #src
+g.basis
 
 # As detailed in the Example 4 of [L09], we can exploit the *sign symmetry* of
 # the polynomial to decompose the large positive semidefinite matrix into smaller ones.
 
 g = sos_check(Sparsity.SignSymmetry())
-monos = [sub.basis.monomials for sub in g.blocks]
+monos = [keys_as_monomials(sub.basis) for sub in g.blocks]
 @test length(monos) == 3 #src
 @test [x[1], x[2]] in monos #src
 @test [x[3]] in monos #src

docs/src/tutorials/Sparsity/term_sparsity.jl

@@ -45,7 +45,7 @@ atomic_measure(ν, 1e-6)
 
 # We can see below that the basis contained 6 monomials hence we needed to use 6x6 PSD matrix variables.
 
-@test ν.basis.monomials == [x[1]*x[2], x[2]*x[3], x[1], x[2], x[3], 1] #src
+@test keys_as_monomials(ν.basis) == [x[1]*x[2], x[2]*x[3], x[1], x[2], x[3], 1] #src
 ν.basis
 
 # Using the monomial/term sparsity method of [WML20a] based on cluster completion, we find the same bound.
@@ -57,7 +57,7 @@ bound
 # Which is not suprising as no sparsity reduction could be performed.
 
 @test length(ν.blocks) == 1 #src
-@test ν.blocks[1].basis.monomials == [x[1]*x[2], x[2]*x[3], x[1], x[2], x[3], 1] #src
+@test keys_as_monomials(ν.blocks[1].basis) == [x[1]*x[2], x[2]*x[3], x[1], x[2], x[3], 1] #src
 [sub.basis for sub in ν.blocks]
 
 # Using the monomial/term sparsity method of [WML20b] based on chordal completion, the lower bound is smaller than 0.
@@ -69,8 +69,8 @@ bound
 # However, this bound was obtained with an SDP with 4 matrices of size 3x3.
 
 @test length(ν.blocks) == 4                                       #src
-@test ν.blocks[1].basis.monomials == [1, x[2]*x[3], x[1]*x[2]]    #src
-@test ν.blocks[2].basis.monomials == [x[3], x[2], x[2]*x[3]]      #src
-@test ν.blocks[3].basis.monomials == [x[2], x[2]*x[3], x[1]*x[2]] #src
-@test ν.blocks[4].basis.monomials == [x[2], x[1], x[1]*x[2]]      #src
+@test keys_as_monomials(ν.blocks[1].basis) == [1, x[2]*x[3], x[1]*x[2]]    #src
+@test keys_as_monomials(ν.blocks[2].basis) == [x[3], x[2], x[2]*x[3]]      #src
+@test keys_as_monomials(ν.blocks[3].basis) == [x[2], x[2]*x[3], x[1]*x[2]] #src
+@test keys_as_monomials(ν.blocks[4].basis) == [x[2], x[1], x[1]*x[2]]      #src
 [sub.basis for sub in ν.blocks]

docs/src/tutorials/Symmetry/cyclic.jl

@@ -101,16 +101,12 @@ solution_summary(model)
 gram = gram_matrix(con_ref).blocks #src
 @test length(gram) == 2                       #src
 @test gram[1].Q ≈ [0 0; 0 2] #src
-polys = gram[1].basis.bases[].elements #src
-@test length(polys) == 2 #src
-@test polys[1] ≈ 1 #src
-@test polys[2] ≈ -sum(x)/√3 #src
+@test gram[1].basis[1].elements[] ≈ 1 #src
+@test gram[1].basis[2].elements[] ≈ -sum(x)/√3 #src
 @test gram[2].Q ≈ [0.5;;] #src
-@test length(gram[2].basis.bases) == 2 #src
-polys = gram[2].basis.bases[1].elements #src
-@test polys[] ≈ (x[1] + x[2] - 2x[3])/√6 #src
-polys = gram[2].basis.bases[2].elements #src
-@test polys[] ≈ (x[1] - x[2])/√2 #src
+@test length(gram[2].basis[1].elements) == 2 #src
+@test gram[2].basis[1].elements[1] ≈ (x[1] + x[2] - 2x[3])/√6 #src
+@test gram[2].basis[1].elements[2] ≈ (x[1] - x[2])/√2 #src
 gram_matrix(con_ref)
 
 # Let's look into more details at the last two elements of the basis.
@@ -150,18 +146,12 @@ solution_summary(model)
 gram = gram_matrix(con_ref).blocks #src
 @test length(gram) == 3                       #src
 @test gram[1].Q ≈ [0 0; 0 2] #src
-polys = gram[1].basis.bases[].elements #src
-@test length(polys) == 2 #src
-@test polys[1] ≈ 1 #src
-@test polys[2] ≈ -sum(x)/√3 #src
+@test gram[1].basis[1].elements[] ≈ 1 #src
+@test gram[1].basis[2].elements[] ≈ -sum(x)/√3 #src
 @test gram[2].Q ≈ [0.5;;] rtol = 1e-6 #src
-polys = gram[2].basis.bases[].elements #src
-@test length(polys) == 1 #src
-@test polys[] ≈ (basis[1] + basis[2] * im) / √2  #src
+@test gram[2].basis[1].elements[] ≈ (basis[1] + basis[2] * im) / √2  #src
 @test gram[3].Q ≈ [0.5;;] rtol = 1e-6 #src
-polys = gram[3].basis.bases[].elements #src
-@test length(polys) == 1 #src
-@test polys[] ≈ (basis[1] - basis[2] * im) / √2 #src
+@test gram[3].basis[1].elements[] ≈ (basis[1] - basis[2] * im) / √2 #src
 gram_matrix(con_ref)
 
 # We can see that the real invariant subspace was in fact coming from two complex conjugate complex invariant subspaces:

docs/src/tutorials/Symmetry/dihedral.jl

@@ -148,17 +148,13 @@ function solve(G)
 
     g = gram_matrix(con_ref).blocks #src
     @test length(g) == 4            #src
-    @test length(g[4].basis.bases) == 2 #src
-    polys = g[4].basis.bases[1].elements #src
-    @test length(polys) == 3 #src
-    @test polys[1] ≈ y^3 #src
-    @test polys[2] ≈ x^2*y #src
-    @test polys[3] ≈ y #src
-    polys = g[4].basis.bases[2].elements #src
-    @test length(polys) == 3 #src
-    @test polys[1] ≈ -x^3 #src
-    @test polys[2] ≈ -x*y^2 #src
-    @test polys[3] ≈ -x #src
+    @test length(g[4].basis[1].elements) == 2 #src
+    @test g[4].basis[1].elements[1] ≈ y^3 #src
+    @test g[4].basis[2].elements[1] ≈ x^2*y #src
+    @test g[4].basis[3].elements[1] ≈ y #src
+    @test g[4].basis[1].elements[2] ≈ -x^3 #src
+    @test g[4].basis[2].elements[2] ≈ -x*y^2 #src
+    @test g[4].basis[3].elements[2] ≈ -x #src
     I = 3:-1:1                      #src
     Q = g[4].Q[I, I]                #src
     @test size(Q) == (3, 3)         #src
@@ -168,20 +164,16 @@ function solve(G)
     @test Q[1, 1] ≈ 25/64 rtol=1e-2 #src
     @test Q[1, 3] ≈ -5/8  rtol=1e-2 #src
     @test Q[3, 3] ≈  1    rtol=1e-2 #src
-    polys = g[1].basis.bases[].elements #src
-    @test length(polys) == 2 #src
-    @test polys[1] ≈ 1.0 #src
-    @test polys[2] ≈ -(√2/2)x^2 - (√2/2)y^2 #src
+    @test g[1].basis[1].elements[] ≈ 1.0 #src
+    @test g[1].basis[2].elements[] ≈ -(√2/2)x^2 - (√2/2)y^2 #src
     @test size(g[1].Q) == (2, 2)             #src
     @test g[1].Q[1, 1] ≈ 7921/4096 rtol=1e-2 #src
     @test g[1].Q[1, 2] ≈ 0.983 rtol=1e-2  #src
     @test g[1].Q[2, 2] ≈ 1/2 rtol=1e-2 #src
-    polys = g[2].basis.bases[].elements #src
-    @test polys[] ≈ x * y #src
+    @test g[2].basis[1].elements[] ≈ x * y #src
     @test size(g[2].Q) == (1, 1)       #src
     @test g[2].Q[1, 1] ≈ 0   atol=1e-2 #src
-    polys = g[3].basis.bases[].elements #src
-    @test polys[] ≈ (√2/2)x^2 - (√2/2)y^2 #src
+    @test g[3].basis[1].elements[] ≈ (√2/2)x^2 - (√2/2)y^2 #src
     @test size(g[3].Q) == (1, 1)       #src
     @test g[3].Q[1, 1] ≈ 0   atol=1e-2 #src
     gram_matrix(con_ref)

docs/src/tutorials/Symmetry/even_reduction.jl

@@ -45,9 +45,7 @@ value(t)
 # We indeed find `-1`, let's verify that symmetry was exploited:
 
 @test length(gram_matrix(con_ref).blocks) == 2 #src
-polys = gram_matrix(con_ref).blocks[1].basis.bases[].elements #src
-@test polys[1] ≈ 1 #src
-@test polys[2] ≈ x^2 #src
-polys = gram_matrix(con_ref).blocks[2].basis.bases[].elements #src
-@test polys[] ≈ x #src
+@test gram_matrix(con_ref).blocks[1].basis[1].elements[] ≈ 1 #src
+@test gram_matrix(con_ref).blocks[1].basis[2].elements[] ≈ x^2 #src
+@test gram_matrix(con_ref).blocks[2].basis[1].elements[] ≈ x #src
 gram_matrix(con_ref)

docs/src/tutorials/Symmetry/permutation_symmetry.jl

@@ -47,26 +47,18 @@ value(t)
 
 gram = gram_matrix(con_ref).blocks    #src
 @test length(gram) == 3                          #src
-polys = gram[1].basis.bases[].elements #src
-@test length(polys) == 2     #src
-@test polys[1] ≈ 1        #src
-@test polys[2] ≈ -0.5 * sum(x) #src
+@test gram[1].basis[1].elements[] ≈ 1        #src
+@test gram[1].basis[2].elements[] ≈ -0.5 * sum(x) #src
 @test size(gram[1].Q) == (2, 2)                  #src
 @test gram[1].Q[1, 1] ≈ 1.0 atol=1e-6            #src
 @test gram[1].Q[1, 2] ≈ -1.0 atol=1e-6            #src
 @test gram[1].Q[2, 2] ≈ 1.0 atol=1e-6            #src
-@test length(gram[2].basis.bases) == 2 #src
-polys = gram[2].basis.bases[1].elements #src
-@test length(polys) == 1     #src
-@test polys[] ≈ (x[2] - x[4]) / √2 #src
-@test size(g[2].Q) == (1, 1)                  #src
-@test g[2].Q[1, 1] ≈ 1.0 atol=1e-6            #src
-polys = gram[2].basis.bases[2].elements #src
-@test length(polys) == 1     #src
-@test polys[1] ≈ (x[1] - x[3]) / √2 #src
-polys = gram[3].basis.bases[].elements #src
-@test length(polys) == 1     #src
-@test polys[] ≈ (x[1] - x[2] + x[3] - x[4]) / 2 #src
+@test length(gram[2].basis[1].elements) == 2 #src
+@test gram[2].basis[1].elements[1] ≈ (x[2] - x[4]) / √2 #src
+@test size(gram[2].Q) == (1, 1)                  #src
+@test gram[2].Q[1, 1] ≈ 1.0 atol=1e-6            #src
+@test gram[2].basis[1].elements[2] ≈ (x[1] - x[3]) / √2 #src
+@test gram[3].basis[1].elements[] ≈ (x[1] - x[2] + x[3] - x[4]) / 2 #src
 @test size(gram[3].Q) == (1, 1)                  #src
 @test gram[3].Q[1, 1] ≈ 1.0 atol=1e-6            #src
 gram_matrix(con_ref)