{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first explicit example of nonnegative polynomial that is not a sum of squares was found by Motzkin in 1967. By the [Arithmetic-geometric mean](https://en.wikipedia.org/wiki/Arithmetic%E2%80%93geometric_mean),\n",
    "$$ \\frac{x^4y^2 + x^2y^4 + 1}{3} \\ge \\sqrt[3]{x^4y^2 \\cdot x^2y^4 \\cdot 1} = x^2y^2 $$\n",
    "hence\n",
    "$$ x^4y^2 + x^2y^4 + 1 - 3x^2y^2 \\ge 0. $$\n",
    "The code belows construct the Motzkin polynomial using [DynamicPolynomials](https://github.com/JuliaAlgebra/DynamicPolynomials.jl)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "x^4y^2 + x^2y^4 - 3x^2y^2 + 1"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using DynamicPolynomials\n",
    "@polyvar x y\n",
    "motzkin = x^4*y^2 + x^2*y^4 + 1 - 3x^2*y^2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Motzkin polynomial is nonnegative but is not a sum of squares as we can verify numerically as follows.\n",
    "We first need to pick an SDP solver, see [here](http://www.juliaopt.org/) for a list of the available choices."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using CSDP\n",
    "solver = CSDPSolver();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Mosek\n",
    "solver = MosekSolver(LOG=0);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\u001b[1m\u001b[33mWarning: \u001b[39m\u001b[22m\u001b[33mNot solved to optimality, status: Infeasible\u001b[39m\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       ":Infeasible"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using SumOfSquares\n",
    "using JuMP\n",
    "m = SOSModel(solver = solver)\n",
    "@constraint m motzkin >= 0 # We constraint `motzkin` to be a sum of squares\n",
    "solve(m) # Returns the status `:Infeasible`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Even if the Motzkin polynomial is not a sum of squares, it can still be certified to be nonnegative using sums of squares.\n",
    "Indeed a polynomial is certified to be nonnegative if it is equal to a fraction of sums of squares.\n",
    "The Motzkin polynomial is equal to a fraction of sums of squares whose denominator is $x^2 + y^2$.\n",
    "This can be verified numerically as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       ":Optimal"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using SumOfSquares\n",
    "using JuMP\n",
    "m = SOSModel(solver = solver)\n",
    "@constraint m (x^2 + y^2) * motzkin >= 0 # We constraint the `(x^2 + y^2) * motzkin` to be a sum of squares\n",
    "solve(m) # Returns the status `:Optimal` which means that it is feasible"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One may consider ourself lucky to have had the intuition that $x^2 + y^2$ would work as denominator.\n",
    "In fact, the search for the denominator can be carried out in parallel to the search of the numerator.\n",
    "In the example below, we search for a denominator with monomials of degrees from 0 to 2.\n",
    "If none is found, we can increase the maximum degree 2 to 4, 6, 8, ...\n",
    "This gives a hierarchy of programs to try in order to certify the nonnegativity of a polynomial by identifying it with a fraction of sum of squares polynomials.\n",
    "In the case of the Motzkin polynomial we now that degree 2 is enough since $x^2 + y^2$ works."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       ":Optimal"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using SumOfSquares\n",
    "using JuMP\n",
    "using MultivariatePolynomials\n",
    "m = SOSModel(solver = solver)\n",
    "X = monomials([x, y], 0:2)\n",
    "# We create a quadratic polynomial that is not necessarily a sum of squares\n",
    "# since this is implied by the next constraint: `deno >= 1`\n",
    "@variable m deno Poly(X)\n",
    "# We want the denominator polynomial to be strictly positive,\n",
    "# this prevents the trivial solution deno = 0 for instance.\n",
    "@constraint m deno >= 1\n",
    "@constraint m deno * motzkin >= 0\n",
    "solve(m)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can check the denominator found by the program using `JuMP.getvalue`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.8994524919313149x^2 - 8.417376223825856e-11xy + 0.8994524979159756y^2 + 6.987367755126592e-16x - 1.0014392160847468e-15y + 1.9999999943329119"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "getvalue(deno)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because a picture is worth a thousand words let's plot the beast.\n",
    "We can easily extend `Plots` by adding a recipe to plot bivariate polynomials."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using RecipesBase\n",
    "using MultivariatePolynomials\n",
    "@recipe function f(x::AbstractVector, y::AbstractVector, p::Polynomial)\n",
    "    x, y, (x, y) -> p(variables(p) => [x, y])\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Plots\n",
    "plot(linspace(-2, 2, 100), linspace(-2, 2, 100), motzkin, st = [:surface], seriescolor=:heat, colorbar=:none, clims = (-10, 80))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.6.4",
   "language": "julia",
   "name": "julia-0.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}