using DynamicPolynomials
@polyvar x y
p = x^3 - x^2 + 2x*y -y^2 + y^3
using SemialgebraicSets
S = @set x >= 0 && y >= 0 && x + y >= 1
p(x=>1, y=>0), p(x=>1//2, y=>1//2), p(x=>0, y=>1)

using JuMP
using Ipopt
model = Model(optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0))
@variable(model, a >= 0)
@variable(model, b >= 0)
@constraint(model, a + b >= 1)
@NLobjective(model, Min, a^3 - a^2 + 2a*b - b^2 + b^3)
optimize!(model)

@show termination_status(model)
@show value(a)
@show value(b)
@show objective_value(model)

using JuMP
using Ipopt
model = Model(optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0))
@variable(model, a >= 0)
@variable(model, b >= 0)
@constraint(model, a + b >= 1)
peval(a, b) = p(x=>a, y=>b)
register(model, :peval, 2, peval, autodiff=true)
@NLobjective(model, Min, peval(a, b))
optimize!(model)

@show termination_status(model)
@show value(a)
@show value(b)
@show objective_value(model)

using SumOfSquares

using CSDP
optimizer = optimizer_with_attributes(CSDP.Optimizer, "printlevel" => 0);

using MosekTools
optimizer = optimizer_with_attributes(Mosek.Optimizer, "QUIET" => true);

model = SOSModel(optimizer)
@variable(model, α)
@objective(model, Max, α)
@constraint(model, c3, p >= α, domain = S)

optimize!(model)
@show termination_status(model)
@show objective_value(model)

using MultivariateMoments
ν3 = moment_matrix(c3)
extractatoms(ν3, 1e-3) # Returns nothing as the dual is not atomic

model = SOSModel(optimizer)
@variable(model, α)
@objective(model, Max, α)
@constraint(model, c5, p >= α, domain = S, maxdegree = 5)

optimize!(model)
@show termination_status(model)
@show objective_value(model)

using MultivariateMoments
ν5 = moment_matrix(c5)
extractatoms(ν5, 1e-3)

ν3 = moment_matrix(c3)
MultivariateMoments.computesupport!(ν3, 1e-3)

ν5 = moment_matrix(c5)
MultivariateMoments.computesupport!(ν5, 1e-3)

SemialgebraicSets.computegröbnerbasis!(ideal(ν5.support))
ν5.support