{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "# Certificate"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "**Contributed by**: Benoît Legat"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Introduction"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Consider the polynomial optimization problem (a variation from [L09, Example 2.2]) of\n",
    "minimizing the polynomial $x^3 - x^2 + 2xy - y^2 + y^3$\n",
    "over the polyhedron defined by the inequalities $x \\ge 0, y \\ge 0$ and $x + y \\geq 1$."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "[L09] Lasserre, J. B.\n",
    "*Moments, positive polynomials and their applications*.\n",
    "World Scientific, **2009**."
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "Basic semialgebraic Set defined by no equality\n3 inequalities\n x ≥ 0\n y ≥ 0\n x^2 + y^2 - 2 ≥ 0\n"
     },
     "metadata": {},
     "execution_count": 1
    }
   ],
   "cell_type": "code",
   "source": [
    "using DynamicPolynomials\n",
    "@polyvar x y\n",
    "p = x^3 - x^2 + 2x*y -y^2 + y^3 + x^3 * y\n",
    "using SumOfSquares\n",
    "S = @set x >= 0 && y >= 0 && x^2 + y^2 >= 2"
   ],
   "metadata": {},
   "execution_count": 1
  },
  {
   "cell_type": "markdown",
   "source": [
    "We will now see how to find the optimal solution using Sum of Squares Programming.\n",
    "We first need to pick an SDP solver, see [here](https://jump.dev/JuMP.jl/v0.21.6/installation/#Supported-solvers) for a list of the available choices."
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "MathOptInterface.OptimizerWithAttributes(CSDP.Optimizer, Pair{MathOptInterface.AbstractOptimizerAttribute, Any}[MathOptInterface.Silent() => true])"
     },
     "metadata": {},
     "execution_count": 2
    }
   ],
   "cell_type": "code",
   "source": [
    "import CSDP\n",
    "solver = optimizer_with_attributes(CSDP.Optimizer, MOI.Silent() => true)"
   ],
   "metadata": {},
   "execution_count": 2
  },
  {
   "cell_type": "markdown",
   "source": [
    "A Sum-of-Squares certificate that $p \\ge \\alpha$ over the domain `S`, ensures that $\\alpha$ is a lower bound to the polynomial optimization problem.\n",
    "The following program searches for the largest upper bound and finds zero."
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "termination_status(model) = MathOptInterface.INFEASIBLE\n",
      "objective_value(model) = -207.16544739193887\n"
     ]
    },
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "-207.16544739193887"
     },
     "metadata": {},
     "execution_count": 3
    }
   ],
   "cell_type": "code",
   "source": [
    "model = SOSModel(solver)\n",
    "@variable(model, α)\n",
    "@objective(model, Max, α)\n",
    "@constraint(model, c, p >= α, domain = S)\n",
    "optimize!(model)\n",
    "@show termination_status(model)\n",
    "@show objective_value(model)"
   ],
   "metadata": {},
   "execution_count": 3
  },
  {
   "cell_type": "markdown",
   "source": [
    "We now define the Schmüdgen's certificate:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using MultivariateBases\n",
    "const MB = MultivariateBases\n",
    "const SOS = SumOfSquares\n",
    "const SOSC = SOS.Certificate\n",
    "struct Schmüdgen{IC <: SOSC.AbstractIdealCertificate, CT <: SOS.SOSLikeCone, BT <: SOS.AbstractPolynomialBasis} <: SOSC.AbstractPreorderCertificate\n",
    "    ideal_certificate::IC\n",
    "    cone::CT\n",
    "    basis::Type{BT}\n",
    "    maxdegree::Int\n",
    "end\n",
    "\n",
    "SOSC.cone(certificate::Schmüdgen) = certificate.cone\n",
    "\n",
    "function SOSC.preprocessed_domain(::Schmüdgen, domain::BasicSemialgebraicSet, p)\n",
    "    return SOSC.with_variables(domain, p)\n",
    "end\n",
    "\n",
    "function SOSC.preorder_indices(::Schmüdgen, domain::SOSC.WithVariables)\n",
    "    n = length(domain.inner.p)\n",
    "    if n >= Sys.WORD_SIZE\n",
    "        error(\"There are $(2^n - 1) products in Schmüdgen's certificate, they cannot even be indexed with `$Int`.\")\n",
    "    end\n",
    "    return map(SOSC.PreorderIndex, 1:(2^n-1))\n",
    "end\n",
    "\n",
    "function SOSC.multiplier_basis(certificate::Schmüdgen, index::SOSC.PreorderIndex, domain::SOSC.WithVariables)\n",
    "    q = SOSC.generator(certificate, index, domain)\n",
    "    return SOSC.maxdegree_gram_basis(certificate.basis, variables(domain), SOSC.multiplier_maxdegree(certificate.maxdegree, q))\n",
    "end\n",
    "function SOSC.multiplier_basis_type(::Type{Schmüdgen{IC, CT, BT}}) where {IC, CT, BT}\n",
    "    return BT\n",
    "end\n",
    "\n",
    "function SOSC.generator(::Schmüdgen, index::SOSC.PreorderIndex, domain::SOSC.WithVariables)\n",
    "    I = [i for i in eachindex(domain.inner.p) if !iszero(index.value & (1 << (i - 1)))]\n",
    "    return prod([domain.inner.p[i] for i in eachindex(domain.inner.p) if !iszero(index.value & (1 << (i - 1)))])\n",
    "end\n",
    "\n",
    "SOSC.ideal_certificate(certificate::Schmüdgen) = certificate.ideal_certificate\n",
    "SOSC.ideal_certificate(::Type{<:Schmüdgen{IC}}) where {IC} = IC\n",
    "\n",
    "SOS.matrix_cone_type(::Type{<:Schmüdgen{IC, CT}}) where {IC, CT} = SOS.matrix_cone_type(CT)"
   ],
   "metadata": {},
   "execution_count": 4
  },
  {
   "cell_type": "markdown",
   "source": [
    "Let's try again with our the Schmüdgen certificate:"
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "termination_status(model) = MathOptInterface.OPTIMAL\n",
      "objective_value(model) = 0.8284271199674382\n"
     ]
    },
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "0.8284271199674382"
     },
     "metadata": {},
     "execution_count": 5
    }
   ],
   "cell_type": "code",
   "source": [
    "model = SOSModel(solver)\n",
    "@variable(model, α)\n",
    "@objective(model, Max, α)\n",
    "ideal_certificate = SOSC.Newton(SOSCone(), MB.MonomialBasis, tuple())\n",
    "certificate = Schmüdgen(ideal_certificate, SOSCone(), MB.MonomialBasis, maxdegree(p))\n",
    "@constraint(model, c, p >= α, domain = S, certificate = certificate)\n",
    "optimize!(model)\n",
    "@show termination_status(model)\n",
    "@show objective_value(model)"
   ],
   "metadata": {},
   "execution_count": 5
  },
  {
   "cell_type": "markdown",
   "source": [
    "---\n",
    "\n",
    "*This notebook was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*"
   ],
   "metadata": {}
  }
 ],
 "nbformat_minor": 3,
 "metadata": {
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.7.2"
  },
  "kernelspec": {
   "name": "julia-1.7",
   "display_name": "Julia 1.7.2",
   "language": "julia"
  }
 },
 "nbformat": 4
}