using DynamicPolynomials
@ncpolyvar x y
p = (x * y + x^2)^2

import CSDP
using JuMP
optimizer_constructor = optimizer_with_attributes(CSDP.Optimizer, MOI.Silent() => true)

using SumOfSquares
model = Model(optimizer_constructor)
con_ref = @constraint(model, p in SOSCone())

optimize!(model)

certificate_basis(con_ref)

gram_matrix(con_ref).Q

sos_decomposition(con_ref)

sos_decomposition(con_ref, 1e-6)

using DynamicPolynomials
@ncpolyvar x y
n = 10
p = (x + x^n * y^(2n) * x^n)^2

using SumOfSquares
model = Model(optimizer_constructor)
con_ref = @constraint(model, p in SOSCone())

optimize!(model)

certificate_basis(con_ref)

gram_matrix(con_ref).Q

sos_decomposition(con_ref, 1e-6)