using DynamicPolynomials
@polyvar x
P = [x^2 - 2x + 2 x
            x     x^2]

using SumOfSquares

using CSDP
factory = with_optimizer(CSDP.Optimizer, printlevel=0);

using MosekTools
factory = with_optimizer(Mosek.Optimizer, QUIET=true);

model = SOSModel(factory)
mat_cref = @constraint(model, P in PSDCone())
optimize!(model)
termination_status(model)

certificate_monomials(mat_cref)

@polyvar y[1:2]
p = vec(y)' * P * vec(y)

model = SOSModel(factory)
cref = @constraint(model, p in SOSCone())

optimize!(model)
termination_status(model)

certificate_monomials(cref)