{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "# Getting started"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "**Adapted from**: SOSTOOLS' SOSDEMO1 (See Section 4.1 of [SOSTOOLS User's Manual](http://sysos.eng.ox.ac.uk/sostools/sostools.pdf)) and Example 2.4 of [PJ08]\n",
    "\n",
    "P. Parrilo and A. Jadbabaie\n",
    "*Approximation of the joint spectral radius using sum of squares*.\n",
    "Linear Algebra and its Applications, Elsevier (2008), 428, 2385-2402"
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "5y⁴ - x²y² + 2x³y + 2x⁴",
      "text/latex": "$$ 5y^{4} - x^{2}y^{2} + 2x^{3}y + 2x^{4} $$"
     },
     "metadata": {},
     "execution_count": 1
    }
   ],
   "cell_type": "code",
   "source": [
    "using DynamicPolynomials\n",
    "@polyvar x y\n",
    "p = 2*x^4 + 2*x^3*y - x^2*y^2 + 5*y^4"
   ],
   "metadata": {},
   "execution_count": 1
  },
  {
   "cell_type": "markdown",
   "source": [
    "We need to pick an SDP solver, see [here](https://jump.dev/JuMP.jl/v1.12/installation/#Supported-solvers) for a list of the available choices.\n",
    "We use `SOSModel` instead of `Model` to be able to use the `>=` syntax for Sum-of-Squares constraints."
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "FEASIBLE_POINT::ResultStatusCode = 1"
     },
     "metadata": {},
     "execution_count": 2
    }
   ],
   "cell_type": "code",
   "source": [
    "using SumOfSquares\n",
    "import CSDP\n",
    "solver = optimizer_with_attributes(CSDP.Optimizer, MOI.Silent() => true)\n",
    "model = SOSModel(solver)\n",
    "con_ref = @constraint(model, p >= 0)\n",
    "optimize!(model)\n",
    "primal_status(model)"
   ],
   "metadata": {},
   "execution_count": 2
  },
  {
   "cell_type": "markdown",
   "source": [
    "We see above that the solver found a feasible solution.\n",
    "We now inspect this solution:"
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "GramMatrix with row/column basis:\n MonomialBasis([y^2, x*y, x^2])\nAnd entries in a 3×3 SymMatrix{Float64}:\n  5.0                     4.4408920985006264e-17  -2.9518518518518517\n  4.4408920985006264e-17  4.903703703703703        1.0\n -2.9518518518518517      1.0                      2.0"
     },
     "metadata": {},
     "execution_count": 3
    }
   ],
   "cell_type": "code",
   "source": [
    "q = gram_matrix(con_ref)"
   ],
   "metadata": {},
   "execution_count": 3
  },
  {
   "cell_type": "markdown",
   "source": [
    "We can get the SOS decomposition from the gram matrix as follows:"
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "(-2.1236350925148186*y^2 + 0.6856027338723613*x*y + 1.4043287476933908*x^2)^2 + (-0.6931517747987983*y^2 - 2.105348694339436*x*y - 0.020343251432300993*x^2)^2 + (0.09856272587979271*y^2 - 0.034050994900062574*x*y + 0.16567112157245342*x^2)^2"
     },
     "metadata": {},
     "execution_count": 4
    }
   ],
   "cell_type": "code",
   "source": [
    "sosdec = SOSDecomposition(q)"
   ],
   "metadata": {},
   "execution_count": 4
  },
  {
   "cell_type": "markdown",
   "source": [
    "We now seek for the SOS decomposition of the following polynomial:"
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "y² + x² - xy² + 4x⁴y⁶",
      "text/latex": "$$ y^{2} + x^{2} - xy^{2} + 4x^{4}y^{6} $$"
     },
     "metadata": {},
     "execution_count": 5
    }
   ],
   "cell_type": "code",
   "source": [
    "p = 4*x^4*y^6 + x^2 - x*y^2 + y^2"
   ],
   "metadata": {},
   "execution_count": 5
  },
  {
   "cell_type": "markdown",
   "source": [
    "We build the same model as previously with this new polynomial.\n",
    "Here we can use `Model` instead of `SOSModel` as we explicitly constrain\n",
    "`p` to belong to the SOS cone with `p in SOSCone()`."
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "FEASIBLE_POINT::ResultStatusCode = 1"
     },
     "metadata": {},
     "execution_count": 6
    }
   ],
   "cell_type": "code",
   "source": [
    "model = Model(solver)\n",
    "con_ref = @constraint(model, p in SOSCone())\n",
    "optimize!(model)\n",
    "primal_status(model)"
   ],
   "metadata": {},
   "execution_count": 6
  },
  {
   "cell_type": "markdown",
   "source": [
    "We can query the SOS decomposition directly from the constraint reference\n",
    "as follows:"
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "(-0.8633544111907671*y + 0.1929775871076749*x*y + 2.8118914006146728e-18*x*y^2 + 1.9748040632461146*x^2*y^3)^2 + (-1.0044005280851495e-16*y + 0.7962855107812749*x - 3.604670114761173e-17*x*y - 1.7452642425911866*x*y^2 - 5.2795306922881863e-17*x^2*y^3)^2 + (0.24383463418983628*y - 1.7639164207923265e-16*x - 1.548639212897575*x*y + 3.447822046557963e-17*x*y^2 + 0.25793362243589024*x^2*y^3)^2 + (-5.040390840687567e-17*y + 0.604920974441955*x + 4.554029379037467e-17*x*y + 0.27599821010522485*x*y^2 - 4.6874560662018034e-17*x^2*y^3)^2 + (-0.4417735074073058*y - 5.425371220674707e-17*x - 0.10009637589813167*x*y + 1.353482119225013e-17*x*y^2 - 0.18335527863613874*x^2*y^3)^2"
     },
     "metadata": {},
     "execution_count": 7
    }
   ],
   "cell_type": "code",
   "source": [
    "sos_decomposition(con_ref)"
   ],
   "metadata": {},
   "execution_count": 7
  },
  {
   "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.10.3"
  },
  "kernelspec": {
   "name": "julia-1.10",
   "display_name": "Julia 1.10.3",
   "language": "julia"
  }
 },
 "nbformat": 4
}