using DynamicPolynomials
@polyvar x y
motzkin = x^4*y^2 + x^2*y^4 + 1 - 3x^2*y^2

using CSDP
solver = CSDPSolver();

using Mosek
solver = MosekSolver();

using SumOfSquares
using JuMP
m = SOSModel(solver = solver)
@constraint m motzkin >= 0 # We constraint `motzkin` to be a sum of squares
solve(m) # Returns the status `:Infeasible`

using SumOfSquares
using JuMP
m = SOSModel(solver = solver)
@constraint m (x^2 + y^2) * motzkin >= 0 # We constraint the `(x^2 + y^2) * motzkin` to be a sum of squares
solve(m) # Returns the status `:Optimal` which means that it is feasible

using SumOfSquares
using PolyJuMP
using JuMP
m = SOSModel(solver = solver)
X = monomials([x, y], 0:2)
# We create a quadratic polynomial that is not necessarily a sum of squares
# since this is implied by the next constraint: `deno >= 1`
@variable m deno Poly(X)
# We want the denominator polynomial to be strictly positive,
# this prevents the trivial solution deno = 0 for instance.
@constraint m deno >= 1
@constraint m deno * motzkin >= 0
solve(m)

getvalue(deno)

using RecipesBase
using MultivariatePolynomials
@recipe function f(x::AbstractVector, y::AbstractVector, p::Polynomial)
    x, y, (x, y) -> p(variables(p) => [x, y])
end

using Plots
plot(linspace(-2, 2, 100), linspace(-2, 2, 100), motzkin, st = [:surface], seriescolor=:heat, colorbar=:none, clims = (-10, 80))