diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 3ba2c1dc6..ce02b523d 100644
--- a/.github/workflows/ci.yml
+++ b/.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:
diff --git a/.github/workflows/documentation.yml b/.github/workflows/documentation.yml
index e64f7e5eb..6eb44fdee 100644
--- a/.github/workflows/documentation.yml
+++ b/.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()
diff --git a/.github/workflows/examples.yml b/.github/workflows/examples.yml
index 9b9ebab5e..d002108a9 100644
--- a/.github/workflows/examples.yml
+++ b/.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()
diff --git a/Project.toml b/Project.toml
index 2128c7566..83754b12e 100644
--- a/Project.toml
+++ b/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"
diff --git a/docs/Project.toml b/docs/Project.toml
index 6d8379355..2abdd3d17 100644
--- a/docs/Project.toml
+++ b/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"
diff --git a/docs/src/tutorials/Extension/certificate.jl b/docs/src/tutorials/Extension/certificate.jl
index 86fe5a980..38e1861c7 100644
--- a/docs/src/tutorials/Extension/certificate.jl
+++ b/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)
diff --git a/docs/src/tutorials/Getting started/sos_decomposition.jl b/docs/src/tutorials/Getting started/sos_decomposition.jl
index c7bb4e949..001f9f50c 100644
--- a/docs/src/tutorials/Getting started/sos_decomposition.jl	
+++ b/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
diff --git a/docs/src/tutorials/Getting started/sum-of-squares_matrices.jl b/docs/src/tutorials/Getting started/sum-of-squares_matrices.jl
index 265b24605..534987c9c 100644
--- a/docs/src/tutorials/Getting started/sum-of-squares_matrices.jl	
+++ b/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
diff --git a/docs/src/tutorials/Getting started/univariate.jl b/docs/src/tutorials/Getting started/univariate.jl
index e08fe64cb..e98876c0f 100644
--- a/docs/src/tutorials/Getting started/univariate.jl	
+++ b/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:
diff --git a/docs/src/tutorials/Polynomial Optimization/polynomial_optimization.jl b/docs/src/tutorials/Polynomial Optimization/polynomial_optimization.jl
index 4cb550869..c24fcf8cd 100644
--- a/docs/src/tutorials/Polynomial Optimization/polynomial_optimization.jl	
+++ b/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
diff --git a/docs/src/tutorials/Sparsity/sign_symmetry.jl b/docs/src/tutorials/Sparsity/sign_symmetry.jl
index 03275a318..c0f99e02b 100644
--- a/docs/src/tutorials/Sparsity/sign_symmetry.jl
+++ b/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
diff --git a/docs/src/tutorials/Sparsity/term_sparsity.jl b/docs/src/tutorials/Sparsity/term_sparsity.jl
index dbe5a6e19..815efee28 100644
--- a/docs/src/tutorials/Sparsity/term_sparsity.jl
+++ b/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]
diff --git a/docs/src/tutorials/Symmetry/cyclic.jl b/docs/src/tutorials/Symmetry/cyclic.jl
index 27f687e96..2d4a20948 100644
--- a/docs/src/tutorials/Symmetry/cyclic.jl
+++ b/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:
diff --git a/docs/src/tutorials/Symmetry/dihedral.jl b/docs/src/tutorials/Symmetry/dihedral.jl
index 94e85867e..946c2140d 100644
--- a/docs/src/tutorials/Symmetry/dihedral.jl
+++ b/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)
diff --git a/docs/src/tutorials/Symmetry/even_reduction.jl b/docs/src/tutorials/Symmetry/even_reduction.jl
index 856878967..8d3a5433e 100644
--- a/docs/src/tutorials/Symmetry/even_reduction.jl
+++ b/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)
diff --git a/docs/src/tutorials/Symmetry/permutation_symmetry.jl b/docs/src/tutorials/Symmetry/permutation_symmetry.jl
index f2dd41f9e..36e9513c7 100644
--- a/docs/src/tutorials/Symmetry/permutation_symmetry.jl
+++ b/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)
diff --git a/docs/src/tutorials/Systems and Control/barrier_certificate.jl b/docs/src/tutorials/Systems and Control/barrier_certificate.jl
index 4f174da9a..74f1ed19a 100644
--- a/docs/src/tutorials/Systems and Control/barrier_certificate.jl	
+++ b/docs/src/tutorials/Systems and Control/barrier_certificate.jl	
@@ -52,8 +52,8 @@ JuMP.primal_status(model)
 @test JuMP.primal_status(model) == MOI.FEASIBLE_POINT #src
 
 # Plot the phase plot with the 0-level set of the barrier function, and the boundary of the initial and unsafe sets
-import DifferentialEquations, Plots, ImplicitPlots
-function phase_plot(f, B, g₁, h₁, quiver_scaling, Δt, X0, solver = DifferentialEquations.Tsit5())
+import OrdinaryDiffEq, Plots, ImplicitPlots
+function phase_plot(f, B, g₁, h₁, quiver_scaling, Δt, X0, solver = OrdinaryDiffEq.Tsit5())
     X₀plot = ImplicitPlots.implicit_plot(h₁; xlims=(-2, 3), ylims=(-2.5, 2.5), resolution = 1000, label="X₀", linecolor=:blue)
     Xᵤplot = ImplicitPlots.implicit_plot!(g₁; xlims=(-2, 3), ylims=(-2.5, 2.5), resolution = 1000, label="Xᵤ", linecolor=:teal)
     Bplot  = ImplicitPlots.implicit_plot!(B; xlims=(-2, 3), ylims=(-2.5, 2.5), resolution = 1000, label="B = 0", linecolor=:red)
@@ -64,8 +64,8 @@ function phase_plot(f, B, g₁, h₁, quiver_scaling, Δt, X0, solver = Differen
     ∇pt(v, p, t) = ∇(v[1], v[2])
     function traj(v0)
         tspan = (0.0, Δt)
-        prob = DifferentialEquations.ODEProblem(∇pt, v0, tspan)
-        return DifferentialEquations.solve(prob, solver, reltol=1e-8, abstol=1e-8)
+        prob = OrdinaryDiffEq.ODEProblem(∇pt, v0, tspan)
+        return OrdinaryDiffEq.solve(prob, solver, reltol=1e-8, abstol=1e-8)
     end
     ticks = -5:0.5:5
     X = repeat(ticks, 1, length(ticks))
diff --git a/docs/src/tutorials/Systems and Control/julia_set.jl b/docs/src/tutorials/Systems and Control/julia_set.jl
index a78cf9ed7..29daa90ee 100644
--- a/docs/src/tutorials/Systems and Control/julia_set.jl	
+++ b/docs/src/tutorials/Systems and Control/julia_set.jl	
@@ -108,8 +108,8 @@ end
 
 # We 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.
 
-import CSDP
-solver = optimizer_with_attributes(CSDP.Optimizer, MOI.Silent() => true)
+import Clarabel
+solver = optimizer_with_attributes(Clarabel.Optimizer, MOI.Silent() => true)
 
 # Let's start with the value of `c` corresponding to the left image of [KHJ14, Figure 3] and with degree 2.
 
diff --git a/docs/src/tutorials/Systems and Control/stabilization_of_nonlinear_systems.jl b/docs/src/tutorials/Systems and Control/stabilization_of_nonlinear_systems.jl
index fdf948912..4e91f71f4 100644
--- a/docs/src/tutorials/Systems and Control/stabilization_of_nonlinear_systems.jl	
+++ b/docs/src/tutorials/Systems and Control/stabilization_of_nonlinear_systems.jl	
@@ -36,14 +36,14 @@ function controller(f, g, b, α, degs)
     return MultivariatePolynomials.map_coefficients(coef -> abs(coef) < 1e-6 ? 0.0 : coef, u)
 end
 
-import DifferentialEquations, Plots
-function phase_plot(f, quiver_scaling, Δt, X0, solver = DifferentialEquations.Tsit5())
+import OrdinaryDiffEq, Plots
+function phase_plot(f, quiver_scaling, Δt, X0, solver = OrdinaryDiffEq.Tsit5())
     ∇(vx, vy) = [fi(x[1] => vx, x[2] => vy) for fi in f]
     ∇pt(v, p, t) = ∇(v[1], v[2])
     function traj(v0)
         tspan = (0.0, Δt)
-        prob = DifferentialEquations.ODEProblem(∇pt, v0, tspan)
-        return DifferentialEquations.solve(prob, solver, reltol=1e-8, abstol=1e-8)
+        prob = OrdinaryDiffEq.ODEProblem(∇pt, v0, tspan)
+        return OrdinaryDiffEq.solve(prob, solver, reltol=1e-8, abstol=1e-8)
     end
     ticks = -5:0.5:5
     X = repeat(ticks, 1, length(ticks))
diff --git a/docs/src/variables.md b/docs/src/variables.md
index 5d17b4698..277a0ddfc 100644
--- a/docs/src/variables.md
+++ b/docs/src/variables.md
@@ -43,23 +43,23 @@ product between `a` and `X`.
 
 Just like with classical JuMP's decision variables, containers of polynomial
 variables can be created as follows:
-```jldoctest variables
+```jldoctest variables; filter = [r"(Matrix|Vector|DenseAxisArray|SparseAxisArray)\{.*\}" => s"\1{…}"]
 julia> @variable(model, [1:3, 1:4], Poly(X))       # Creates a Matrix
-3×4 Matrix{StarAlgebras.AlgebraElement{MultivariateBases.Algebra{SubBasis{Monomial, DynamicPolynomials.Monomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, Graded{LexOrder}}, MonomialVector{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, Graded{LexOrder}}}, Monomial, DynamicPolynomials.Monomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, Graded{LexOrder}}}, VariableRef, Vector{VariableRef}}}:
+3×4 Matrix{…}:
  (_[7])·1 + (_[8])·y + (_[9])·x + (_[10])·y² + (_[11])·xy + (_[12])·x²     …  (_[61])·1 + (_[62])·y + (_[63])·x + (_[64])·y² + (_[65])·xy + (_[66])·x²
  (_[13])·1 + (_[14])·y + (_[15])·x + (_[16])·y² + (_[17])·xy + (_[18])·x²     (_[67])·1 + (_[68])·y + (_[69])·x + (_[70])·y² + (_[71])·xy + (_[72])·x²
  (_[19])·1 + (_[20])·y + (_[21])·x + (_[22])·y² + (_[23])·xy + (_[24])·x²     (_[73])·1 + (_[74])·y + (_[75])·x + (_[76])·y² + (_[77])·xy + (_[78])·x²
 
 julia> @variable(model, [[:a, :b], -2:2], Poly(X)) # Creates a DenseAxisArray
-2-dimensional DenseAxisArray{StarAlgebras.AlgebraElement{MultivariateBases.Algebra{SubBasis{Monomial, DynamicPolynomials.Monomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}, DynamicPolynomials.MonomialVector{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}}, Monomial, DynamicPolynomials.Monomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}}, VariableRef, Vector{VariableRef}},2,...} with index sets:
+2-dimensional DenseAxisArray{…} with index sets:
     Dimension 1, [:a, :b]
     Dimension 2, -2:2
-And data, a 2×5 Matrix{StarAlgebras.AlgebraElement{MultivariateBases.Algebra{SubBasis{Monomial, DynamicPolynomials.Monomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}, DynamicPolynomials.MonomialVector{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}}, Monomial, DynamicPolynomials.Monomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}}, VariableRef, Vector{VariableRef}}}:
+And data, a 2×5 Matrix{…}:
  (_[79])·1 + (_[80])·y + (_[81])·x + (_[82])·y² + (_[83])·xy + (_[84])·x²  …  (_[127])·1 + (_[128])·y + (_[129])·x + (_[130])·y² + (_[131])·xy + (_[132])·x²
  (_[85])·1 + (_[86])·y + (_[87])·x + (_[88])·y² + (_[89])·xy + (_[90])·x²     (_[133])·1 + (_[134])·y + (_[135])·x + (_[136])·y² + (_[137])·xy + (_[138])·x²
 
 julia> @variable(model, [i=1:3, j=i:3], Poly(X))   # Creates a SparseAxisArray
-JuMP.Containers.SparseAxisArray{StarAlgebras.AlgebraElement{MultivariateBases.Algebra{SubBasis{Monomial, DynamicPolynomials.Monomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, Graded{LexOrder}}, MonomialVector{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, Graded{LexOrder}}}, Monomial, DynamicPolynomials.Monomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, Graded{LexOrder}}}, VariableRef, Vector{VariableRef}}, 2, Tuple{Int64, Int64}} with 6 entries:
+JuMP.Containers.SparseAxisArray{…} with 6 entries:
   [1, 1]  =  (_[139])·1 + (_[140])·y + (_[141])·x + (_[142])·y^2 + (_[143])·x*y + (_[144])·x^2
   [1, 2]  =  (_[145])·1 + (_[146])·y + (_[147])·x + (_[148])·y^2 + (_[149])·x*y + (_[150])·x^2
   [1, 3]  =  (_[151])·1 + (_[152])·y + (_[153])·x + (_[154])·y^2 + (_[155])·x*y + (_[156])·x^2
@@ -94,9 +94,9 @@ In order to create a sum-of-squares polynomial variable, the syntax is exactly
 the same except `SOSPoly` should be used instead of `Poly`.
 For instance, the following code creates a ``3 \times 4`` matrix of
 sum-of-squares polynomial variables:
-```jldoctest variables
+```jldoctest variables; filter = [r"(Matrix|Vector|DenseAxisArray|SparseAxisArray)\{.*\}" => s"\1{…}"]
 julia> @variable(model, [1:2], SOSPoly(X))
-2-element Vector{GramMatrix{VariableRef, SubBasis{Monomial, DynamicPolynomials.Monomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, Graded{LexOrder}}, MonomialVector{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, Graded{LexOrder}}}, AffExpr, SymMatrix{VariableRef}}}:
+2-element Vector{…}:
  GramMatrix with row/column basis:
  SubBasis{Monomial}([1, y, x, y^2, x*y, x^2])
 And entries in a 6×6 SymMatrix{VariableRef}:
diff --git a/examples/Project.toml b/examples/Project.toml
index 41cf159a4..55a44372b 100644
--- a/examples/Project.toml
+++ b/examples/Project.toml
@@ -2,7 +2,7 @@
 CSDP = "0a46da34-8e4b-519e-b418-48813639ff34"
 Cyclotomics = "da8f5974-afbb-4dc8-91d8-516d5257c83b"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
-DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 DynamicPolynomials = "7c1d4256-1411-5781-91ec-d7bc3513ac07"
 GroupsCore = "d5909c97-4eac-4ecc-a3dc-fdd0858a4120"
 HomotopyContinuation = "f213a82b-91d6-5c5d-acf7-10f1c761b327"
diff --git a/examples/run_examples.jl b/examples/run_examples.jl
index 6afbc1d4b..efe7c5152 100644
--- a/examples/run_examples.jl
+++ b/examples/run_examples.jl
@@ -29,9 +29,7 @@ end
 
 @testset "run_examples.jl" begin
     @testset "$dir" for dir in readdir(_TUTORIAL_DIR)
-        if dir != "Symmetry"
-            run_examples(dir)
-        end
+        run_examples(dir)
     end
     @testset "Chordal" begin
         include("chordal_sparsity.jl")
diff --git a/src/Bridges/Constraint/image.jl b/src/Bridges/Constraint/image.jl
index 41dd9e87e..2ecdce6fe 100644
--- a/src/Bridges/Constraint/image.jl
+++ b/src/Bridges/Constraint/image.jl
@@ -81,17 +81,17 @@ function MOI.Bridges.Constraint.bridge_constraint(
     @assert MOI.output_dimension(g) == length(set.basis)
     scalars = MOI.Utilities.scalarize(g)
     k = 0
-    found = Dict{eltype(set.basis.monomials),Int}()
+    found = Dict{MP.monomial_type(set.basis),Int}()
     first = Union{Nothing,Int}[nothing for _ in eachindex(scalars)]
     variables = MOI.VariableIndex[]
     constraints = MOI.ConstraintIndex{F}[]
     for (gram_basis, weight) in zip(set.gram_bases, set.weights)
         cone = SOS.matrix_cone(M, length(gram_basis))
         f = MOI.Utilities.zero_with_output_dimension(F, MOI.dimension(cone))
-        for j in eachindex(gram_basis.monomials)
+        for j in eachindex(gram_basis)
             for i in 1:j
                 k += 1
-                mono = gram_basis.monomials[i] * gram_basis.monomials[j]
+                mono = MP.monomial(gram_basis[i]) * MP.monomial(gram_basis[j])
                 is_diag = i == j
                 if haskey(found, mono)
                     var = MOI.add_variable(model)
@@ -119,7 +119,7 @@ function MOI.Bridges.Constraint.bridge_constraint(
                     MOI.Utilities.operate_output_index!(-, T, k, f, var)
                 else
                     found[mono] = k
-                    t = MB.monomial_index(set.basis, mono)
+                    t = SA.key_index(set.basis, MP.exponents(mono))
                     if !isnothing(t)
                         first[t] = k
                         if is_diag
@@ -305,8 +305,8 @@ function MOI.get(
 ) where {T}
     dual =
         MOI.get(model, MOI.ConstraintDual(attr.result_index), bridge.constraint)
-    output = similar(dual, length(bridge.set.monomials))
-    for i in eachindex(bridge.set.monomials)
+    output = similar(dual, length(bridge.set.basis))
+    for i in eachindex(bridge.set.basis)
         output[i] = dual[bridge.first[i]]
     end
     return output
diff --git a/src/Bridges/Constraint/sos_polynomial.jl b/src/Bridges/Constraint/sos_polynomial.jl
index c3910c755..79b90cc36 100644
--- a/src/Bridges/Constraint/sos_polynomial.jl
+++ b/src/Bridges/Constraint/sos_polynomial.jl
@@ -42,8 +42,10 @@ end
 
 function _poly(coeffs, basis::MB.MonomialIndexedBasis{B,M}) where {B,M}
     return MB.algebra_element(
-        MB.sparse_coefficients(MP.polynomial(coeffs, basis.monomials)),
-        MB.FullBasis{B,M}(),
+        MB.sparse_coefficients(
+            MP.polynomial(coeffs, MB.keys_as_monomials(basis)),
+        ),
+        MB.implicit_basis(basis),
     )
 end
 
@@ -65,7 +67,7 @@ function MOI.Bridges.Constraint.bridge_constraint(
         # MOI does not modify the coefficients of the functions so we can modify `p`.
         # without altering `f`.
         # The basis may be copied by MA however so we need to copy it.
-        _poly(MOI.Utilities.scalarize(func), copy(set.basis)),
+        _poly(MOI.Utilities.scalarize(func), set.basis), # TODO use `copy(set.basis)` once https://github.com/JuliaAlgebra/StarAlgebras.jl/pull/83 is done
         domain,
     )
     gram_basis = SOS.Certificate.gram_basis(
@@ -73,7 +75,7 @@ function MOI.Bridges.Constraint.bridge_constraint(
         SOS.Certificate.with_variables(poly, set.domain),
     )
     gram_bases = [gram_basis]
-    weights = [MB.constant_algebra_element(typeof(SA.basis(poly)), T)]
+    weights = [MB.constant_algebra_element(SA.basis(poly), T)]
     flat_gram_bases, flat_weights, flat_indices = _flatten(gram_bases, weights)
     new_basis = SOS.Certificate.zero_basis(
         set.certificate,
@@ -115,7 +117,10 @@ function MOI.Bridges.Constraint.concrete_bridge_type(
     # promotes VectorOfVariables into VectorAffineFunction, it should be enough
     # for most use cases
     M = SOS.matrix_cone_type(CT)
-    W = SOS.Certificate._weight_type(T, BT)
+    W = MB.constant_algebra_element_type(
+        MA.promote_operation(MB.implicit_basis, BT),
+        T,
+    )
     B = MA.promote_operation(
         SOS.Certificate.zero_basis,
         CT,
diff --git a/src/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl b/src/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl
index de3bbe59b..79e8efb98 100644
--- a/src/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl
+++ b/src/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl
@@ -22,29 +22,23 @@ struct SOSPolynomialInSemialgebraicSetBridge{
     },
     UMST,
     MCT,
-    MT<:MP.AbstractMonomial,
-    MVT<:AbstractVector{MT},
+    BT<:MB.SubBasis{MB.Monomial},
 } <: MOI.Bridges.Constraint.AbstractBridge
     lagrangian_bases::Vector{B}
     lagrangian_variables::Vector{
         Union{Vector{MOI.VariableIndex},Vector{Vector{MOI.VariableIndex}}},
     }
     lagrangian_constraints::Vector{UMCT}
-    constraint::MOI.ConstraintIndex{
-        F,
-        SOS.SOSPolynomialSet{DT,MB.SubBasis{MB.Monomial,MT,MVT},CT},
-    }
-    monomials::MVT
+    constraint::MOI.ConstraintIndex{F,SOS.SOSPolynomialSet{DT,BT,CT}}
+    basis::BT
 end
 
 function MOI.Bridges.Constraint.bridge_constraint(
-    ::Type{
-        SOSPolynomialInSemialgebraicSetBridge{T,F,DT,CT,B,UMCT,UMST,MCT,MT,MVT},
-    },
+    ::Type{SOSPolynomialInSemialgebraicSetBridge{T,F,DT,CT,B,UMCT,UMST,MCT,BT}},
     model::MOI.ModelLike,
     f::MOI.AbstractVectorFunction,
     set::SOS.SOSPolynomialSet{<:SemialgebraicSets.BasicSemialgebraicSet},
-) where {T,F,DT,CT,B,UMCT,UMST,MCT,MT,MVT}
+) where {T,F,DT,CT,B,UMCT,UMST,MCT,BT}
     @assert MOI.output_dimension(f) == length(set.basis)
     # MOI does not modify the coefficients of the functions so we can modify `p`.
     # without altering `f`.
@@ -52,8 +46,9 @@ function MOI.Bridges.Constraint.bridge_constraint(
     # TODO remove `collect` when `DynamicPolynomials.MonomialVector` can be used as keys
     p = MB.algebra_element(
         SA.SparseCoefficients(
-            copy(collect(set.basis.monomials)),
+            copy(collect(set.basis.keys)),
             MOI.Utilities.scalarize(f),
+            SA.comparable(MB.implicit_basis(set.basis)),
         ),
         MB.implicit_basis(set.basis),
     )
@@ -96,10 +91,7 @@ function MOI.Bridges.Constraint.bridge_constraint(
     MA.operate!(SA.canonical, SA.coeffs(p))
     new_set = SOS.SOSPolynomialSet(
         set.domain.V,
-        # For terms, `monomials` is `OneOrZeroElementVector`
-        # so we convert it with `monomial_vector`
-        # Later, we'll use `MP.MonomialBasis` which is going to do that anyway
-        MB.SubBasis{MB.Monomial}(MP.monomial_vector(SA.keys(SA.coeffs(p)))),
+        MB.explicit_basis(p),
         Certificate.ideal_certificate(set.certificate),
     )
     constraint = MOI.add_constraint(
@@ -108,23 +100,12 @@ function MOI.Bridges.Constraint.bridge_constraint(
         new_set,
     )
 
-    return SOSPolynomialInSemialgebraicSetBridge{
-        T,
-        F,
-        DT,
-        CT,
-        B,
-        UMCT,
-        UMST,
-        MCT,
-        MT,
-        MVT,
-    }(
+    return SOSPolynomialInSemialgebraicSetBridge{T,F,DT,CT,B,UMCT,UMST,MCT,BT}(
         λ_bases,
         λ_variables,
         λ_constraints,
         constraint,
-        set.basis.monomials,
+        set.basis,
     )
 end
 
@@ -141,11 +122,9 @@ function MOI.Bridges.added_constrained_variable_types(
     return constrained_variable_types(SOS.matrix_cone_type(CT))
 end
 function MOI.Bridges.added_constraint_types(
-    ::Type{
-        SOSPolynomialInSemialgebraicSetBridge{T,F,DT,CT,B,UMCT,UMST,MCT,MT,MVT},
-    },
-) where {T,F,DT,CT,B,UMCT,UMST,MCT,MT,MVT}
-    return [(F, SOS.SOSPolynomialSet{DT,MB.SubBasis{MB.Monomial,MT,MVT},CT})]
+    ::Type{SOSPolynomialInSemialgebraicSetBridge{T,F,DT,CT,B,UMCT,UMST,MCT,BT}},
+) where {T,F,DT,CT,B,UMCT,UMST,MCT,BT}
+    return [(F, SOS.SOSPolynomialSet{DT,BT,CT})]
 end
 function MOI.Bridges.Constraint.concrete_bridge_type(
     ::Type{<:SOSPolynomialInSemialgebraicSetBridge{T}},
@@ -153,32 +132,22 @@ function MOI.Bridges.Constraint.concrete_bridge_type(
     ::Type{
         <:SOS.SOSPolynomialSet{
             SemialgebraicSets.BasicSemialgebraicSet{S,PS,AT},
-            MB.SubBasis{MB.Monomial,MT,MVT},
+            BT,
             CT,
         },
     },
-) where {T,S,PS,AT,CT,MT,MVT}
+) where {T,S,PS,AT,CT,BT<:MB.MonomialIndexedBasis{MB.Monomial}}
 
     # promotes VectorOfVariables into VectorAffineFunction, it should be enough
     # for most use cases
     G = MOI.Utilities.promote_operation(-, T, F, MOI.VectorOfVariables)
     MCT = SOS.matrix_cone_type(CT)
+    MT = MP.monomial_type(BT)
     B = Certificate.multiplier_basis_type(CT, MT)
     UMCT = union_constraint_types(MCT)
     UMST = union_set_types(MCT)
     IC = Certificate.ideal_certificate(CT)
-    return SOSPolynomialInSemialgebraicSetBridge{
-        T,
-        G,
-        AT,
-        IC,
-        B,
-        UMCT,
-        UMST,
-        MCT,
-        MT,
-        MVT,
-    }
+    return SOSPolynomialInSemialgebraicSetBridge{T,G,AT,IC,B,UMCT,UMST,MCT,BT}
 end
 
 # Attributes, Bridge acting as an model
@@ -209,7 +178,7 @@ function MOI.get(
     )
 end
 function MOI.get(
-    b::SOSPolynomialInSemialgebraicSetBridge{T,F,DT,CT,B,UMCT,UMST},
+    ::SOSPolynomialInSemialgebraicSetBridge{T,F,DT,CT,B,UMCT,UMST},
     ::MOI.ListOfConstraintIndices{MOI.VectorOfVariables,S},
 ) where {T,F,DT,CT,B,UMCT,UMST,S<:UMST}
     C = MOI.ConstraintIndex{MOI.VectorOfVariables,S}
@@ -218,21 +187,15 @@ function MOI.get(
     ])
 end
 function MOI.get(
-    ::SOSPolynomialInSemialgebraicSetBridge{T,F,DT,CT,B,UMCT,UMST,MCT,MT,MVT},
-    ::MOI.NumberOfConstraints{
-        F,
-        SOS.SOSPolynomialSet{DT,MB.SubBasis{MB.Monomial,MT,MVT},CT},
-    },
-) where {T,F,DT,CT,B,UMCT,UMST,MCT,MT,MVT}
+    ::SOSPolynomialInSemialgebraicSetBridge{T,F,DT,CT,B,UMCT,UMST,MCT,BT},
+    ::MOI.NumberOfConstraints{F,SOS.SOSPolynomialSet{DT,BT,CT}},
+) where {T,F,DT,CT,B,UMCT,UMST,MCT,BT}
     return 1
 end
 function MOI.get(
-    b::SOSPolynomialInSemialgebraicSetBridge{T,F,DT,CT,B,UMCT,UMST,MCT,MT,MVT},
-    ::MOI.ListOfConstraintIndices{
-        F,
-        SOS.SOSPolynomialSet{DT,MB.SubBasis{MB.Monomial,MT,MVT},CT},
-    },
-) where {T,F,DT,CT,B,UMCT,UMST,MCT,MT,MVT}
+    b::SOSPolynomialInSemialgebraicSetBridge{T,F,DT,CT,B,UMCT,UMST,MCT,BT},
+    ::MOI.ListOfConstraintIndices{F,SOS.SOSPolynomialSet{DT,BT,CT}},
+) where {T,F,DT,CT,B,UMCT,UMST,MCT,BT}
     return [b.constraint]
 end
 
@@ -280,7 +243,8 @@ function MOI.get(
     dual = MOI.get(model, attr, bridge.constraint)
     set = MOI.get(model, MOI.ConstraintSet(), bridge.constraint)
     μ = MultivariateMoments.moment_vector(dual, set.basis)
-    return [dot(mono, μ) for mono in bridge.monomials]
+    monos = MB.keys_as_monomials(bridge.basis)
+    return [dot(mono, μ) for mono in monos]
 end
 function MOI.get(
     model::MOI.ModelLike,
diff --git a/src/Certificate/Certificate.jl b/src/Certificate/Certificate.jl
index e75a279b9..a450e93ff 100644
--- a/src/Certificate/Certificate.jl
+++ b/src/Certificate/Certificate.jl
@@ -54,7 +54,7 @@ _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(_vec(MP.variables(a))))
-    comm = is_commutative(vars)
+    comm = MP.is_commutative(vars)
     return all(vars) do v
         return _degree(a, v, comm) <= _degree(b, v, comm)
     end
@@ -74,34 +74,52 @@ function within_bounds(mono, bounds)
 end
 
 function maxdegree_gram_basis(
-    ::MB.FullBasis{B},
+    full::MB.FullBasis{B},
     bounds::DegreeBounds,
 ) where {B<:MB.AbstractMonomial}
     variables = MP.variables(bounds.variablewise_maxdegree)
-    return MB.SubBasis{B}(
-        MP.monomials(
-            variables,
-            bounds.mindegree:bounds.maxdegree,
-            Base.Fix2(within_variablewise_bounds, bounds),
-        ),
+    monos = MP.monomials(
+        variables,
+        bounds.mindegree:bounds.maxdegree,
+        Base.Fix2(within_variablewise_bounds, bounds),
     )
+    sub = MB.SubBasis{B}(monos)
+    new_sub, new_full = SA.promote_bases(sub, full)
+    @assert new_full === full
+    return new_sub
 end
 
 function maxdegree_gram_basis(basis::SA.AbstractBasis, ::Nothing)
-    return MB.empty_basis(MB.explicit_basis_type(typeof(basis)))
+    return SA.SubBasis(basis, SA.key_type(basis)[])
 end
 
 function maxdegree_gram_basis(basis::SA.AbstractBasis, bounds::DegreeBounds)
     # TODO use bounds here too
-    variables = MP.variables(bounds.variablewise_maxdegree)
-    return MB.maxdegree_basis(basis, variables, bounds.maxdegree)
+    @assert MP.variables(basis) == MP.variables(bounds.variablewise_maxdegree)
+    return MB.maxdegree_basis(basis, bounds.maxdegree)
 end
+
 function maxdegree_gram_basis(
     basis::SA.AbstractBasis,
     variables,
     maxdegree::Int,
 )
-    return MB.maxdegree_basis(basis, variables, fld(maxdegree, 2))
+    if MP.variables(basis) == variables
+        return MB.maxdegree_basis(basis, fld(maxdegree, 2))
+    else
+        return _maxdegree_gram_basis(basis, variables, fld(maxdegree, 2))
+    end
+end
+
+function _maxdegree_gram_basis(
+    full::MB.FullBasis{B},
+    variables,
+    halfdegree::Int,
+) where {B}
+    monos = MP.monomials(variables, 0:halfdegree)
+    sub = MB.SubBasis{B}(monos)
+    new_sub, _ = SA.promote_bases(sub, full)
+    return new_sub
 end
 
 include("ideal.jl")
diff --git a/src/Certificate/Sparsity/ChordalExtensionGraph.jl b/src/Certificate/Sparsity/ChordalExtensionGraph.jl
index 3f0a1d31a..f0a1fe691 100644
--- a/src/Certificate/Sparsity/ChordalExtensionGraph.jl
+++ b/src/Certificate/Sparsity/ChordalExtensionGraph.jl
@@ -187,7 +187,7 @@ Return a cluster completion of `G` and the corresponding maximal cliques.
 """
 function completion(G::Graph, ::ClusterCompletion)
     H = copy(G)
-    union_find = DataStructures.IntDisjointSets(num_nodes(G))
+    union_find = DataStructures.IntDisjointSet(num_nodes(G))
     for from in 1:num_nodes(G)
         for to in neighbors(G, from)
             DataStructures.union!(union_find, from, to)
diff --git a/src/Certificate/Sparsity/ideal.jl b/src/Certificate/Sparsity/ideal.jl
index e3d8500a2..071e87b36 100644
--- a/src/Certificate/Sparsity/ideal.jl
+++ b/src/Certificate/Sparsity/ideal.jl
@@ -63,6 +63,21 @@ function Ideal(
     )
 end
 
+# `clique has a subset of variables`
+function _maxdegree_basis(
+    ambient_basis::MB.FullBasis{B},
+    clique,
+    maxdegree,
+) where {B}
+    basis = SumOfSquares.Certificate.maxdegree_gram_basis(
+        MB.FullBasis{B}(clique),
+        clique,
+        maxdegree,
+    )
+    full_basis, _ = SA.promote_bases(basis, ambient_basis)
+    return full_basis
+end
+
 function sparsity(
     poly,
     ::Variable,
@@ -71,7 +86,7 @@ function sparsity(
     basis = MB.explicit_basis(poly)
     H, cliques = chordal_csp_graph(basis, SemialgebraicSets.FullSpace())
     return map(cliques) do clique
-        return SumOfSquares.Certificate.maxdegree_gram_basis(
+        return _maxdegree_basis(
             certificate.gram_basis,
             clique,
             certificate.maxdegree,
@@ -83,7 +98,9 @@ function sparsity(
     sp::Union{SignSymmetry,Monomial},
     gram_basis::MB.SubBasis{MB.Monomial},
 )
-    return MB.SubBasis{MB.Monomial}.(sparsity(monos, sp, gram_basis.monomials))
+    return MB.SubBasis{MB.Monomial}.(
+        sparsity(monos, sp, MB.keys_as_monomials(gram_basis)),
+    )
 end
 function sparsity(
     poly,
@@ -91,7 +108,7 @@ function sparsity(
     certificate::SumOfSquares.Certificate.AbstractIdealCertificate,
 )
     return sparsity(
-        MB.explicit_basis(poly).monomials,
+        MB.keys_as_monomials(MB.explicit_basis(poly)),
         sp,
         SumOfSquares.Certificate.gram_basis(certificate, poly),
     )
diff --git a/src/Certificate/Sparsity/monomial.jl b/src/Certificate/Sparsity/monomial.jl
index 598fa72fd..64a9a9493 100644
--- a/src/Certificate/Sparsity/monomial.jl
+++ b/src/Certificate/Sparsity/monomial.jl
@@ -183,7 +183,7 @@ function sparsity(
     return _monomial_vector(cliques)
 end
 # This also checks that it is indeed a monomial basis
-_monos(basis::MB.SubBasis{MB.Monomial}) = basis.monomials
+_monos(basis::MB.SubBasis{MB.Monomial}) = MB.keys_as_monomials(basis)
 function _gram_monos(vars, certificate::SumOfSquares.Certificate.MaxDegree)
     return _monos(
         SumOfSquares.Certificate.maxdegree_gram_basis(
@@ -215,8 +215,7 @@ struct DummyPolynomial{M}
 end
 function SumOfSquares.Certificate._algebra_element(p::DummyPolynomial)
     return MB.algebra_element(
-        SA.SparseCoefficients(p.monomials, ones(length(p.monomials))),
-        MB.FullBasis{MB.Monomial,eltype(p.monomials)}(),
+        MP.polynomial(ones(length(p.monomials)), p.monomials),
     )
 end
 MP.monomials(p::DummyPolynomial) = p.monomials
@@ -253,14 +252,14 @@ function sparsity(
             SumOfSquares.Certificate.ideal_certificate(certificate),
             DummyPolynomial(
                 _ideal_monos(
-                    MB.explicit_basis(poly).monomials,
+                    MB.keys_as_monomials(MB.explicit_basis(poly)),
                     multiplier_generator_monos,
                 ),
             ),
         ),
     )
     cliques, multiplier_cliques = sparsity(
-        MB.explicit_basis(poly).monomials,
+        MB.keys_as_monomials(MB.explicit_basis(poly)),
         sp,
         gram_monos,
         multiplier_generator_monos,
diff --git a/src/Certificate/Sparsity/variable.jl b/src/Certificate/Sparsity/variable.jl
index f6b81a0d3..580461fe2 100644
--- a/src/Certificate/Sparsity/variable.jl
+++ b/src/Certificate/Sparsity/variable.jl
@@ -10,8 +10,8 @@ struct Variable <: Pattern end
 const CEG = ChordalExtensionGraph
 
 function csp_graph(basis::MB.SubBasis{MB.Monomial}, ::FullSpace)
-    G = CEG.LabelledGraph{MP.variable_union_type(eltype(basis.monomials))}()
-    for mono in basis.monomials
+    G = CEG.LabelledGraph{MP.variable_union_type(MP.monomial_type(basis))}()
+    for mono in MB.keys_as_monomials(basis)
         CEG.add_clique!(G, MP.effective_variables(mono))
     end
     return G
@@ -54,7 +54,7 @@ function sparsity(
 )
     basis = MB.explicit_basis(poly)
     H, cliques = chordal_csp_graph(basis, domain)
-    function bases(q)
+    function clique_bases(q, vars)
         return [
             SumOfSquares.Certificate.maxdegree_gram_basis(
                 certificate.multipliers_certificate.gram_basis,
@@ -63,8 +63,11 @@ function sparsity(
                     certificate.maxdegree,
                     q,
                 ),
-            ) for clique in cliques if MP.variables(q) ⊆ clique
+            ) for clique in cliques if vars ⊆ clique
         ]
     end
-    return bases(basis), map(bases, domain.p)
+    ideal_bases = clique_bases(basis, MP.variables(basis))
+    preorder_bases =
+        map(q -> clique_bases(q, MP.effective_variables(q)), domain.p)
+    return ideal_bases, preorder_bases
 end
diff --git a/src/Certificate/Symmetry/block_diag.jl b/src/Certificate/Symmetry/block_diag.jl
index 3a2a04f7b..ab1ad2fbf 100644
--- a/src/Certificate/Symmetry/block_diag.jl
+++ b/src/Certificate/Symmetry/block_diag.jl
@@ -253,7 +253,7 @@ function block_diag(As, d)
         iZ = Z'
         @assert iZ ≈ inv(Z)
         n = LinearAlgebra.checksquare(A)
-        union_find = DataStructures.IntDisjointSets(n)
+        union_find = DataStructures.IntDisjointSet(n)
         Bs = [iZ * A * Z for A in As]
         for B in Bs
             merge_sparsity!(union_find, B)
@@ -288,7 +288,7 @@ function block_diag(As, d)
 end
 
 function merge_sparsity!(
-    union_find::DataStructures.IntDisjointSets,
+    union_find::DataStructures.IntDisjointSet,
     A,
     tol = 1e-8,
 )
diff --git a/src/Certificate/Symmetry/wedderburn.jl b/src/Certificate/Symmetry/wedderburn.jl
index bb50e8b66..4934695bc 100644
--- a/src/Certificate/Symmetry/wedderburn.jl
+++ b/src/Certificate/Symmetry/wedderburn.jl
@@ -47,7 +47,7 @@ function SymbolicWedderburn.action(
     el,
     p::MB.Polynomial{MB.Monomial},
 )
-    res = SymbolicWedderburn.action(a, el, p.monomial)
+    res = SymbolicWedderburn.action(a, el, MP.monomial(p))
     if res isa MP.AbstractMonomial
         return MB.Polynomial{MB.Monomial}(res)
     else
@@ -95,10 +95,12 @@ function SumOfSquares.matrix_cone_type(::Type{<:Ideal{C}}) where {C}
     return SumOfSquares.matrix_cone_type(C)
 end
 
-function _multi_basis_type(::Type{<:MB.SubBasis{B,M}}, ::Type{T}) where {B,M,T}
-    SC = SA.SparseCoefficients{M,T,Vector{M},Vector{T}}
-    AE = SA.AlgebraElement{MB.Algebra{MB.FullBasis{B,M},B,M},T,SC}
-    return Vector{MB.SemisimpleBasis{AE,Int,MB.FixedBasis{B,M,T,SC}}}
+function _multi_basis_type(::Type{BT}, ::Type{T}) where {BT<:SA.SubBasis,T}
+    AE = MB.algebra_element_type(
+        Vector{T},
+        MA.promote_operation(MB.implicit_basis, BT),
+    )
+    return Vector{MB.SemisimpleBasis{AE,MB.SimpleBasis{AE}}}
 end
 function MA.promote_operation(
     ::typeof(SumOfSquares.Certificate.gram_basis),
@@ -172,13 +174,14 @@ function MA.promote_operation(
 end
 
 function matrix_reps(pattern, R, basis, ::Type{T}, form) where {T}
-    polys = R * basis.monomials
+    monos = MB.keys_as_monomials(basis)
+    polys = R * monos
     return map(SymbolicWedderburn.gens(pattern.group)) do g
         S = Matrix{T}(undef, length(polys), length(polys))
         for i in eachindex(polys)
             p = polys[i]
             q = SymbolicWedderburn.action(pattern.action, g, p)
-            coefs = MP.coefficients(q, basis.monomials)
+            coefs = MP.coefficients(q, MB.keys_as_monomials(basis))
             S[:, i] = _linsolve(R, coefs, form)
         end
         return S
@@ -195,9 +198,14 @@ function SumOfSquares.Certificate.gram_basis(cert::Ideal, poly)
 end
 
 function _fixed_basis(F, basis)
-    return MB.FixedBasis([
-        MB.implicit(MB.algebra_element(row, basis)) for row in eachrow(F)
-    ])
+    return MB.SimpleBasis(
+        map(eachrow(F)) do row
+            ae = MB.implicit(MB.algebra_element(collect(row), basis))
+            # Copy coefficients to avoid sharing the basis.keys vector,
+            # which would be mutated by `canonical` (called during hash/==)
+            return SA.AlgebraElement(copy(SA.coeffs(ae)), Base.parent(ae))
+        end,
+    )
 end
 
 function _gram_basis(pattern::Pattern, basis, ::Type{T}) where {T}
diff --git a/src/Certificate/ideal.jl b/src/Certificate/ideal.jl
index b7d9813cb..4def1276a 100644
--- a/src/Certificate/ideal.jl
+++ b/src/Certificate/ideal.jl
@@ -15,12 +15,19 @@ Base.:+(::Number, a::_NonZero) = a
 Base.:+(::_NonZero, a::_NonZero) = a
 
 function _combine_with_gram(
-    basis::MB.SubBasis{B,M},
+    basis::MB.SubBasis{B},
     gram_bases::AbstractVector{<:SA.ExplicitBasis},
     weights,
-) where {B,M}
-    p = zero(_NonZero, MB.algebra(MB.FullBasis{B,M}()))
-    cache = zero(_NonZero, MB.algebra(MB.FullBasis{B,M}()))
+) where {B}
+    for gram_basis in gram_bases
+        if parent(basis) != parent(gram_basis)
+            error(
+                "Bases $(parent(basis)) and $(parent(gram_basis)) are incompatible",
+            )
+        end
+    end
+    p = zero(_NonZero, MB.algebra(parent(basis)))
+    cache = zero(_NonZero, MB.algebra(parent(basis)))
     for mono in basis
         MA.operate!(
             SA.UnsafeAddMul(*),
@@ -38,7 +45,7 @@ function _combine_with_gram(
         MA.operate!(SA.UnsafeAddMul(*), p, cache, weight)
     end
     MA.operate!(SA.canonical, SA.coeffs(p))
-    return MB.SubBasis{B}(keys(SA.coeffs(p)))
+    return SA.sub_basis(parent(basis), keys(SA.coeffs(p)))
 end
 
 function _reduce_with_domain(basis::MB.SubBasis, zero_basis, ::FullSpace)
@@ -52,6 +59,7 @@ end
 function __reduce_with_domain(_, _, _)
     return error("Only Monomial basis support with an equalities in domain")
 end
+
 function __reduce_with_domain(
     basis::MB.SubBasis{MB.Monomial},
     ::MB.FullBasis{MB.Monomial},
@@ -59,10 +67,13 @@ function __reduce_with_domain(
 )
     I = ideal(domain)
     # set of standard monomials that are hit
-    standard = Set{eltype(basis.monomials)}()
-    for mono in basis.monomials
+    standard = Set{MP.monomial_type(basis)}()
+    for mono in MB.keys_as_monomials(basis)
         r = rem(mono, I)
-        union!(standard, MP.monomials(r))
+        # `FixedVariablesSet`, use substitutions which makes us
+        # loose variables
+        s, _ = SA.promote_bases(r, mono)
+        union!(standard, MP.monomials(s))
     end
     return MB.QuotientBasis(
         MB.SubBasis{MB.Monomial}(MP.monomial_vector(collect(standard))),
@@ -244,9 +255,12 @@ struct Remainder{C<:AbstractIdealCertificate} <: AbstractIdealCertificate
 end
 
 function _rem(coeffs, basis::MB.FullBasis{MB.Monomial}, I)
-    poly = MP.polynomial(SA.values(coeffs), SA.keys(coeffs))
+    poly = MP.polynomial(MB.algebra_element(coeffs, basis))
     r = convert(typeof(poly), rem(poly, I))
-    return MB.algebra_element(MB.sparse_coefficients(r), basis)
+    # `FixedVariablesSet`, use substitutions which makes us
+    # loose variables
+    s, _ = SA.promote_bases(r, poly)
+    return MB.algebra_element(MB.sparse_coefficients(s), basis)
 end
 
 function reduced_polynomial(::Remainder, a::SA.AlgebraElement, domain)
diff --git a/src/Certificate/newton_polytope.jl b/src/Certificate/newton_polytope.jl
index fc7d525ae..742513f0f 100644
--- a/src/Certificate/newton_polytope.jl
+++ b/src/Certificate/newton_polytope.jl
@@ -1,7 +1,3 @@
-function is_commutative(vars)
-    return length(vars) < 2 || prod(vars[1:2]) == prod(reverse(vars[1:2]))
-end
-
 abstract type AbstractNewtonPolytopeApproximation end
 
 # `maxdegree` and `mindegree` on each variable separately
@@ -82,7 +78,7 @@ function min_degree(p::SA.AlgebraElement, v)
 end
 function min_degree(p::MB.Polynomial{B}, v) where {B}
     if _is_monomial_basis(B)
-        return MP.degree(p.monomial, v)
+        return MP.degree(MP.monomial(p), v)
     else
         error("TODO $B")
     end
@@ -98,7 +94,7 @@ end
 function max_degree(p::SA.AlgebraElement, v)
     return mapreduce(Base.Fix2(max_degree, v), max, SA.supp(p))
 end
-max_degree(p::MB.Polynomial, v) = MP.degree(p.monomial, v)
+max_degree(p::MB.Polynomial, v) = MP.degree(MP.monomial(p), v)
 
 function min_shift(d, shift)
     return max(0, d - shift)
@@ -244,12 +240,12 @@ end
 #_maxdegree(basis::MB.SubBasis{MB.Monomial}) = MP.maxdegree(basis.monomials)
 function _mindegree(p::MB.Polynomial{B}, vars) where {B}
     if _is_monomial_basis(B)
-        _sum_degree(p.monomial, vars)
+        _sum_degree(MP.monomial(p), vars)
     else
         error("TODO $B")
     end
 end
-_maxdegree(p::MB.Polynomial, vars) = _sum_degree(p.monomial, vars)
+_maxdegree(p::MB.Polynomial, vars) = _sum_degree(MP.monomial(p), vars)
 function _degree(mono, var::MP.AbstractVariable, comm::Bool)
     if comm
         return MP.degree(mono, var)
@@ -264,17 +260,20 @@ function _degree(mono, var::MP.AbstractVariable, comm::Bool)
     end
 end
 function _sum_degree(mono, vars)
-    comm = is_commutative(vars)
+    comm = MP.is_commutative(vars)
     return sum(var -> _degree(mono, var, comm), vars)
 end
 
 _is_monomial_basis(::Type{<:MB.AbstractMonomialIndexed}) = false
 _is_monomial_basis(::Type{<:Union{MB.Monomial,MB.ScaledMonomial}}) = true
 function _is_monomial_basis(
-    ::Type{<:SA.AlgebraElement{<:MB.Algebra{BT,B}}},
-) where {BT,B}
+    ::Type{<:Union{MB.FullBasis{B},MB.SubBasis{B}}},
+) where {B}
     return _is_monomial_basis(B)
 end
+function _is_monomial_basis(::Type{A}) where {A<:SA.AlgebraElement}
+    return _is_monomial_basis(MA.promote_operation(SA.basis, A))
+end
 
 # Minimum degree of a gram basis for a gram matrix `s`
 # such that the minimum degree of `s * g` is at least `mindegree`.
@@ -390,19 +389,16 @@ function putinar_degree_bounds(
     )
 end
 
-function multiplier_basis(
-    g::SA.AlgebraElement{<:MB.Algebra{BT,B}},
-    bounds::DegreeBounds,
-) where {BT,B}
+function multiplier_basis(g::SA.AlgebraElement, bounds::DegreeBounds)
     return maxdegree_gram_basis(
-        MB.FullBasis{B,MP.monomial_type(typeof(g))}(),
+        MB.implicit_basis(SA.basis(g)),
         _half(minus_shift(bounds, g)),
     )
 end
 
 # Cartesian product of the newton polytopes of the different parts
 function _cartesian_product(bases::Vector{<:MB.SubBasis{B}}, bounds) where {B}
-    monos = [b.monomials for b in bases]
+    monos = [MB.keys_as_monomials(b) for b in bases]
     basis = MB.SubBasis{B}(
         vec([prod(monos) for monos in Iterators.product(monos...)]),
     )
@@ -412,19 +408,22 @@ function _cartesian_product(bases::Vector{<:MB.SubBasis{B}}, bounds) where {B}
         return MB.empty_basis(typeof(basis))
     else
         return MB.SubBasis{B}(
-            filter(Base.Fix2(within_bounds, _half(bounds)), basis.monomials),
+            filter(
+                Base.Fix2(within_bounds, _half(bounds)),
+                MB.keys_as_monomials(basis),
+            ),
         )
     end
 end
 
 function half_newton_polytope(
-    p::SA.AlgebraElement{<:MB.Algebra{BT,B,M}},
+    p::SA.AlgebraElement,
     gs::AbstractVector{<:SA.AlgebraElement},
     vars,
     maxdegree,
     newton::NewtonDegreeBounds,
-) where {BT,B,M}
-    if !is_commutative(vars)
+)
+    if !MP.is_commutative(vars)
         throw(
             ArgumentError(
                 "Multipartite Newton polytope not supported with noncommutative variables.",
@@ -481,16 +480,16 @@ function half_newton_polytope(
 end
 
 function half_newton_polytope(
-    p::SA.AlgebraElement{<:MB.Algebra{BT,B,M}},
+    p::SA.AlgebraElement,
     gs::AbstractVector{<:SA.AlgebraElement},
     vars,
     maxdegree,
     ::NewtonDegreeBounds{Tuple{}},
-) where {BT,B,M}
-    if is_commutative(vars)
+)
+    full = MB.implicit_basis(SA.basis(p))
+    if MP.is_commutative(vars)
         # TODO take `variable_groups` into account
         bounds = putinar_degree_bounds(p, gs, vars, maxdegree)
-        full = MB.FullBasis{B,M}()
         return maxdegree_gram_basis(full, _half(bounds)),
         MB.explicit_basis_type(typeof(full))[
             multiplier_basis(g, bounds) for g in gs
@@ -509,18 +508,19 @@ function half_newton_polytope(
         # Berlin: Springer, 2016.
         vars = unique!(sort(vars))
         bounds = _half(putinar_degree_bounds(p, gs, vars, maxdegree))
-        monos = MP.monomial_type(typeof(p))[]
-        for mono in MB.explicit_basis(p).monomials
+        keys = SA.key_type(full)[]
+        for mono in MB.keys_as_monomials(MB.explicit_basis(p))
             if _is_hermitian_square(mono)
                 for i in 1:div(MP.degree(mono), 2)
                     w = _chip(mono, -i)
                     if within_bounds(w, bounds)
-                        push!(monos, w)
+                        push!(keys, MP.exponents(w, MP.variables(full)))
                     end
                 end
             end
         end
-        basis = MB.SubBasis{B}(monos)
+        unique!(keys)
+        basis = SA.SubBasis(full, keys)
         return basis, typeof(basis)[]
     end
 end
@@ -543,17 +543,20 @@ end
 
 function half_newton_polytope(
     p::SA.AlgebraElement,
-    gs::AbstractVector{<:SA.AlgebraElement{<:MB.Algebra{B},T}},
+    gs::AbstractVector{G},
     vars,
     maxdegree,
     filter::NewtonFilter{<:NewtonDegreeBounds},
-) where {B,T}
+) where {G<:SA.AlgebraElement}
     basis, multipliers_bases =
         half_newton_polytope(p, gs, vars, maxdegree, filter.outer_approximation)
     bases = copy(multipliers_bases)
     push!(bases, basis)
     gs = copy(gs)
-    push!(gs, MB.constant_algebra_element(B, T))
+    push!(
+        gs,
+        MB.constant_algebra_element(MB.implicit_basis(SA.basis(p)), eltype(G)),
+    )
     filtered_bases = post_filter(p, gs, bases)
     # The last one will be recomputed by the ideal certificate
     return filtered_bases[end], filtered_bases[1:(end-1)]
@@ -668,12 +671,12 @@ function SA.unsafe_push!(c::_DictCoefficients{K}, key::K, value) where {K}
     return c
 end
 
-function _term(α, p::MB.Polynomial{B,M}) where {B,M}
-    return MB.algebra_element(MP.term(α, p.monomial), MB.FullBasis{B,M}())
+function _term(α, p::MB.Polynomial)
+    return MB.term_element(α, p)
 end
 
-function _term_constant_monomial(α, ::MB.Polynomial{B,M}) where {B,M}
-    return _term(α, MB.Polynomial{B}(MP.constant_monomial(M)))
+function _term_constant_monomial(α, p::MB.Polynomial{B,M}) where {B,M}
+    return _term(α, MB.Polynomial{B}(MP.constant_monomial(MP.monomial(p))))
 end
 
 # If `mono` is such that there is no other way to have `mono^2` by multiplying
@@ -707,8 +710,14 @@ function post_filter(
     # We use `_DictCoefficients` instead `SA.SparseCoefficients` because
     # we need to keep it canonicalized (without duplicate actually)
     # and don't care about the list of monomials being ordered
+    # TODO Base.promote_op is a bit hacky, MA.promote_operation or MP.exponents_type would be better
     counter = MB.algebra_element(
-        _DictCoefficients(Dict{MP.monomial_type(typeof(poly)),SignCount}()),
+        _DictCoefficients(
+            Dict{
+                Base.promote_op(MP.exponents, MP.monomial_type(typeof(poly))),
+                SignCount,
+            }(),
+        ),
         MB.implicit_basis(SA.basis(poly)),
     )
     cache = zero(SignCount, MB.algebra(MB.implicit_basis(SA.basis(poly))))
@@ -721,6 +730,16 @@ function post_filter(
         )
     end
     for (mult, gram_monos) in zip(generators, multipliers_gram_monos)
+        if MP.variables(mult) != MP.variables(poly)
+            error(
+                "Generator has variables $(MP.variables(mult)) and poly $(MP.variables(poly))",
+            )
+        end
+        if MP.variables(gram_monos) != MP.variables(poly)
+            error(
+                "Multipliers has variables $(MP.variables(gram_monos)) and poly $(MP.variables(poly))",
+            )
+        end
         MA.operate_to!(
             cache,
             +,
@@ -791,22 +810,18 @@ function post_filter(
         end
     end
     return [
-        _sub(gram_monos, keep) for
+        SA.sub_basis(gram_monos, keep) for
         (keep, gram_monos) in zip(keep, multipliers_gram_monos)
     ]
 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{
+        T,
         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),
@@ -830,12 +845,18 @@ end
 
 function half_newton_polytope(basis::MB.SubBasis, args...)
     a = MB.algebra_element(
-        SA.SparseCoefficients(basis.monomials, ones(length(basis))),
+        SA.SparseCoefficients(
+            basis.keys,
+            ones(length(basis)),
+            SA.comparable(MB.implicit_basis(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
+    return MB.keys_as_monomials(
+        half_newton_polytope(MB.SubBasis{MB.Monomial}(monos), args...),
+    )
 end
diff --git a/src/Certificate/preorder.jl b/src/Certificate/preorder.jl
index 4dc6994d0..5064b1b5b 100644
--- a/src/Certificate/preorder.jl
+++ b/src/Certificate/preorder.jl
@@ -32,6 +32,8 @@ end
 
 cone(certificate::Putinar) = cone(certificate.multipliers_certificate)
 
+# SemialgebraicSets.FullSpace does not implement MP.variables so we need this
+# It would be cleaner to make FullSpace have a list of variables though
 struct WithVariables{S,V}
     inner::S
     variables::V
@@ -50,39 +52,33 @@ struct WithFixedBases{S,B}
     bases::Vector{B}
 end
 
-_merge_sorted(a::Vector, ::Tuple{}) = a
-function _merge_sorted(a::Vector, b::Vector)
-    vars = sort!(vcat(a, b), rev = true)
-    unique!(vars)
-    return vars
-end
-_merge_sorted(a::Tuple{}, ::Tuple{}) = a
-_merge_sorted(a::Tuple, ::Tuple{}) = a
-_merge_sorted(::Tuple{}, b::Tuple) = b
-function _merge_sorted(a::Tuple, b::Tuple)
-    v = first(a)
-    w = first(b)
-    if v == w
-        return (v, _merge_sorted(Base.tail(a), Base.tail(b))...)
-    elseif v > w
-        return (v, _merge_sorted(Base.tail(a), b)...)
-    else
-        return (w, _merge_sorted(a, Base.tail(b))...)
-    end
+# TODO temporary workaround because SS doesn't support AlgebraElement yet
+function with_variables(p::SA.AlgebraElement, ::FullSpace)
+    return WithVariables(p, MP.variables(p))
 end
 
-_vars(::SemialgebraicSets.FullSpace) = tuple()
-function _vars(x::SA.AlgebraElement)
-    if SA.basis(x) isa SA.ImplicitBasis
-        return MP.variables(SA.coeffs(x))
-    else
-        return MP.variables(SA.basis(x))
-    end
+function with_variables(p::SA.AlgebraElement, domain)
+    _, q = SumOfSquares._promote_bases(domain, p)
+    return WithVariables(q, MP.variables(q))
+end
+
+# TODO not needed
+function with_variables(p::MP.AbstractPolynomialLike, ::FullSpace)
+    return WithVariables(p, MP.variables(p))
+end
+
+function with_variables(p::MP.AbstractPolynomialLike, domain)
+    inner, outer = SumOfSquares._promote_bases(p, domain)
+    return WithVariables(inner, MP.variables(outer))
+end
+
+function with_variables(domain, p::WithVariables)
+    return with_variables(domain, p.inner)
 end
-_vars(x) = MP.variables(x)
 
-function with_variables(inner, outer)
-    return WithVariables(inner, _merge_sorted(_vars(inner), _vars(outer)))
+function with_variables(domain, p)
+    inner, outer = SumOfSquares._promote_bases(domain, p)
+    return WithVariables(inner, MP.variables(outer))
 end
 
 function with_fixed_basis(
@@ -96,7 +92,7 @@ function with_fixed_basis(
         v.inner,
         half_newton_polytope(
             _algebra_element(p),
-            SemialgebraicSets.inequalities(domain),
+            SemialgebraicSets.inequalities(v.inner),
             v.variables,
             maxdegree,
             newton,
diff --git a/src/attributes.jl b/src/attributes.jl
index ec7dee798..8d5c7dd23 100644
--- a/src/attributes.jl
+++ b/src/attributes.jl
@@ -144,8 +144,8 @@ function MOI.Bridges.unbridged_function(
         Vector{<:GramMatrix{T}},
         BlockDiagonalGramMatrix{T},
         Vector{<:BlockDiagonalGramMatrix{T}},
-        SOSDecomposition{<:Any,T},
-        SOSDecompositionWithDomain{<:Any,T},
+        SOSDecomposition{T},
+        SOSDecompositionWithDomain{T},
         MultivariateMoments.MomentMatrix{T},
         MultivariateMoments.BlockDiagonalMomentMatrix{T},
         MultivariateMoments.MomentVector{T},
@@ -158,8 +158,8 @@ end
 # This is type piracy but we tolerate it.
 const ObjectWithoutIndex = Union{
     AbstractGramMatrix{<:MOI.Utilities.ObjectWithoutIndex},
-    SOSDecomposition{<:Any,<:MOI.Utilities.ObjectWithoutIndex},
-    SOSDecompositionWithDomain{<:Any,<:MOI.Utilities.ObjectWithoutIndex},
+    SOSDecomposition{<:MOI.Utilities.ObjectWithoutIndex},
+    SOSDecompositionWithDomain{<:MOI.Utilities.ObjectWithoutIndex},
 }
 const ObjectOrTupleWithoutIndex =
     Union{ObjectWithoutIndex,Tuple{Vararg{ObjectWithoutIndex}}}
diff --git a/src/constraints.jl b/src/constraints.jl
index 9f6b096f8..22383c10f 100644
--- a/src/constraints.jl
+++ b/src/constraints.jl
@@ -376,7 +376,11 @@ function _default_basis(
         return coeffs, basis
     else
         return _default_basis(
-            SA.SparseCoefficients(basis.monomials, coeffs),
+            SA.SparseCoefficients(
+                basis.keys,
+                coeffs,
+                SA.comparable(MB.implicit_basis(basis)),
+            ),
             MB.implicit_basis(basis),
             gram_basis,
         )
@@ -391,11 +395,11 @@ function _default_basis(
     if B === G
         return _default_basis(
             collect(SA.values(p)),
-            MB.SubBasis{B}(collect(SA.keys(p))),
+            SA.sub_basis(basis, collect(SA.keys(p))),
             gram_basis,
         )
     else
-        new_basis = MB.FullBasis{G,MP.monomial_type(typeof(basis))}()
+        new_basis = MB.FullBasis{G}(MP.variables(basis))
         return _default_basis(
             SA.coeffs(p, basis, new_basis),
             new_basis,
@@ -404,18 +408,15 @@ function _default_basis(
     end
 end
 
-function _default_gram_basis(
-    ::MB.MonomialIndexedBasis{B,M},
-    ::Nothing,
-) where {B,M}
-    return MB.FullBasis{B,M}()
+function _default_gram_basis(b::MB.MonomialIndexedBasis, ::Nothing)
+    return MB.implicit_basis(b)
 end
 
 function _default_gram_basis(
-    ::MB.MonomialIndexedBasis{_B,M},
+    b::MB.MonomialIndexedBasis{_B,M},
     ::Type{B},
 ) where {_B,B,M}
-    return MB.FullBasis{B,M}()
+    return MB.FullBasis{B}(MP.variables(b))
 end
 
 function _default_gram_basis(_, basis::MB.MonomialIndexedBasis)
@@ -432,7 +433,7 @@ function _default_basis(p::MP.AbstractPolynomialLike, basis)
     return _default_basis(
         MB.algebra_element(
             MB.sparse_coefficients(p),
-            MB.FullBasis{MB.Monomial,MP.monomial_type(p)}(),
+            MB.FullBasis{MB.Monomial}(p),
         ),
         basis,
     )
@@ -443,20 +444,31 @@ function _default_zero_basis(basis, nodes::MB.AbstractNodes)
     return MB.ImplicitLagrangeBasis(MP.variables(basis), nodes)
 end
 
+_promote_bases(domain, p) = SA.promote_bases(domain, p)
+_promote_bases(domain::FullSpace, p::SA.AlgebraElement) = domain, p
+function _promote_bases(domain, p::SA.AlgebraElement)
+    _domain, _ = SA.promote_bases(domain, MP.polynomial(p))
+    _, _p = SA.promote_bases(MB.algebra_element(prod(MP.variables(domain))), p)
+    return _domain, _p
+end
+
 function JuMP.build_constraint(
     _error::Function,
     p,
     cone::SOSLikeCone;
     basis = nothing,
     zero_basis = nothing,
+    domain::AbstractSemialgebraicSet = FullSpace(),
     kws...,
 )
+    domain, p = _promote_bases(domain, p)
     __coefs, basis, gram_basis = _default_basis(p, basis)
     set = JuMP.moi_set(
         cone,
         basis,
         gram_basis,
         _default_zero_basis(basis, zero_basis);
+        domain,
         kws...,
     )
     _coefs = PolyJuMP.non_constant(__coefs)
@@ -556,7 +568,7 @@ Return the monomials of [`certificate_basis`](@ref). If the basis if not
 function certificate_monomials(cref::JuMP.ConstraintRef)
     return basis_monomials(certificate_basis(cref))
 end
-basis_monomials(basis::MB.SubBasis) = basis.monomials
+basis_monomials(basis::MB.SubBasis) = keys_as_monomials(basis)
 function basis_monomials(basis::SA.AbstractBasis)
     return error(
         "`certificate_monomials` is not supported with `$(typeof(basis))`, use `certificate_basis` instead.",
diff --git a/src/gram_matrix.jl b/src/gram_matrix.jl
index 9cd892caa..42a684372 100644
--- a/src/gram_matrix.jl
+++ b/src/gram_matrix.jl
@@ -159,7 +159,8 @@ function gram_operate(
     p::GramMatrix{S,B,US,SymMatrix{S}},
     q::GramMatrix{T,B,UT,SymMatrix{T}},
 ) where {S,US,T,UT,B}
-    basis, Ip, Iq = MultivariateBases.merge_bases(p.basis, q.basis)
+    pb, qb = SA.promote_bases(p.basis, q.basis)
+    basis, Ip, Iq = SA.merge_bases_with_maps(pb, qb)
     U = MA.promote_operation(+, S, T)
     n = length(basis)
     Qvec = Vector{U}(undef, div(n * (n + 1), 2))
@@ -211,13 +212,29 @@ end
 struct BlockDiagonalGramMatrix{T,B,U,MT} <: AbstractGramMatrix{T,B,U}
     blocks::Vector{GramMatrix{T,B,U,MT}}
 end
+function BlockDiagonalGramMatrix(blocks::Vector{Any})
+    if isempty(blocks)
+        return BlockDiagonalGramMatrix(
+            GramMatrix{Float64,Nothing,Float64,Matrix{Float64}}[],
+        )
+    end
+    b1 = first(blocks)::GramMatrix
+    T, B, U, MT = typeof(b1).parameters
+    return BlockDiagonalGramMatrix(
+        convert(Vector{GramMatrix{T,B,U,MT}}, blocks),
+    )
+end
 
 function MB.implicit_basis(g::BlockDiagonalGramMatrix)
     return MB.implicit_basis(first(g.blocks))
 end
 
-function MultivariateMoments.block_diagonal(blocks::Vector{<:GramMatrix})
-    return BlockDiagonalGramMatrix(blocks)
+function MultivariateMoments.block_diagonal(
+    blocks::Vector{<:GramMatrix{T,B,U,MT}},
+) where {T,B,U,MT}
+    return BlockDiagonalGramMatrix(
+        convert(Vector{GramMatrix{T,B,U,MT}}, blocks),
+    )
 end
 
 function _sparse_type(::Type{GramMatrix{T,B,U,MT}}) where {T,B,U,MT}
diff --git a/src/sosdec.jl b/src/sosdec.jl
index 427e0af56..088adf35e 100644
--- a/src/sosdec.jl
+++ b/src/sosdec.jl
@@ -1,44 +1,47 @@
 export SOSDecomposition, SOSDecompositionWithDomain, sos_decomposition
 
 """
-    struct SOSDecomposition{A,T,V,U}
+    struct SOSDecomposition{T,A,V,U}
 
 Represents a Sum-of-Squares decomposition without domain.
 """
-struct SOSDecomposition{A,T,V,U} <: AbstractDecomposition{U}
-    ps::Vector{SA.AlgebraElement{A,T,V}} # TODO rename `elements`
-    function SOSDecomposition{A,T,V,U}(
-        ps::Vector{SA.AlgebraElement{A,T,V}},
-    ) where {A,T,V,U}
+struct SOSDecomposition{T,A,V,U} <: AbstractDecomposition{U}
+    ps::Vector{SA.AlgebraElement{T,A,V}} # TODO rename `elements`
+    function SOSDecomposition{T,A,V,U}(
+        ps::Vector{SA.AlgebraElement{T,A,V}},
+    ) where {T,A,V,U}
         return new(ps)
     end
 end
 
 function SOSDecomposition(
-    elements::Vector{SA.AlgebraElement{A,T,V}},
-) where {A,T,V}
-    return SOSDecomposition{A,T,V,_promote_add_mul(T)}(elements)
+    elements::Vector{SA.AlgebraElement{T,A,V}},
+) where {T,A,V}
+    return SOSDecomposition{T,A,V,_promote_add_mul(T)}(elements)
 end
 function MP.polynomial_type(
-    ::Union{SOSDecomposition{A,T,V,U},Type{SOSDecomposition{A,T,V,U}}},
+    ::Union{SOSDecomposition{T,A,V,U},Type{SOSDecomposition{T,A,V,U}}},
 ) where {A,T,V,U}
-    return MP.polynomial_type(MP.polynomial_type(SA.AlgebraElement{A,T,V}), U)
+    return MP.polynomial_type(MP.polynomial_type(SA.AlgebraElement{T,A,V}), U)
 end
 
-function GramMatrix(p::SOSDecomposition{A,T}) where {A,T}
-    basis = mapreduce(SA.basis, (b1, b2) -> MB.merge_bases(b1, b2)[1], p.ps)
+function GramMatrix(p::SOSDecomposition{T}) where {T}
+    sub_bases = [MB.explicit_basis(a) for a in p.ps]
+    basis = reduce(
+        (b1, b2) -> SA.merge_bases(SA.promote_bases(b1, b2)...),
+        sub_bases,
+    )
     m = length(p.ps)
     n = length(basis)
     Q = zeros(T, m, n)
     for i in 1:m
-        j = 1
         for (k, v) in SA.nonzero_pairs(SA.coeffs(p.ps[i]))
-            poly = SA.basis(p.ps[i])[k]
-            while j in eachindex(basis) && basis[j] != poly
-                j += 1
+            # k is the key (exponent vector) in the algebra element's basis
+            # Find position j in merged basis with the same key
+            idx = findfirst(==(k), basis.keys)
+            if idx !== nothing
+                Q[i, idx] = v
             end
-            Q[i, j] = v
-            j += 1
         end
     end
     return GramMatrix(Q' * Q, basis)
@@ -102,26 +105,26 @@ function Base.isapprox(p::SOSDecomposition, q::SOSDecomposition; kwargs...)
 end
 
 function Base.promote_rule(
-    ::Type{SOSDecomposition{A,T1,V1,U1}},
-    ::Type{SOSDecomposition{A,T2,V2,U2}},
-) where {A,T1,T2,V1,V2,U1,U2}
+    ::Type{SOSDecomposition{T1,A,V1,U1}},
+    ::Type{SOSDecomposition{T2,A,V2,U2}},
+) where {T1,T2,A,V1,V2,U1,U2}
     T = promote_type(T1, T2)
     V = SA.similar_type(V1, T)
-    return SOSDecomposition{A,T,V,_promote_add_mul(T)}
+    return SOSDecomposition{T,A,V,_promote_add_mul(T)}
 end
 
 function Base.convert(
-    ::Type{SOSDecomposition{A,T,V,U}},
+    ::Type{SOSDecomposition{T,A,V,U}},
     p::SOSDecomposition,
-) where {A,T,V,U}
-    return SOSDecomposition(convert(Vector{SA.AlgebraElement{A,T,V}}, p.ps))
+) where {T,A,V,U}
+    return SOSDecomposition(convert(Vector{SA.AlgebraElement{T,A,V}}, p.ps))
 end
 
 function MP.polynomial(d::SOSDecomposition)
     return MP.polynomial(MB.algebra_element(d))
 end
 
-function MB.algebra_element(decomp::SOSDecomposition{A,T,V,U}) where {A,T,V,U}
+function MB.algebra_element(decomp::SOSDecomposition{T,A,V,U}) where {T,A,V,U}
     basis = MB.implicit_basis(SA.basis(first(decomp.ps)))
     res = zero(U, MB.algebra(basis))
     for p in decomp.ps
@@ -141,20 +144,20 @@ end
 
 Represents a Sum-of-Squares decomposition on a basic semi-algebraic domain.
 """
-struct SOSDecompositionWithDomain{A,T,V,U,S<:AbstractSemialgebraicSet}
-    sos::SOSDecomposition{A,T,V,U}
-    sosj::Vector{SOSDecomposition{A,T,V,U}}
+struct SOSDecompositionWithDomain{T,A,V,U,S<:AbstractSemialgebraicSet}
+    sos::SOSDecomposition{T,A,V,U}
+    sosj::Vector{SOSDecomposition{T,A,V,U}}
     domain::S
 end
 
 function SOSDecompositionWithDomain(
-    ps::SOSDecomposition{A1,T1,V1,U1},
-    vps::Vector{SOSDecomposition{A2,T2,V2,U2}},
+    ps::SOSDecomposition{T1,A1,V1,U1},
+    vps::Vector{SOSDecomposition{T2,A2,V2,U2}},
     set::AbstractSemialgebraicSet,
 ) where {A1,A2,T1,T2,V1,V2,U1,U2}
     ptype = promote_type(
-        SOSDecomposition{A1,T1,V1,U1},
-        SOSDecomposition{A2,T2,V2,U2},
+        SOSDecomposition{T1,A1,V1,U1},
+        SOSDecomposition{T2,A2,V2,U2},
     )
     return SOSDecompositionWithDomain(
         convert(ptype, ps),
diff --git a/test/Project.toml b/test/Project.toml
index b821dc5f9..c674ebf11 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -1,4 +1,5 @@
 [deps]
+Clarabel = "61c947e1-3e6d-4ee4-985a-eec8c727bd6e"
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
 DynamicPolynomials = "7c1d4256-1411-5781-91ec-d7bc3513ac07"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
@@ -15,6 +16,10 @@ StarAlgebras = "0c0c59c1-dc5f-42e9-9a8b-b5dc384a6cd1"
 SumOfSquares = "4b9e565b-77fc-50a5-a571-1244f986bda1"
 SymbolicWedderburn = "858aa9a9-4c7c-4c62-b466-2421203962a2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+TypedPolynomials = "afbbf031-7a57-5f58-a1b9-b774a0fad08d"
 
 [compat]
-DynamicPolynomials = "0.5,0.6"
+DynamicPolynomials = "0.6"
+
+[sources]
+SumOfSquares = {path = ".."}
diff --git a/test/Tests/BPT12e399.jl b/test/Tests/BPT12e399.jl
index aba307047..4a8233ae7 100644
--- a/test/Tests/BPT12e399.jl
+++ b/test/Tests/BPT12e399.jl
@@ -1,4 +1,5 @@
 # Adapted from Example 3.99 of [BPT12].
+import MultivariatePolynomials as MP
 #
 # [BPT12] Blekherman, G.; Parrilo, P. & Thomas, R.
 # Semidefinite Optimization and Convex Algebraic Geometry
@@ -50,14 +51,14 @@ function BPT12e399_test(optimizer, config::MOI.Test.Config, remainder::Bool)
     p = gram_matrix(cref)
     if remainder
         @test value_matrix(p) ≈ [9 -3; -3 1] atol = atol rtol = rtol
-        @test p.basis.monomials == [1, y]
+        @test MB.keys_as_monomials(p.basis) == [1, y]
     else
         @test value_matrix(p) ≈ [
             5 -5 0
             -5 5 0
             0 0 4
         ] atol = atol rtol = rtol
-        @test p.basis.monomials == [1, y, x]
+        @test MB.keys_as_monomials(p.basis) == [1, y, x]
     end
 
     @test dual_status(model) == MOI.FEASIBLE_POINT
@@ -66,35 +67,35 @@ function BPT12e399_test(optimizer, config::MOI.Test.Config, remainder::Bool)
     @test length(moments(μ)) == 3
     @test moment_value(moments(μ)[1]) ≈ (remainder ? 1 / 3 : 1.0) atol = atol rtol =
         rtol
-    @test moments(μ)[1].polynomial.monomial == 1
+    @test MP.monomial(moments(μ)[1].polynomial) == 1
     @test moment_value(moments(μ)[2]) ≈ 1 atol = atol rtol = rtol
-    @test moments(μ)[2].polynomial.monomial == y
+    @test MP.monomial(moments(μ)[2].polynomial) == y
     @test moment_value(moments(μ)[3]) ≈ (remainder ? -8 / 3 : 0.0) atol = atol rtol =
         rtol
-    @test moments(μ)[3].polynomial.monomial == x^2
+    @test MP.monomial(moments(μ)[3].polynomial) == x^2
 
     μ = moments(cref)
     @test μ isa AbstractMeasure{Float64}
     if remainder
         @test length(moments(μ)) == 3
         @test moment_value(moments(μ)[1]) ≈ 1 / 3 atol = atol rtol = rtol
-        @test moments(μ)[1].polynomial.monomial == 1
+        @test MP.monomial(moments(μ)[1].polynomial) == 1
         @test moment_value(moments(μ)[2]) ≈ 1 atol = atol rtol = rtol
-        @test moments(μ)[2].polynomial.monomial == y
+        @test MP.monomial(moments(μ)[2].polynomial) == y
         @test moment_value(moments(μ)[3]) ≈ 3 atol = atol rtol = rtol
-        @test moments(μ)[3].polynomial.monomial == y^2
+        @test MP.monomial(moments(μ)[3].polynomial) == y^2
     else
         @test length(moments(μ)) == 5
         @test moment_value(moments(μ)[1]) ≈ 1 atol = atol rtol = rtol
-        @test moments(μ)[1].polynomial.monomial == 1
+        @test MP.monomial(moments(μ)[1].polynomial) == 1
         @test moment_value(moments(μ)[2]) ≈ 1 atol = atol rtol = rtol
-        @test moments(μ)[2].polynomial.monomial == y
+        @test MP.monomial(moments(μ)[2].polynomial) == y
         @test moment_value(moments(μ)[3]) ≈ 0 atol = atol rtol = rtol
-        @test moments(μ)[3].polynomial.monomial == x
+        @test MP.monomial(moments(μ)[3].polynomial) == x
         @test moment_value(moments(μ)[4]) ≈ 1 atol = atol rtol = rtol
-        @test moments(μ)[4].polynomial.monomial == y^2
+        @test MP.monomial(moments(μ)[4].polynomial) == y^2
         @test moment_value(moments(μ)[5]) ≈ 0 atol = atol rtol = rtol
-        @test moments(μ)[5].polynomial.monomial == x * y
+        @test MP.monomial(moments(μ)[5].polynomial) == x * y
     end
 
     @objective(model, Min, α)
@@ -111,14 +112,14 @@ function BPT12e399_test(optimizer, config::MOI.Test.Config, remainder::Bool)
     p = gram_matrix(cref)
     if remainder
         @test value_matrix(p) ≈ [9 3; 3 1] atol = atol rtol = rtol
-        @test p.basis.monomials == [1, y]
+        @test MB.keys_as_monomials(p.basis) == [1, y]
     else
         @test value_matrix(p) ≈ [
             5 5 0
             5 5 0
             0 0 4
         ] atol = atol rtol = rtol
-        @test p.basis.monomials == [1, y, x]
+        @test MB.keys_as_monomials(p.basis) == [1, y, x]
     end
 
     @test dual_status(model) == MOI.FEASIBLE_POINT
@@ -127,35 +128,35 @@ function BPT12e399_test(optimizer, config::MOI.Test.Config, remainder::Bool)
     @test length(moments(μ)) == 3
     @test moment_value(moments(μ)[1]) ≈ (remainder ? 1 / 3 : 1.0) atol = atol rtol =
         rtol
-    @test moments(μ)[1].polynomial.monomial == 1
+    @test MP.monomial(moments(μ)[1].polynomial) == 1
     @test moment_value(moments(μ)[2]) ≈ -1 atol = atol rtol = rtol
-    @test moments(μ)[2].polynomial.monomial == y
+    @test MP.monomial(moments(μ)[2].polynomial) == y
     @test moment_value(moments(μ)[3]) ≈ (remainder ? -8 / 3 : 0.0) atol = atol rtol =
         rtol
-    @test moments(μ)[3].polynomial.monomial == x^2
+    @test MP.monomial(moments(μ)[3].polynomial) == x^2
 
     μ = moments(cref)
     @test μ isa AbstractMeasure{Float64}
     if remainder
         @test length(moments(μ)) == 3
         @test moment_value(moments(μ)[1]) ≈ 1 / 3 atol = atol rtol = rtol
-        @test moments(μ)[1].polynomial.monomial == 1
+        @test MP.monomial(moments(μ)[1].polynomial) == 1
         @test moment_value(moments(μ)[2]) ≈ -1 atol = atol rtol = rtol
-        @test moments(μ)[2].polynomial.monomial == y
+        @test MP.monomial(moments(μ)[2].polynomial) == y
         @test moment_value(moments(μ)[3]) ≈ 3 atol = atol rtol = rtol
-        @test moments(μ)[3].polynomial.monomial == y^2
+        @test MP.monomial(moments(μ)[3].polynomial) == y^2
     else
         @test length(moments(μ)) == 5
         @test moment_value(moments(μ)[1]) ≈ 1 atol = atol rtol = rtol
-        @test moments(μ)[1].polynomial.monomial == 1
+        @test MP.monomial(moments(μ)[1].polynomial) == 1
         @test moment_value(moments(μ)[2]) ≈ -1 atol = atol rtol = rtol
-        @test moments(μ)[2].polynomial.monomial == y
+        @test MP.monomial(moments(μ)[2].polynomial) == y
         @test moment_value(moments(μ)[3]) ≈ 0 atol = atol rtol = rtol
-        @test moments(μ)[3].polynomial.monomial == x
+        @test MP.monomial(moments(μ)[3].polynomial) == x
         @test moment_value(moments(μ)[4]) ≈ 1 atol = atol rtol = rtol
-        @test moments(μ)[4].polynomial.monomial == y^2
+        @test MP.monomial(moments(μ)[4].polynomial) == y^2
         @test moment_value(moments(μ)[5]) ≈ 0 atol = atol rtol = rtol
-        @test moments(μ)[5].polynomial.monomial == x * y
+        @test MP.monomial(moments(μ)[5].polynomial) == x * y
     end
 
     return model
diff --git a/test/Tests/quadratic.jl b/test/Tests/quadratic.jl
index 074b1df30..e138118df 100644
--- a/test/Tests/quadratic.jl
+++ b/test/Tests/quadratic.jl
@@ -1,4 +1,5 @@
 using Test
+import MultivariatePolynomials as MP
 import MultivariateBases
 using DynamicPolynomials
 
@@ -6,7 +7,7 @@ function _test_moments(test_values, μ, monos)
     @test μ isa AbstractMeasure{Float64}
     @test length(moments(μ)) == length(monos)
     test_values(moment_value.(moments(μ)))
-    @test [m.polynomial.monomial for m in moments(μ)] == monos
+    @test [MP.monomial(m.polynomial) for m in moments(μ)] == monos
 end
 
 function quadratic_test(
@@ -66,7 +67,7 @@ function quadratic_test(
     else
         @test value_matrix(p) ≈ ones(2, 2) atol = atol rtol = rtol
     end
-    @test p.basis.monomials == cert_monos
+    @test MB.keys_as_monomials(p.basis) == cert_monos
 
     μ = moments(dual(cref))
     a = moment_value.(μ)
@@ -81,7 +82,7 @@ function quadratic_test(
         atol rtol = rtol
     @test b[1] + b[3] ≈ 2.0 atol = atol rtol = rtol
     @test μ[2].polynomial == MB.Polynomial{basis}(
-        bivariate ? (basis === MB.Chebyshev ? y^2 : x * y) : x^1,
+        bivariate ? (basis === MB.Chebyshev ? x^0*y^2 : x * y) : x^1,
     )
 
     @test dual_status(model) == MOI.FEASIBLE_POINT
@@ -122,20 +123,17 @@ function quadratic_test(
         b[1] off
         off b[3]
     ] atol = atol rtol = rtol
-    @test ν.basis.monomials == cert_monos
+    @test MB.keys_as_monomials(ν.basis) == cert_monos
 
+    _FB = typeof(MB.FullBasis{basis}(x))
+    _SB = MB.explicit_basis_type(_FB)
     N = SumOfSquares.Certificate.NewtonFilter{
         SumOfSquares.Certificate.NewtonDegreeBounds{Tuple{}},
     }
     S = SumOfSquares.SOSPolynomialSet{
         SumOfSquares.FullSpace,
-        MB.SubBasis{basis,monomial_type(x),monomial_vector_type(x)},
-        SumOfSquares.Certificate.Newton{
-            typeof(cone),
-            MB.FullBasis{basis,monomial_type(x)},
-            MB.FullBasis{basis,monomial_type(x)},
-            N,
-        },
+        _SB,
+        SumOfSquares.Certificate.Newton{typeof(cone),_FB,_FB,N},
     }
     @test list_of_constraint_types(model) == [(Vector{AffExpr}, S)]
     return test_delete_bridge(
diff --git a/test/Tests/quartic_constant.jl b/test/Tests/quartic_constant.jl
index 7e4b34537..48033fa5e 100644
--- a/test/Tests/quartic_constant.jl
+++ b/test/Tests/quartic_constant.jl
@@ -34,16 +34,13 @@ function quartic_constant_test(
     @test p isa SumOfSquares.GramMatrix
     @test value_matrix(p) ≈ ones(1, 1) atol = atol rtol = rtol
     @test p.basis isa MB.SubBasis{MB.Monomial}
-    @test p.basis.monomials == [x^2]
+    @test MB.keys_as_monomials(p.basis) == [x^2]
 
+    _FB = typeof(MB.FullBasis{MB.Monomial}(x))
     S = SumOfSquares.SOSPolynomialSet{
         SumOfSquares.FullSpace,
-        MB.SubBasis{MB.Monomial,monomial_type(x),monomial_vector_type(x)},
-        SumOfSquares.Certificate.FixedBasis{
-            typeof(cone),
-            typeof(p.basis),
-            MB.FullBasis{MB.Monomial,monomial_type(x)},
-        },
+        MB.explicit_basis_type(_FB),
+        SumOfSquares.Certificate.FixedBasis{typeof(cone),typeof(p.basis),_FB},
     }
     @test list_of_constraint_types(model) == [(Vector{AffExpr}, S)]
     return test_delete_bridge(
diff --git a/test/Tests/quartic_ideal.jl b/test/Tests/quartic_ideal.jl
index b9ac96169..669aeb0a5 100644
--- a/test/Tests/quartic_ideal.jl
+++ b/test/Tests/quartic_ideal.jl
@@ -1,4 +1,5 @@
 using Test
+import MultivariatePolynomials as MP
 using SumOfSquares
 using DynamicPolynomials
 
@@ -33,14 +34,14 @@ function quartic_ideal_test(
         @test termination_status(model) == MOI.OPTIMAL
         @test primal_status(model) == MOI.FEASIBLE_POINT
         μ = dual(cref)
-        @test moments(μ)[1].polynomial.monomial == 1
-        @test moments(μ)[2].polynomial.monomial == x^2
-        @test moments(μ)[3].polynomial.monomial == x^4
+        @test MP.monomial(moments(μ)[1].polynomial) == 1
+        @test MP.monomial(moments(μ)[2].polynomial) == x^2
+        @test MP.monomial(moments(μ)[3].polynomial) == x^4
         μ = moments(cref)
-        @test moments(μ)[1].polynomial.monomial == 1
-        @test moments(μ)[2].polynomial.monomial == x^1
-        @test moments(μ)[3].polynomial.monomial == x^2
-        @test moment_matrix(cref).basis.monomials == [1, x, x^2]
+        @test MP.monomial(moments(μ)[1].polynomial) == 1
+        @test MP.monomial(moments(μ)[2].polynomial) == x^1
+        @test MP.monomial(moments(μ)[3].polynomial) == x^2
+        @test MB.keys_as_monomials(moment_matrix(cref).basis) == [1, x, x^2]
     end
 end
 function quartic_ideal_test(optimizer, config)
diff --git a/test/Tests/rearrangement.jl b/test/Tests/rearrangement.jl
index 126328cd4..163c9dce9 100644
--- a/test/Tests/rearrangement.jl
+++ b/test/Tests/rearrangement.jl
@@ -39,34 +39,34 @@ function rearrangement_test(optimizer, config::MOI.Test.Config)
         0 1 -1
         0 -1 1
     ] atol = 9atol rtol = 9rtol
-    @test p.blocks[1].basis.monomials == [1, y, x]
+    @test MB.keys_as_monomials(p.blocks[1].basis) == [1, y, x]
     @test value_matrix(p.blocks[2]) ≈ [
         0 0 0
         0 1 -1
         0 -1 1
     ] atol = 9atol rtol = 9rtol
-    @test p.blocks[2].basis.monomials == [1, z, y]
+    @test MB.keys_as_monomials(p.blocks[2].basis) == [1, z, y]
 
     λ = lagrangian_multipliers(cref)
     @test length(λ) == 4
     @test λ[1] isa SumOfSquares.BlockDiagonalGramMatrix
     @test length(λ[1].blocks) == 1
-    @test λ[1].blocks[1].basis.monomials == [x, y, 1]
+    @test MB.keys_as_monomials(λ[1].blocks[1].basis) == [x, y, 1]
     @test λ[2] isa SumOfSquares.BlockDiagonalGramMatrix
     @test length(λ[2].blocks) == 1
-    @test λ[2].blocks[1].basis.monomials == [y, z, 1]
+    @test MB.keys_as_monomials(λ[2].blocks[1].basis) == [y, z, 1]
     @test λ[3] isa SumOfSquares.BlockDiagonalGramMatrix
     @test length(λ[3].blocks) == 1
-    @test λ[3].blocks[1].basis.monomials == [x, y, 1]
+    @test MB.keys_as_monomials(λ[3].blocks[1].basis) == [x, y, 1]
     @test λ[4] isa SumOfSquares.BlockDiagonalGramMatrix
     @test length(λ[4].blocks) == 1
-    @test λ[4].blocks[1].basis.monomials == [y, z, 1]
+    @test MB.keys_as_monomials(λ[4].blocks[1].basis) == [y, z, 1]
 
     ν = moment_matrix(cref)
     @test ν isa BlockDiagonalMomentMatrix
     @test length(ν.blocks) == 2
-    @test ν.blocks[1].basis.monomials == [1, y, x]
-    @test ν.blocks[2].basis.monomials == [1, z, y]
+    @test MB.keys_as_monomials(ν.blocks[1].basis) == [1, y, x]
+    @test MB.keys_as_monomials(ν.blocks[2].basis) == [1, z, y]
 
     return model
 end
diff --git a/test/Tests/simple_matrix.jl b/test/Tests/simple_matrix.jl
index b925e41d0..4dbf5aa7b 100644
--- a/test/Tests/simple_matrix.jl
+++ b/test/Tests/simple_matrix.jl
@@ -20,7 +20,7 @@ function simple_matrix_test(optimizer, ::MOI.Test.Config)
     JuMP.optimize!(mat_model)
     @test JuMP.termination_status(mat_model) == MOI.OPTIMAL
     @test JuMP.primal_status(mat_model) == MOI.FEASIBLE_POINT
-    @test length(gram_matrix(mat_cref).basis.monomials) == 3
+    @test length(gram_matrix(mat_cref).basis) == 3
 
     # Example 3.79
     @polyvar y[1:2]
@@ -30,6 +30,6 @@ function simple_matrix_test(optimizer, ::MOI.Test.Config)
     JuMP.optimize!(model)
     @test JuMP.termination_status(model) == MOI.OPTIMAL
     @test JuMP.primal_status(model) == MOI.FEASIBLE_POINT
-    @test length(gram_matrix(cref).basis.monomials) == 3
+    @test length(gram_matrix(cref).basis) == 3
 end
 sd_tests["simple_matrix"] = simple_matrix_test
diff --git a/test/Tests/term.jl b/test/Tests/term.jl
index 991d2e60c..668f43773 100644
--- a/test/Tests/term.jl
+++ b/test/Tests/term.jl
@@ -1,5 +1,6 @@
 using Test
 import MultivariateBases as MB
+import MultivariatePolynomials as MP
 using SumOfSquares
 using DynamicPolynomials
 
@@ -32,32 +33,29 @@ function term_test(
 
     p = gram_matrix(cref)
     @test value_matrix(p) ≈ zeros(1, 1) atol = atol rtol = rtol
-    @test p.basis.monomials == [x]
+    @test MB.keys_as_monomials(p.basis) == [x]
 
     @test dual_status(model) == MOI.FEASIBLE_POINT
     for μ in [dual(cref), moments(cref)]
         @test μ isa AbstractMeasure{Float64}
         @test length(moments(μ)) == 1
         @test moment_value(moments(μ)[1]) ≈ 1.0 atol = atol rtol = rtol
-        @test moments(μ)[1].polynomial.monomial == x^2
+        @test MP.monomial(moments(μ)[1].polynomial) == x^2
     end
 
     ν = moment_matrix(cref)
     @test value_matrix(ν) ≈ ones(1, 1) atol = atol rtol = rtol
-    @test ν.basis.monomials == [x]
+    @test MB.keys_as_monomials(ν.basis) == [x]
 
+    _FB = typeof(MB.FullBasis{MB.Monomial}(x))
+    _SB = MB.explicit_basis_type(_FB)
     N = SumOfSquares.Certificate.NewtonFilter{
         SumOfSquares.Certificate.NewtonDegreeBounds{Tuple{}},
     }
     S = SumOfSquares.SOSPolynomialSet{
         SumOfSquares.FullSpace,
-        SubBasis{MB.Monomial,monomial_type(x),monomial_vector_type(x)},
-        SumOfSquares.Certificate.Newton{
-            typeof(cone),
-            FullBasis{MB.Monomial,monomial_type(x)},
-            FullBasis{MB.Monomial,monomial_type(x)},
-            N,
-        },
+        _SB,
+        SumOfSquares.Certificate.Newton{typeof(cone),_FB,_FB,N},
     }
     @test list_of_constraint_types(model) == [(Vector{VariableRef}, S)]
     return test_delete_bridge(
@@ -70,11 +68,7 @@ function term_test(
                 MOI.VectorAffineFunction{Float64},
                 SumOfSquares.PolyJuMP.ZeroPolynomialSet{
                     SumOfSquares.FullSpace,
-                    SubBasis{
-                        MB.Monomial,
-                        monomial_type(x),
-                        monomial_vector_type(x),
-                    },
+                    _SB,
                 },
                 0,
             ),
diff --git a/test/Tests/term_fixed.jl b/test/Tests/term_fixed.jl
index 05c4a4d3b..a62e6a587 100644
--- a/test/Tests/term_fixed.jl
+++ b/test/Tests/term_fixed.jl
@@ -1,4 +1,6 @@
 using Test
+import MultivariatePolynomials as MP
+import MultivariateBases as MB
 using SumOfSquares
 using DynamicPolynomials
 
@@ -32,33 +34,30 @@ function term_fixed_test(
 
     p = gram_matrix(cref)
     @test value_matrix(p) ≈ zeros(1, 1) atol = atol rtol = rtol
-    @test p.basis.monomials == [x]
+    @test MB.keys_as_monomials(p.basis) == [x]
 
     @test dual_status(model) == MOI.FEASIBLE_POINT
     for (m, μ) in [(x^2, moments(cref))]
         @test μ isa AbstractMeasure{Float64}
         @test length(moments(μ)) == 1
         @test moment_value(moments(μ)[1]) ≈ 1.0 atol = atol rtol = rtol
-        @test moments(μ)[1].polynomial.monomial == m
+        @test MP.monomial(moments(μ)[1].polynomial) == m
     end
 
     ν = moment_matrix(cref)
     @test value_matrix(ν) ≈ ones(1, 1) atol = atol rtol = rtol
-    @test ν.basis.monomials == [x]
+    @test MB.keys_as_monomials(ν.basis) == [x]
 
+    _FB = typeof(MB.FullBasis{MB.Monomial}(x))
+    _SB = MB.explicit_basis_type(_FB)
     N = SumOfSquares.Certificate.NewtonFilter{
         SumOfSquares.Certificate.NewtonDegreeBounds{Tuple{}},
     }
     S = SumOfSquares.SOSPolynomialSet{
         typeof(set),
-        SubBasis{MB.Monomial,monomial_type(x),monomial_vector_type(x)},
+        _SB,
         SumOfSquares.Certificate.Remainder{
-            SumOfSquares.Certificate.Newton{
-                typeof(cone),
-                FullBasis{MB.Monomial,monomial_type(x)},
-                FullBasis{MB.Monomial,monomial_type(x)},
-                N,
-            },
+            SumOfSquares.Certificate.Newton{typeof(cone),_FB,_FB,N},
         },
     }
     @test list_of_constraint_types(model) == [(Vector{JuMP.AffExpr}, S)]
@@ -70,14 +69,7 @@ function term_fixed_test(
             (MOI.VectorOfVariables, MOI.Nonnegatives, 0),
             (
                 MOI.VectorAffineFunction{Float64},
-                SumOfSquares.PolyJuMP.ZeroPolynomialSet{
-                    typeof(set),
-                    SubBasis{
-                        MB.Monomial,
-                        monomial_type(x),
-                        monomial_vector_type(x),
-                    },
-                },
+                SumOfSquares.PolyJuMP.ZeroPolynomialSet{typeof(set),_SB},
                 0,
             ),
         ),
diff --git a/test/Tests/univariate_sum.jl b/test/Tests/univariate_sum.jl
index 1060a84b0..4f56b3f73 100644
--- a/test/Tests/univariate_sum.jl
+++ b/test/Tests/univariate_sum.jl
@@ -30,20 +30,18 @@ function univariate_sum_test(
     @test p isa SumOfSquares.BlockDiagonalGramMatrix
     @test length(p.blocks) == 2
     @test value_matrix(p.blocks[1]) ≈ [1 -1; -1 1] atol = atol rtol = rtol
-    @test p.blocks[1].basis.monomials == [1, x]
+    @test MB.keys_as_monomials(p.blocks[1].basis) == [1, x]
     @test value_matrix(p.blocks[2]) ≈ ones(2, 2) atol = atol rtol = rtol
-    @test p.blocks[2].basis.monomials == [1, y]
+    @test MB.keys_as_monomials(p.blocks[2].basis) == [1, y]
 
+    _FB = typeof(MB.FullBasis{MB.Monomial}(x))
+    _SB = MB.explicit_basis_type(_FB)
     S = SumOfSquares.SOSPolynomialSet{
         SumOfSquares.FullSpace,
-        MB.SubBasis{MB.Monomial,monomial_type(x),monomial_vector_type(x)},
+        _SB,
         SumOfSquares.Certificate.Sparsity.Ideal{
             Sparsity.Variable,
-            SumOfSquares.Certificate.MaxDegree{
-                typeof(cone),
-                MB.FullBasis{MB.Monomial,monomial_type(x)},
-                MB.FullBasis{MB.Monomial,monomial_type(x)},
-            },
+            SumOfSquares.Certificate.MaxDegree{typeof(cone),_FB,_FB},
         },
     }
     @test list_of_constraint_types(model) == [(Vector{AffExpr}, S)]
diff --git a/test/ceg_test.jl b/test/ceg_test.jl
index e60cf6ec7..b1e2b29f1 100644
--- a/test/ceg_test.jl
+++ b/test/ceg_test.jl
@@ -58,10 +58,12 @@ const CEG = SumOfSquares.Certificate.Sparsity.ChordalExtensionGraph
         CEG.add_edge!(G, 1, 4)
 
         H, cliques = CEG.chordal_extension(G, CEG.GreedyFillIn())
-        @test cliques == [[4, 1, 3], [2, 1, 3], [5]]
+        @test Set(Set.(cliques)) ==
+              Set([Set([4, 1, 3]), Set([2, 1, 3]), Set([5])])
 
         CEG.add_edge!(H, 1, 5)
         I, cliques = CEG.chordal_extension(H, CEG.GreedyFillIn())
-        @test cliques == [[5, 1], [4, 1, 3], [2, 1, 3]]
+        @test Set(Set.(cliques)) ==
+              Set([Set([5, 1]), Set([4, 1, 3]), Set([2, 1, 3])])
     end
 end
diff --git a/test/certificate.jl b/test/certificate.jl
index c36422079..109b8bf69 100644
--- a/test/certificate.jl
+++ b/test/certificate.jl
@@ -1,3 +1,4 @@
+using LinearAlgebra
 import StarAlgebras as SA
 import MutableArithmetics as MA
 import MultivariatePolynomials as MP
@@ -5,19 +6,14 @@ import MultivariateBases as MB
 
 const SOS = SumOfSquares
 
-@testset "_merge_sorted" begin
-    @test SumOfSquares.Certificate._merge_sorted([4, 1], [3, 0]) == [4, 3, 1, 0]
-    @test SumOfSquares.Certificate._merge_sorted((4, 1), (3, 0)) == (4, 3, 1, 0)
-    @test SumOfSquares.Certificate._merge_sorted([4, 1], [3, 2]) == [4, 3, 2, 1]
-    @test SumOfSquares.Certificate._merge_sorted((4, 1), (3, 2)) == (4, 3, 2, 1)
-end
-
 @testset "with_variables" begin
     @polyvar x y z
     p = x + z
-    v = SumOfSquares.Certificate.with_variables(p, y)
+    v = SumOfSquares.Certificate.with_variables(p, FullSpace())
     @test v.inner === p
-    @test MP.variables(v.inner) == [x, z]
+    v = SumOfSquares.Certificate.with_variables(p, y)
+    @test v.inner == p
+    @test MP.variables(v.inner) == [x, y, z]
     @test MP.variables(v) == [x, y, z]
 
     @ncpolyvar a b
@@ -64,29 +60,29 @@ end
             SOS.Certificate.monomials_half_newton_polytope([x, y], uni),
         )
         @test SOS.Certificate.monomials_half_newton_polytope([x^2, y^2], uni) ==
-              [x, y]
+              [y, x]
         @test SOS.Certificate.monomials_half_newton_polytope(
             [x^2, y^2],
             Certificate.NewtonDegreeBounds(([x, y],)),
-        ) == [x, y]
+        ) == [y, x]
         @test SOS.Certificate.monomials_half_newton_polytope(
             [x^2, y^2],
             Certificate.NewtonDegreeBounds(([y, x],)),
-        ) == [x, y]
+        ) == [y, x]
         @test SOS.Certificate.monomials_half_newton_polytope(
             [x^2, x^3 * y^2, x^4 * y^4],
             uni,
-        ) == [x^2 * y^2, x, x * y, x^2, x * y^2, x^2 * y]
+        ) == MP.monomial_vector([x^2 * y^2, x, x * y, x^2, x * y^2, x^2 * y])
         @test SOS.Certificate.monomials_half_newton_polytope(
             [x^2, x^3 * y^2, x^4 * y^4],
             Certificate.NewtonFilter(uni),
-        ) == [x^2 * y^2, x]
+        ) == [x, x^2 * y^2]
     end
     @testset "Non-commutative" begin
         @test SOS.Certificate.monomials_half_newton_polytope(
             [a^4, a^3 * b, a * b * a^2, a * b * a * b],
             uni,
-        ) == [a^2, a * b]
+        ) == [a * b, a^2]
         @test SOS.Certificate.monomials_half_newton_polytope(
             [
                 a^2,
@@ -95,13 +91,13 @@ end
                 a^10 * b^20 * a^20 * b^20 * a^10,
             ],
             Certificate.NewtonFilter(uni),
-        ) == [a^10 * b^20 * a^10, a]
+        ) == [a, a^10 * b^20 * a^10]
     end
     @testset "Multipartite" begin
         # In the part [y, z], the degree is between 0 and 2
         X = [x^4, x^2 * y^2, x^2 * z^2, x^2 * y * z, y * z]
         @test SOS.Certificate.monomials_half_newton_polytope(X, uni) ==
-              [x^2, x * y, x * z, y * z, x, y, z]
+              sort([x^2, x * y, x * z, y * z, x, y, z])
         function full_test(X, Y, part1, part2)
             @test SOS.Certificate.monomials_half_newton_polytope(
                 X,
@@ -141,7 +137,7 @@ end
         @test SOS.Certificate.monomials_half_newton_polytope(
             [x^4, x^2 * y^2, x^2 * z^2, x^2 * y * z, y * z],
             Certificate.NewtonDegreeBounds(([x], [y], [z])),
-        ) == [x^2, x * y, x * z, y * z, x, y, z]
+        ) == sort([x^2, x * y, x * z, y * z, x, y, z])
     end
 end
 
@@ -172,7 +168,8 @@ function _basis_check_each(basis::SA.ExplicitBasis, basis_type)
         # `dot(X, Q * X)` which does not work (`promote_operation` calls `zero(eltype(X))`
         # which gives `Polynomial{true, Int}` which then tries to multiply a
         # `ScalarAffineFunction{Float64}` with an `Int`).
-        monos = basis.monomials
+        monos = MB.keys_as_monomials(basis)
+        # FIXME not the `monomial_vector` anymore
         @test typeof(monos) == typeof(monomial_vector(monos))
         @test issorted(monos)
     end
@@ -189,14 +186,13 @@ function _basis_check(basis, basis_type)
     end
 end
 
-function certificate_api(certificate::Certificate.AbstractIdealCertificate)
+function certificate_api(certificate::Certificate.AbstractIdealCertificate, x)
     _certificate_api(certificate)
-    @polyvar x
     poly = x + 1
     domain = @set x == 1
     a = MB.algebra_element(
         MB.sparse_coefficients(poly),
-        MB.FullBasis{MB.Monomial,MP.monomial_type(poly)}(),
+        MB.FullBasis{MB.Monomial}(x^2),
     )
     @test Certificate.reduced_polynomial(certificate, a, domain) isa
           SA.AlgebraElement
@@ -228,14 +224,16 @@ function certificate_api(certificate::Certificate.AbstractIdealCertificate)
     )
 end
 
-function certificate_api(certificate::Certificate.AbstractPreorderCertificate)
+function certificate_api(
+    certificate::Certificate.AbstractPreorderCertificate,
+    x,
+)
     _certificate_api(certificate)
-    @polyvar x
     poly = x + 1
     domain = @set x >= 1
     a = MB.algebra_element(
         MB.sparse_coefficients(poly),
-        MB.FullBasis{MB.Monomial,MP.monomial_type(poly)}(),
+        MB.FullBasis{MB.Monomial}(x),
     )
     processed = Certificate.preprocessed_domain(
         certificate,
@@ -262,10 +260,10 @@ end
     @polyvar x
     cone = SumOfSquares.SOSCone()
     B = MB.Monomial
-    full_basis = MB.FullBasis{B,MP.monomial_type(x)}()
+    full_basis = MB.FullBasis{B}(x^2)
     maxdegree = 2
     function _test(certificate::Certificate.AbstractIdealCertificate)
-        certificate_api(certificate)
+        certificate_api(certificate, x)
         mult_cert = certificate
         if mult_cert isa Certificate.Sparsity.Ideal
             mult_cert = mult_cert.certificate
@@ -278,13 +276,16 @@ end
                 Certificate.MaxDegree(cone, full_basis, full_basis, maxdegree)
         end
         preorder = Certificate.Putinar(mult_cert, certificate, maxdegree)
-        certificate_api(preorder)
+        certificate_api(preorder, x)
         sparsities = Sparsity.Pattern[Sparsity.Variable()]
         if certificate isa Certificate.MaxDegree
             push!(sparsities, Sparsity.Monomial(ChordalCompletion(), 1))
         end
         @testset "$(typeof(sparsity))" for sparsity in sparsities
-            certificate_api(Certificate.Sparsity.Preorder(sparsity, preorder))
+            certificate_api(
+                Certificate.Sparsity.Preorder(sparsity, preorder),
+                x,
+            )
         end
     end
     basis = MB.SubBasis{B}(monomial_vector([x^2, x]))
@@ -293,6 +294,8 @@ end
         Certificate.FixedBasis(cone, basis, full_basis),
         Certificate.Newton(cone, full_basis, full_basis, tuple()),
     ]
+        certificate =
+            Certificate.MaxDegree(cone, full_basis, full_basis, maxdegree)
         _test(certificate)
         _test(Certificate.Remainder(certificate))
         if certificate isa Certificate.MaxDegree
@@ -315,14 +318,13 @@ end
 
 function test_putinar_ijk(i, j, k, default::Bool, post_filter::Bool = default)
     v = @polyvar x y
-    poly = x^(2i) + y^(2j + 1)
-    domain = @set y^(2k + 1) >= 0
+    domain, poly = SA.promote_bases(@set(y^(2k + 1) >= 0), x^(2i) + y^(2j + 1))
     if default
         certificate = JuMP.moi_set(
             SOSCone(),
             MB.SubBasis{MB.Monomial}(monomials(poly)),
-            MB.FullBasis{MB.Monomial,MP.monomial_type(poly)}(),
-            MB.FullBasis{MB.Monomial,MP.monomial_type(poly)}();
+            MB.FullBasis{MB.Monomial}(poly),
+            MB.FullBasis{MB.Monomial}(poly);
             domain,
         ).certificate
     else
@@ -332,20 +334,21 @@ function test_putinar_ijk(i, j, k, default::Bool, post_filter::Bool = default)
         end
         cert = Certificate.Newton(
             SOSCone(),
-            MB.FullBasis{MB.Monomial,MP.monomial_type(x * y)}(),
-            MB.FullBasis{MB.Monomial,MP.monomial_type(x * y)}(),
+            MB.FullBasis{MB.Monomial}(x * y),
+            MB.FullBasis{MB.Monomial}(x * y),
             newton,
         )
         certificate = Certificate.Putinar(cert, cert, max(2i, 2j + 1, 2k + 1))
     end
     alg_el = MB.algebra_element(
         MB.sparse_coefficients(poly),
-        MB.FullBasis{MB.Monomial,MP.monomial_type(poly)}(),
+        MB.FullBasis{MB.Monomial}(poly),
     )
     processed = Certificate.preprocessed_domain(certificate, domain, alg_el)
     for idx in Certificate.preorder_indices(certificate, processed)
-        monos =
-            Certificate.multiplier_basis(certificate, idx, processed).monomials
+        monos = MB.keys_as_monomials(
+            Certificate.multiplier_basis(certificate, idx, processed),
+        )
         if k > j
             @test isempty(monos)
         else
diff --git a/test/gram_matrix.jl b/test/gram_matrix.jl
index 3c90ef82d..571936b41 100644
--- a/test/gram_matrix.jl
+++ b/test/gram_matrix.jl
@@ -1,16 +1,6 @@
 using LinearAlgebra, Test, SumOfSquares
 
-function _algebra_element(poly)
-    return MB.algebra_element(
-        SA.SparseCoefficients(
-            collect(MP.monomials(poly)),
-            collect(MP.coefficients(poly)),
-        ),
-        MB.FullBasis{MB.Monomial,MP.monomial_type(poly)}(),
-    )
-end
-
-_sos_dec(polys) = SOSDecomposition(_algebra_element.(polys))
+_sos_dec(polys) = SOSDecomposition(MB.algebra_element.(polys))
 
 @testset "GramMatrix tests" begin
     @testset "GramMatrix" begin
@@ -42,8 +32,8 @@ _sos_dec(polys) = SOSDecomposition(_algebra_element.(polys))
             GramMatrix([1 2; 2 4], monomial_vector([y, x])),
         )
             @test P.Q.Q == [1, 2, 4]
-            @test P.basis.monomials[1] == y
-            @test P.basis.monomials[2] == x
+            @test MB.keys_as_monomials(P.basis)[1] == y
+            @test MB.keys_as_monomials(P.basis)[2] == x
         end
         P = GramMatrix{Int}(
             (i, j) -> ((i, j) == (1, 1) ? 2 : 0),
@@ -91,14 +81,14 @@ _sos_dec(polys) = SOSDecomposition(_algebra_element.(polys))
             q = GramMatrix(5 * ones(1, 1), [x])
             r = @inferred gram_operate(/, q, 5)
             @test r.Q == ones(1, 1)
-            @test r.basis.monomials == [x]
+            @test MB.keys_as_monomials(r.basis) == [x]
             r = @inferred gram_operate(+, p, q)
             @test r.Q == [2 3; 3 7]
-            @test r.basis.monomials == [x, y]
+            @test MB.keys_as_monomials(r.basis) == [x, y]
             q = GramMatrix(5 * ones(1, 1), [y])
             r = @inferred gram_operate(+, p, q)
             @test r.Q == [7 3; 3 2]
-            @test r.basis.monomials == [x, y]
+            @test MB.keys_as_monomials(r.basis) == [x, y]
             q = GramMatrix([5.0 7; 7 9], [x * y, 1])
             r = @inferred gram_operate(+, p, q)
             @test r.Q == [
@@ -107,7 +97,7 @@ _sos_dec(polys) = SOSDecomposition(_algebra_element.(polys))
                 0 3 2 0
                 7 0 0 5
             ]
-            @test r.basis.monomials == [x * y, x, y, 1]
+            @test MB.keys_as_monomials(r.basis) == [1, y, x, x * y]
         end
         @testset "With VariableIndex" begin
             a = MOI.VariableIndex(1)
@@ -126,7 +116,7 @@ _sos_dec(polys) = SOSDecomposition(_algebra_element.(polys))
         #       ps = [1, x + y, x^2, x*y, 1 + x + x^2]
         #       P = GramMatrix(SOSDecomposition(ps))
         #       P.Q == [2 0 1 0 1; 0 1 0 0 0; 1 0 2 1 1; 0 0 1 1 0; 1 0 1 0 2]
-        #       P.basis.monomials == [x^2, x*y, x, y, 1]
+        #       MB.keys_as_monomials(P.basis) == [x^2, x*y, x, y, 1]
         #       @test P == P
         #       @test isapprox(GramMatrix(SOSDecomposition(P)), P)
         P = GramMatrix{Int}((i, j) -> i + j, [x^2, x * y, y^2])
@@ -166,15 +156,9 @@ _sos_dec(polys) = SOSDecomposition(_algebra_element.(polys))
         ps = _sos_dec([x + y, x - y])
         ps1 = _sos_dec([x])
         ps2 = _sos_dec([y])
-        M = typeof(x * y)
-        @test [ps, ps1] isa Vector{
-            SOSDecomposition{
-                MB.Algebra{MB.FullBasis{MB.Monomial,M},MB.Monomial,M},
-                Int,
-                SA.SparseCoefficients{M,Int,Vector{M},Vector{Int64}},
-                Int,
-            },
-        }
+        A = typeof(SA.parent(first(ps.ps)))
+        V = typeof(SA.coeffs(first(ps.ps)))
+        @test [ps, ps1] isa Vector{SOSDecomposition{Int,A,V,Int}}
         @test sprint(show, SOSDecompositionWithDomain(ps, [ps1, ps2], K)) ==
               "(1·y + 1·x)^2 + (-1·y + 1·x)^2 + (1·x)^2 * (1 - x^2) + (1·y)^2 * (1 - y^2)"
 
diff --git a/test/nc.jl b/test/nc.jl
new file mode 100644
index 000000000..4166e631f
--- /dev/null
+++ b/test/nc.jl
@@ -0,0 +1,185 @@
+module TestNoncommutative
+
+using Test
+using SumOfSquares
+using DynamicPolynomials
+using JuMP
+import Clarabel
+import StarAlgebras as SA
+import MultivariateBases as MB
+import MultivariatePolynomials as MP
+using MultivariateMoments: SymMatrix
+
+const nc_optimizer =
+    optimizer_with_attributes(Clarabel.Optimizer, "verbose" => false)
+
+function test_GramMatrix_with_NC_variables()
+    @ncpolyvar x y
+    # NC polynomial algebra_element and polynomial recovery
+    p = x * y + y * x + x^2
+    a = MB.algebra_element(p)
+    @test MP.polynomial(a) == p
+    # NC monomials are distinct: xy ≠ yx
+    @test x * y != y * x
+    @test length(monomials(p)) == 3
+end
+
+function test_SOS_with_NC() # (xy + x^2)^2 [BKP16 Example 2.11]
+    @ncpolyvar x y
+    p = (x * y + x^2)^2
+    model = Model(nc_optimizer)
+    cref = @constraint(model, p in SOSCone())
+    optimize!(model)
+    @test termination_status(model) == OPTIMAL
+
+    # Gram matrix and polynomial recovery
+    g = gram_matrix(cref)
+    @test polynomial(g) ≈ p atol = 1e-5
+
+    # Newton polytope should select [xy, x^2] as basis
+    basis_monos = MB.keys_as_monomials(g.basis)
+    @test Set(basis_monos) == Set([x * y, x^2])
+
+    # SOS decomposition should have 1 term
+    dec = sos_decomposition(cref, 1e-6)
+    @test length(dec.ps) == 1
+    # The decomposition q should satisfy q' * q ≈ p
+    q = polynomial(dec.ps[1])
+    @test q' * q ≈ p atol = 1e-5
+end
+
+function test_large_Newton_polytope() # [BKP16 Example 2.2]
+    @ncpolyvar x y
+    n = 5
+    p = (x + x^n * y^(2n) * x^n)^2
+    model = Model(nc_optimizer)
+    cref = @constraint(model, p in SOSCone())
+    optimize!(model)
+    @test termination_status(model) == OPTIMAL
+
+    # Newton chip should give only 2 basis monomials
+    g = gram_matrix(cref)
+    basis_monos = MB.keys_as_monomials(g.basis)
+    @test length(basis_monos) == 2
+    @test Set(basis_monos) == Set([x, x^n * y^(2n) * x^n])
+    @test polynomial(g) ≈ p atol = 1e-4
+end
+
+function test_Commutator_squared_is_SOS()
+    @ncpolyvar x y
+    comm = x * y - y * x
+    p = comm^2  # = xyxy + yxyx - xy^2x - yx^2y
+    model = Model(nc_optimizer)
+    cref = @constraint(model, p in SOSCone())
+    optimize!(model)
+    @test termination_status(model) == OPTIMAL
+    @test polynomial(gram_matrix(cref)) ≈ p atol = 1e-4
+end
+
+function test_Higher_degree_commutator_squared()
+    @ncpolyvar x y
+    p = (x^2 * y - y * x^2)^2
+    model = Model(nc_optimizer)
+    cref = @constraint(model, p in SOSCone())
+    optimize!(model)
+    @test termination_status(model) == OPTIMAL
+    dec = sos_decomposition(cref, 1e-6)
+    @test length(dec.ps) >= 1
+    @test polynomial(gram_matrix(cref)) ≈ p atol = 1e-4
+end
+
+function test_Commutator_is_not_SOS()
+    @ncpolyvar x y
+    model = Model(nc_optimizer)
+    @constraint(model, x * y - y * x in SOSCone())
+    optimize!(model)
+    @test termination_status(model) in [INFEASIBLE, INFEASIBLE_OR_UNBOUNDED]
+end
+
+function test_SOS_optimization()
+    @ncpolyvar x y
+    # max t s.t. x^2 + y^2 - t(xy + yx) is SOS
+    # Gram matrix Q = [[1, -t], [-t, 1]] on basis [x, y]
+    # PSD iff t <= 1
+    model = Model(nc_optimizer)
+    @variable(model, t)
+    @constraint(model, x^2 + y^2 - t * (x * y + y * x) in SOSCone())
+    @objective(model, Max, t)
+    optimize!(model)
+    @test termination_status(model) == OPTIMAL
+    @test value(t) ≈ 1.0 atol = 1e-4
+end
+
+function test_quartic_optimization()
+    @ncpolyvar x y
+    # max t s.t. x^4 + y^4 + xyxy + yxyx - t(x^4 + x^2y^2 + y^2x^2 + y^4) SOS
+    # This tests optimization with NC degree-4 polynomials
+    model = Model(nc_optimizer)
+    @variable(model, t)
+    p = x^4 + y^4 + x * y * x * y + y * x * y * x
+    q = x^4 + x^2 * y^2 + y^2 * x^2 + y^4
+    @constraint(model, p - t * q in SOSCone())
+    @objective(model, Max, t)
+    optimize!(model)
+    @test termination_status(model) == OPTIMAL
+    @test value(t) ≈ 0.5 atol = 1e-4
+end
+
+function test_Hermitian_SOS()
+    @ncpolyvar x y
+    p = (x + im * y) * (x - im * y)  # = x^2 + y^2 + i(yx - xy)
+    model = Model(nc_optimizer)
+    cone = NonnegPolyInnerCone{MOI.HermitianPositiveSemidefiniteConeTriangle}()
+    cref = @constraint(model, p in cone)
+    optimize!(model)
+    @test termination_status(model) == OPTIMAL
+    dec = sos_decomposition(cref, 1e-6)
+    @test length(dec.ps) == 1
+    # Verify decomposition: p* * p should give original
+    q = polynomial(dec.ps[1])
+    @test q' * q ≈ p atol = 1e-5
+end
+
+function test_NC__not_SOS()
+    @ncpolyvar x y
+    # In NC, x^2*y^2 + y^2*x^2 is NOT a sum of squares
+    p = x^2 * y^2 + y^2 * x^2
+    model = Model(nc_optimizer)
+    @constraint(model, p in SOSCone())
+    optimize!(model)
+    @test termination_status(model) in [INFEASIBLE, INFEASIBLE_OR_UNBOUNDED]
+end
+
+function test_NC_Newton_polytope()
+    @ncpolyvar a b
+    uni = Certificate.NewtonDegreeBounds(tuple())
+    # Test from certificate.jl
+    @test SumOfSquares.Certificate.monomials_half_newton_polytope(
+        [a^4, a^3 * b, a * b * a^2, a * b * a * b],
+        uni,
+    ) == [a * b, a^2]
+    # With NewtonFilter
+    @test SumOfSquares.Certificate.monomials_half_newton_polytope(
+        [
+            a^2,
+            a^10 * b^20 * a^11,
+            a^11 * b^20 * a^10,
+            a^10 * b^20 * a^20 * b^20 * a^10,
+        ],
+        Certificate.NewtonFilter(uni),
+    ) == [a, a^10 * b^20 * a^10]
+end
+
+function runtests()
+    for name in names(@__MODULE__; all = true)
+        if startswith("$name", "test_")
+            @testset "$(name)" begin
+                getfield(@__MODULE__, name)()
+            end
+        end
+    end
+end
+
+end
+
+TestNoncommutative.runtests()
diff --git a/test/runtests.jl b/test/runtests.jl
index e3d5bf0d1..9486b290b 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -30,6 +30,9 @@ include("certificate.jl")
 
 include("gram_matrix.jl")
 
+include("nc.jl")
+include("sets.jl")
+
 include("variable.jl")
 include("constraint.jl")
 include("Bridges/Bridges.jl")
diff --git a/test/sets.jl b/test/sets.jl
new file mode 100644
index 000000000..a199ddfe5
--- /dev/null
+++ b/test/sets.jl
@@ -0,0 +1,46 @@
+using Test
+using SumOfSquares
+import MultivariateBases as MB
+import MultivariatePolynomials as MP
+import StarAlgebras as SA
+import MathOptInterface as MOI
+using DynamicPolynomials
+using LinearAlgebra
+
+@testset "ScaledMonomial gram matrix" begin
+    @polyvar x y
+
+    # For gram basis [y, x] the polynomial is p = b^T Q b.
+    # With Q = I (identity), p should be the same regardless of basis.
+
+    # --- Monomial gram basis ---
+    # b = [y, x], Q = I ⟹ p = y² + 2·1·xy + x² = y² + 2xy + x²
+    # Vectorized Q = [Q₁₁, Q₁₂, Q₂₂] = [1, 1, 1] for the monomial basis
+    p_expected = y^2 + 2x * y + x^2
+
+    Q_mono = [1.0, 1.0, 1.0]
+    gb_mono = MB.SubBasis{MB.Monomial}([y, x])
+    gram_mono = SumOfSquares.build_gram_matrix(
+        Q_mono,
+        gb_mono,
+        MOI.PositiveSemidefiniteConeTriangle,
+        Float64,
+    )
+    p_mono = MP.polynomial(gram_mono)
+    @test p_mono ≈ p_expected
+
+    # --- ScaledMonomial gram basis ---
+    # With the correct MStruct, ŷ·x̂ has coefficient 1/√2 in the algebra,
+    # so b̃^T Q̃ b̃ = Q̃₁₁ ŷ² + 2·(1/√2)·Q̃₁₂ ŝ(xy) + Q̃₂₂ x̂²
+    # For Q̃ = I: p = ŷ² + √2·ŝ(xy) + x̂² = y² + √2·√2·xy + x² = y² + 2xy + x²
+    Q_scaled = [1.0, 1.0, 1.0]
+    gb_scaled = MB.SubBasis{MB.ScaledMonomial}([y, x])
+    gram_scaled = SumOfSquares.build_gram_matrix(
+        Q_scaled,
+        gb_scaled,
+        MOI.PositiveSemidefiniteConeTriangle,
+        Float64,
+    )
+    p_scaled = MP.polynomial(gram_scaled)
+    @test p_scaled ≈ p_expected
+end
diff --git a/test/sparsity.jl b/test/sparsity.jl
index 1c9f0043d..e84621b43 100644
--- a/test/sparsity.jl
+++ b/test/sparsity.jl
@@ -4,9 +4,13 @@ using SumOfSquares
 import MultivariateBases as MB
 
 function _algebra_element(monos)
+    basis = MB.FullBasis{MB.Monomial}(MP.variables(monos))
     return MB.algebra_element(
-        SA.SparseCoefficients(monos, ones(length(monos))),
-        MB.FullBasis{MB.Monomial,eltype(monos)}(),
+        SA.SparseCoefficients(
+            [basis.inverse_map(m) for m in monos],
+            ones(length(monos)),
+        ),
+        basis,
     )
 end
 
@@ -36,7 +40,10 @@ function xor_complement_test()
 end
 
 function set_monos(bases::Vector{<:MB.SubBasis})
-    return Set([basis.monomials for basis in bases])
+    return Set([MB.keys_as_monomials(basis) for basis in bases])
+end
+function set_monos(basis::MB.SubBasis)
+    return Set([MB.keys_as_monomials(basis)])
 end
 
 """
@@ -50,15 +57,15 @@ function wml19()
     @polyvar x[1:3]
     certificate = Certificate.Newton(
         SOSCone(),
-        MB.FullBasis{MB.Monomial,typeof(prod(x))}(),
-        MB.FullBasis{MB.Monomial,typeof(prod(x))}(),
+        MB.FullBasis{MB.Monomial}(x),
+        MB.FullBasis{MB.Monomial}(x),
         tuple(),
     )
     @testset "Example 4.2" begin
         f = 1 + x[1]^4 + x[2]^4 + x[3]^4 + prod(x) + x[2]
         alg_el = MB.algebra_element(
             MB.sparse_coefficients(f),
-            MB.FullBasis{MB.Monomial,MP.monomial_type(f)}(),
+            MB.FullBasis{MB.Monomial}(f),
         )
         with_var = SumOfSquares.Certificate.WithVariables(alg_el, x)
         expected_1_false = Set(
@@ -150,22 +157,22 @@ function wml19()
                                                                     [
         Certificate.MaxDegree(
             SOSCone(),
-            MB.FullBasis{MB.Monomial,typeof(prod(x[1:2]))}(),
-            MB.FullBasis{MB.Monomial,typeof(prod(x[1:2]))}(),
+            MB.FullBasis{MB.Monomial}(prod(x[1:2])),
+            MB.FullBasis{MB.Monomial}(prod(x[1:2])),
             4,
         ),
         Certificate.Newton(
             SOSCone(),
-            MB.FullBasis{MB.Monomial,typeof(prod(x[1:2]))}(),
-            MB.FullBasis{MB.Monomial,typeof(prod(x[1:2]))}(),
+            MB.FullBasis{MB.Monomial}(prod(x[1:2])),
+            MB.FullBasis{MB.Monomial}(prod(x[1:2])),
             tuple(),
         ),
     ]
         preorder_certificate = Certificate.Putinar(
             Certificate.MaxDegree(
                 SOSCone(),
-                MB.FullBasis{MB.Monomial,typeof(prod(x[1:2]))}(),
-                MB.FullBasis{MB.Monomial,typeof(prod(x[1:2]))}(),
+                MB.FullBasis{MB.Monomial}(x[1:2]),
+                MB.FullBasis{MB.Monomial}(x[1:2]),
                 4,
             ),
             ideal_certificate,
@@ -174,7 +181,7 @@ function wml19()
         f = x[1]^4 + x[2]^4 + x[1] * x[2]
         alg_el = MB.algebra_element(
             MB.sparse_coefficients(f),
-            MB.FullBasis{MB.Monomial,MP.monomial_type(f)}(),
+            MB.FullBasis{MB.Monomial}(f),
         )
         K = @set 1 - 2x[1]^2 - x[2]^2 >= 0
         @testset "$(nameof(typeof(completion))) $k $use_all_monomials" for completion in
@@ -241,9 +248,9 @@ function wml19()
             x[1] * x[2]^2 - 3x[1]^2 * x[2]^2
         alg_el = MB.algebra_element(
             MB.sparse_coefficients(f),
-            MB.FullBasis{MB.Monomial,MP.monomial_type(f)}(),
+            MB.FullBasis{MB.Monomial}(f),
         )
-        with_var = SumOfSquares.Certificate.WithVariables(alg_el, x)
+        with_var = SumOfSquares.Certificate.WithVariables(alg_el, x[1:2])
         basis = MB.SubBasis{MB.Monomial}(MP.monomials(f))
         @testset "$(nameof(typeof(completion))) $k $use_all_monomials" for completion in
                                                                            [
@@ -286,13 +293,10 @@ function wml19()
                 certificate,
             ),
         ) == Set(
-            monomial_vector.([[
-                x[1]^2 * x[2]^2,
-                x[1] * x[2]^2,
-                1,
-                x[1]^2 * x[2],
-                x[1] * x[2],
-            ],]),
+            monomial_vector.([
+                [x[1]^2 * x[2]^2, x[1] * x[2]^2, 1],
+                [x[1]^2 * x[2], x[1] * x[2]],
+            ]),
         )
     end
 end
@@ -305,18 +309,18 @@ Examples of [MWL19].
 [L09] Lofberg, Johan. "Pre-and post-processing sum-of-squares programs in practice." IEEE transactions on automatic control 54.5 (2009): 1007-1011.
 """
 function l09()
-    @polyvar x[1:3]
-    certificate = Certificate.Newton(
-        SOSCone(),
-        MB.FullBasis{MB.Monomial,typeof(prod(x[1:2]))}(),
-        MB.FullBasis{MB.Monomial,typeof(prod(x[1:2]))}(),
-        tuple(),
-    )
     @testset "Example 1 and 2" begin
+        @polyvar x[1:2]
+        certificate = Certificate.Newton(
+            SOSCone(),
+            MB.FullBasis{MB.Monomial}(x),
+            MB.FullBasis{MB.Monomial}(x),
+            tuple(),
+        )
         f = 1 + x[1]^4 * x[2]^2 + x[1]^2 * x[2]^4
         alg_el = MB.algebra_element(
             MB.sparse_coefficients(f),
-            MB.FullBasis{MB.Monomial,MP.monomial_type(f)}(),
+            MB.FullBasis{MB.Monomial}(f),
         )
         with_var = SumOfSquares.Certificate.WithVariables(alg_el, x)
         newt = Certificate.NewtonDegreeBounds(tuple())
@@ -356,7 +360,8 @@ function l09()
         expected = Set(
             monomial_vector.([
                 [x[1]^2 * x[2]],
-                [x[1] * x[2]^2, constant_monomial(x[1] * x[2])],
+                [x[1] * x[2]^2],
+                [constant_monomial(x[1] * x[2])],
             ]),
         )
         @test set_monos(
@@ -368,10 +373,17 @@ function l09()
         ) == expected
     end
     @testset "Example 3 and 4" begin
+        @polyvar x[1:3]
+        certificate = Certificate.Newton(
+            SOSCone(),
+            MB.FullBasis{MB.Monomial}(x),
+            MB.FullBasis{MB.Monomial}(x),
+            tuple(),
+        )
         f = 1 + x[1]^4 + x[1] * x[2] + x[2]^4 + x[3]^2
         alg_el = MB.algebra_element(
             MB.sparse_coefficients(f),
-            MB.FullBasis{MB.Monomial,MP.monomial_type(f)}(),
+            MB.FullBasis{MB.Monomial}(f),
         )
         with_var = SumOfSquares.Certificate.WithVariables(alg_el, x)
         @testset "$k $use_all_monomials" for k in 0:2,
@@ -458,16 +470,16 @@ function square_domain(ideal_certificate, d)
     @polyvar x y
     mult_cert = Certificate.MaxDegree(
         SOSCone(),
-        MB.FullBasis{MB.Monomial,typeof(x * y)}(),
-        MB.FullBasis{MB.Monomial,typeof(x * y)}(),
+        MB.FullBasis{MB.Monomial}(x * y),
+        MB.FullBasis{MB.Monomial}(x * y),
         d,
     )
     preorder_certificate =
-        Certificate.Putinar(mult_cert, ideal_certificate(typeof(x * y)), d)
+        Certificate.Putinar(mult_cert, ideal_certificate(x * y), d)
     f = x^2 * y^4 + x^4 * y^2 - 3 * x^2 * y * 2 + 1
     alg_el = MB.algebra_element(
         MB.sparse_coefficients(f),
-        MB.FullBasis{MB.Monomial,MP.monomial_type(f)}(),
+        MB.FullBasis{MB.Monomial}(f),
     )
     K = @set(1 - x^2 >= 0 && 1 - y^2 >= 0)
     @testset "Square domain $k $use_all_monomials" for k in 0:4,
@@ -577,14 +589,14 @@ function sum_square(n)
         @polyvar x[1:(2n)]
         certificate = Certificate.Newton(
             SOSCone(),
-            MB.FullBasis{MB.Monomial,typeof(prod(x))}(),
-            MB.FullBasis{MB.Monomial,typeof(prod(x))}(),
+            MB.FullBasis{MB.Monomial}(x),
+            MB.FullBasis{MB.Monomial}(x),
             tuple(),
         )
         f = sum((x[1:2:(2n-1)] .- x[2:2:(2n)]) .^ 2)
         alg_el = MB.algebra_element(
             MB.sparse_coefficients(f),
-            MB.FullBasis{MB.Monomial,MP.monomial_type(f)}(),
+            MB.FullBasis{MB.Monomial}(f),
         )
         expected = Set(
             monomial_vector.([
@@ -597,8 +609,8 @@ function sum_square(n)
                 Sparsity.Variable(),
                 Certificate.MaxDegree(
                     SOSCone(),
-                    MB.FullBasis{MB.Monomial,typeof(prod(x))}(),
-                    MB.FullBasis{MB.Monomial,typeof(prod(x))}(),
+                    MB.FullBasis{MB.Monomial}(x),
+                    MB.FullBasis{MB.Monomial}(x),
                     2,
                 ),
             ),
@@ -620,8 +632,8 @@ function drop_monomials()
         alg_el = _algebra_element([x^2])
         certificate = Certificate.MaxDegree(
             SOSCone(),
-            MB.FullBasis{MB.Monomial,typeof(x^2)}(),
-            MB.FullBasis{MB.Monomial,typeof(x^2)}(),
+            MB.FullBasis{MB.Monomial}(x^2),
+            MB.FullBasis{MB.Monomial}(x^2),
             2,
         )
         @testset "$k $use_all_monomials" for k in 0:2,
@@ -649,22 +661,22 @@ function drop_monomials()
                                                             [
             Certificate.MaxDegree(
                 SOSCone(),
-                MB.FullBasis{MB.Monomial,typeof(x^2)}(),
-                MB.FullBasis{MB.Monomial,typeof(x^2)}(),
+                MB.FullBasis{MB.Monomial}(x^2),
+                MB.FullBasis{MB.Monomial}(x^2),
                 4,
             ),
             Certificate.Newton(
                 SOSCone(),
-                MB.FullBasis{MB.Monomial,typeof(x^2)}(),
-                MB.FullBasis{MB.Monomial,typeof(x^2)}(),
+                MB.FullBasis{MB.Monomial}(x^2),
+                MB.FullBasis{MB.Monomial}(x^2),
                 tuple(),
             ),
         ]
             preorder_certificate = Certificate.Putinar(
                 Certificate.MaxDegree(
                     SOSCone(),
-                    MB.FullBasis{MB.Monomial,typeof(x^2)}(),
-                    MB.FullBasis{MB.Monomial,typeof(x^2)}(),
+                    MB.FullBasis{MB.Monomial}(x^2),
+                    MB.FullBasis{MB.Monomial}(x^2),
                     3,
                 ),
                 ideal_certificate,
@@ -723,19 +735,19 @@ end
     wml19()
     l09()
     square_domain(
-        M -> Certificate.MaxDegree(
+        m -> Certificate.MaxDegree(
             SOSCone(),
-            MB.FullBasis{MB.Monomial,M}(),
-            MB.FullBasis{MB.Monomial,M}(),
+            MB.FullBasis{MB.Monomial}(m),
+            MB.FullBasis{MB.Monomial}(m),
             6,
         ),
         6,
     )
     square_domain(
-        M -> Certificate.Newton(
+        m -> Certificate.Newton(
             SOSCone(),
-            MB.FullBasis{MB.Monomial,M}(),
-            MB.FullBasis{MB.Monomial,M}(),
+            MB.FullBasis{MB.Monomial}(m),
+            MB.FullBasis{MB.Monomial}(m),
             tuple(),
         ),
         6,
diff --git a/test/variable.jl b/test/variable.jl
index a9be83cc0..bb7f40b35 100644
--- a/test/variable.jl
+++ b/test/variable.jl
@@ -25,8 +25,8 @@ end
     X = monomials([x, y], 0:2)
     for cone in [SOSPoly(X), SDSOSPoly(X), DSOSPoly(X)]
         p = @variable(model, [1:2], cone)
-        @test p[1].basis.monomials == X
-        @test p[2].basis.monomials == X
+        @test MB.keys_as_monomials(p[1].basis) == X
+        @test MB.keys_as_monomials(p[2].basis) == X
         @test eltype(p) <: GramMatrix{JuMP.VariableRef}
     end
 end