{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "# Hypercube"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "**Contributed by**: Benoît Legat"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Given a Sum-of-Squares constraint on an algebraic set:\n",
    "$$\n",
    "g_1(x) = , \\ldots, g_m(x) = 0 \\Rightarrow p(x) \\ge 0.\n",
    "$$\n",
    "We can either use the certificate:\n",
    "$$\n",
    "p(x) = s(x) + \\lambda_1(x) g_1(x) + \\cdots + \\lambda_m(x) g_m(x), s_0(x) \\text{ is SOS},\n",
    "$$\n",
    "or\n",
    "$$\n",
    "p(x) \\equiv s(x) \\pmod{\\la g_1(x), \\ldots, g_m(x) \\ra}, s_0(x) \\text{ is SOS}.\n",
    "$$\n",
    "the second one leads to a *simpler* SDP but needs to compute a *Gr\\\"obner* basis:\n",
    " * SemialgebraicSets implements Buchberger's algorithm.\n",
    " * The `@set` macro recognizes variable fix, e.g., `x = 1`\n",
    "    and provides shortcut.\n",
    " * If you know a \\alert{better} way to take modulo,\n",
    "   better create your \\alert{own} type of algebraic set!"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "We illustrate this in this example."
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "Algebraic Set defined by 3 equalities\n -1//1 + x[1]^2 = 0\n -1//1 + x[2]^2 = 0\n -1//1 + x[3]^2 = 0\n"
     },
     "metadata": {},
     "execution_count": 1
    }
   ],
   "cell_type": "code",
   "source": [
    "using DynamicPolynomials\n",
    "@polyvar x[1:3]\n",
    "p = sum(x)^2\n",
    "using SumOfSquares\n",
    "S = algebraicset([xi^2 - 1 for xi in x])"
   ],
   "metadata": {},
   "execution_count": 1
  },
  {
   "cell_type": "markdown",
   "source": [
    "We will now search for the minimum of `x` over `S` using Sum of Squares Programming.\n",
    "We first 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."
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "termination_status(model) = MathOptInterface.OPTIMAL\n",
      "objective_value(model) = -4.404406006575101e-10\n"
     ]
    },
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "-4.404406006575101e-10"
     },
     "metadata": {},
     "execution_count": 2
    }
   ],
   "cell_type": "code",
   "source": [
    "import CSDP\n",
    "solver = optimizer_with_attributes(CSDP.Optimizer, MOI.Silent() => true)\n",
    "\n",
    "function min_algebraic(S)\n",
    "    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)\n",
    "end\n",
    "\n",
    "min_algebraic(S)"
   ],
   "metadata": {},
   "execution_count": 2
  },
  {
   "cell_type": "markdown",
   "source": [
    "Note that the minimum is in fact `1`.\n",
    "Indeed, since each variables is odd (it is either `-1` or `1`)\n",
    "and there is an odd number of variables, their sum is odd.\n",
    "Therefore it cannot be zero!"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "We can see that the Gröbner basis of `S` was computed"
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "S.I.gröbner_basis = false\n"
     ]
    },
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "Buchberger(1.4901161193847656e-8, SemialgebraicSets.presort!, SemialgebraicSets.normal_selection)"
     },
     "metadata": {},
     "execution_count": 3
    }
   ],
   "cell_type": "code",
   "source": [
    "@show S.I.gröbner_basis\n",
    "S.I.algo"
   ],
   "metadata": {},
   "execution_count": 3
  },
  {
   "cell_type": "markdown",
   "source": [
    "The Gröbner basis is simple to compute in this case as the vector\n",
    "of `xi^2 - 1` is already a Gröbner basis.\n",
    "However, we still need to divide polynomials by the Gröbner basis\n",
    "which can be simplified in this case."
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "termination_status(model) = MathOptInterface.OPTIMAL\n",
      "objective_value(model) = -4.404406006575101e-10\n"
     ]
    },
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "-4.404406006575101e-10"
     },
     "metadata": {},
     "execution_count": 4
    }
   ],
   "cell_type": "code",
   "source": [
    "const MP = MultivariatePolynomials\n",
    "const SS = SemialgebraicSets\n",
    "struct HypercubeIdeal{V} <: SS.AbstractPolynomialIdeal\n",
    "    variables::Vector{V}\n",
    "end\n",
    "struct HypercubeSet{V} <: SS.AbstractAlgebraicSet\n",
    "    ideal::HypercubeIdeal{V}\n",
    "end\n",
    "MP.variables(set::HypercubeSet) = MP.variables(set.ideal)\n",
    "MP.variables(ideal::HypercubeIdeal) = ideal.variables\n",
    "Base.similar(set::HypercubeSet, ::Type) = set\n",
    "SS.ideal(set::HypercubeSet) = set.ideal\n",
    "function Base.rem(p, set::HypercubeIdeal)\n",
    "    return MP.polynomial(map(MP.terms(p)) do term\n",
    "        mono = MP.monomial(term)\n",
    "        new_mono = one(mono)\n",
    "        for (var, exp) in powers(mono)\n",
    "            if var in set.variables\n",
    "                exp = rem(exp, 2)\n",
    "            end\n",
    "            new_mono *= var^exp\n",
    "        end\n",
    "        MP.coefficient(term) * new_mono\n",
    "    end)\n",
    "end\n",
    "\n",
    "H = HypercubeSet(HypercubeIdeal(x))\n",
    "\n",
    "min_algebraic(H)"
   ],
   "metadata": {},
   "execution_count": 4
  },
  {
   "cell_type": "markdown",
   "source": [
    "Let's now try to find the correct lower bound:"
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "min_algebraic_rational (generic function with 1 method)"
     },
     "metadata": {},
     "execution_count": 5
    }
   ],
   "cell_type": "code",
   "source": [
    "function min_algebraic_rational(S, d)\n",
    "    model = SOSModel(solver)\n",
    "    @variable(model, q, SOSPoly(MP.monomials(x, 0:d)))\n",
    "    deno = q + 1\n",
    "    @constraint(model, c, deno * p >= deno, domain = S)\n",
    "    optimize!(model)\n",
    "    @show termination_status(model)\n",
    "end"
   ],
   "metadata": {},
   "execution_count": 5
  },
  {
   "cell_type": "markdown",
   "source": [
    "With `d = 0`, it's the same as previously"
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "termination_status(model) = MathOptInterface.INFEASIBLE\n"
     ]
    },
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "INFEASIBLE::TerminationStatusCode = 2"
     },
     "metadata": {},
     "execution_count": 6
    }
   ],
   "cell_type": "code",
   "source": [
    "min_algebraic_rational(H, 0)"
   ],
   "metadata": {},
   "execution_count": 6
  },
  {
   "cell_type": "markdown",
   "source": [
    "But with `d = 1`, we can find the correct lower bound"
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "termination_status(model) = MathOptInterface.OPTIMAL\n"
     ]
    },
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "OPTIMAL::TerminationStatusCode = 1"
     },
     "metadata": {},
     "execution_count": 7
    }
   ],
   "cell_type": "code",
   "source": [
    "min_algebraic_rational(H, 1)"
   ],
   "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.9.3"
  },
  "kernelspec": {
   "name": "julia-1.9",
   "display_name": "Julia 1.9.3",
   "language": "julia"
  }
 },
 "nbformat": 4
}