{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "# Dihedral symmetry of the Robinson form"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "**Adapted from**: Example 5.4 of [GP04]\n",
    "\n",
    "[GP04] Gatermann, Karin and Parrilo, Pablo A.\n",
    "*Symmetry groups, semidefinite programs, and sums of squares*.\n",
    "Journal of Pure and Applied Algebra 192.1-3 (2004): 95-128."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "We start by defining the Dihedral group of order 8.\n",
    "This group is isomorphic to the following permutation group:"
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "Permutation group on 2 generators generated by\n (1,3)\n (1,2,3,4)"
     },
     "metadata": {},
     "execution_count": 1
    }
   ],
   "cell_type": "code",
   "source": [
    "using PermutationGroups\n",
    "d = perm\"(1, 2, 3, 4)\"\n",
    "c = perm\"(1, 3)\"\n",
    "G = PermGroup([c, d])"
   ],
   "metadata": {},
   "execution_count": 1
  },
  {
   "cell_type": "markdown",
   "source": [
    "We could rely on this isomorphism to define this group.\n",
    "However, in order to illustrate how to do symmetry reduction with a custom group,\n",
    "we show in this example what should be implemented to define a new group."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "import GroupsCore\n",
    "\n",
    "struct DihedralGroup <: GroupsCore.Group\n",
    "    n::Int\n",
    "end\n",
    "\n",
    "struct DihedralElement <: GroupsCore.GroupElement\n",
    "    n::Int\n",
    "    reflection::Bool\n",
    "    id::Int\n",
    "end"
   ],
   "metadata": {},
   "execution_count": 2
  },
  {
   "cell_type": "markdown",
   "source": [
    "Implementing GroupsCore API:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "Base.one(G::DihedralGroup) = DihedralElement(G.n, false, 0)\n",
    "\n",
    "Base.eltype(::Type{DihedralGroup}) = DihedralElement\n",
    "Base.IteratorSize(::Type{DihedralGroup}) = Base.HasLength()\n",
    "\n",
    "function Base.iterate(G::DihedralGroup, prev::DihedralElement=DihedralElement(G.n, false, -1))\n",
    "    if prev.id + 1 >= G.n\n",
    "        if prev.reflection\n",
    "            return nothing\n",
    "        else\n",
    "            next = DihedralElement(G.n, true, 0)\n",
    "        end\n",
    "    else\n",
    "        next = DihedralElement(G.n, prev.reflection, prev.id + 1)\n",
    "    end\n",
    "    return next, next\n",
    "end\n",
    "\n",
    "GroupsCore.order(::Type{T}, G::DihedralGroup) where {T} = convert(T, 2G.n)\n",
    "GroupsCore.gens(G::DihedralGroup) = [DihedralElement(G.n, false, 1), DihedralElement(G.n, true, 0)]"
   ],
   "metadata": {},
   "execution_count": 3
  },
  {
   "cell_type": "markdown",
   "source": [
    "Base.rand not needed for our purposes here"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "Base.parent(g::DihedralElement) = DihedralGroup(g.n)\n",
    "function Base.:(==)(g::DihedralElement, h::DihedralElement)\n",
    "    return g.n == h.n && g.reflection == h.reflection && g.id == h.id\n",
    "end\n",
    "\n",
    "function Base.inv(el::DihedralElement)\n",
    "    if el.reflection || iszero(el.id)\n",
    "        return el\n",
    "    else\n",
    "        return DihedralElement(el.n, false, el.n - el.id)\n",
    "    end\n",
    "end\n",
    "function Base.:*(a::DihedralElement, b::DihedralElement)\n",
    "    a.n == b.n || error(\"Cannot multiply elements from different Dihedral groups\")\n",
    "    id = mod(a.reflection ? a.id - b.id : a.id + b.id, a.n)\n",
    "    return DihedralElement(a.n, a.reflection != b.reflection, id)\n",
    "end\n",
    "\n",
    "Base.copy(a::DihedralElement) = DihedralElement(a.n, a.reflection, a.id)"
   ],
   "metadata": {},
   "execution_count": 4
  },
  {
   "cell_type": "markdown",
   "source": [
    "optional functions:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "function GroupsCore.order(el::DihedralElement)\n",
    "    if el.reflection\n",
    "        return 2\n",
    "    else\n",
    "        if iszero(el.id)\n",
    "            return 1\n",
    "        else\n",
    "            return div(el.n, gcd(el.n, el.id))\n",
    "        end\n",
    "    end\n",
    "end"
   ],
   "metadata": {},
   "execution_count": 5
  },
  {
   "cell_type": "markdown",
   "source": [
    "The Robinson form is invariant under the following action of the Dihedral group on monomials:\n",
    "The action of each element of the groups is to map the variables `x, y` to:\n",
    "\n",
    "| id | rotation | reflection |\n",
    "|----|----------|------------|\n",
    "| 0  | x, y     | y, x       |\n",
    "| 1  | -y, x    | -x, y      |\n",
    "| 2  | -x, -y   | -y, -x     |\n",
    "| 3  | y, -x    | x, -y      |"
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "x⁶ - x⁴y² - x²y⁴ + y⁶ - x⁴ + 3x²y² - y⁴ - x² - y² + 1",
      "text/latex": "$$ x^{6} - x^{4}y^{2} - x^{2}y^{4} + y^{6} - x^{4} + 3x^{2}y^{2} - y^{4} - x^{2} - y^{2} + 1 $$"
     },
     "metadata": {},
     "execution_count": 6
    }
   ],
   "cell_type": "code",
   "source": [
    "using SumOfSquares\n",
    "using DynamicPolynomials\n",
    "@polyvar x y\n",
    "struct DihedralAction <: Symmetry.OnMonomials end\n",
    "import SymbolicWedderburn\n",
    "SymbolicWedderburn.coeff_type(::DihedralAction) = Float64\n",
    "function SymbolicWedderburn.action(::DihedralAction, el::DihedralElement, mono::AbstractMonomial)\n",
    "    if iseven(el.reflection + el.id)\n",
    "        var_x, var_y = x, y\n",
    "    else\n",
    "        var_x, var_y = y, x\n",
    "    end\n",
    "    sign_x = 1 <= el.id <= 2 ? -1 : 1\n",
    "    sign_y = 2 <= el.id ? -1 : 1\n",
    "    return mono([x, y] => [sign_x * var_x, sign_y * var_y])\n",
    "end\n",
    "\n",
    "poly = x^6 + y^6 - x^4 * y^2 - y^4 * x^2 - x^4 - y^4 - x^2 - y^2 + 3x^2 * y^2 + 1"
   ],
   "metadata": {},
   "execution_count": 6
  },
  {
   "cell_type": "markdown",
   "source": [
    "We can verify that `poly` is indeed invariant under the action of each element of the group as follows."
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SymbolicWedderburn.action(DihedralAction(), g, poly) = x⁶ - x⁴y² - x²y⁴ + y⁶ - x⁴ + 3x²y² - y⁴ - x² - y² + 1\n",
      "SymbolicWedderburn.action(DihedralAction(), g, poly) = x⁶ - x⁴y² - x²y⁴ + y⁶ - x⁴ + 3x²y² - y⁴ - x² - y² + 1\n",
      "SymbolicWedderburn.action(DihedralAction(), g, poly) = x⁶ - x⁴y² - x²y⁴ + y⁶ - x⁴ + 3x²y² - y⁴ - x² - y² + 1\n",
      "SymbolicWedderburn.action(DihedralAction(), g, poly) = x⁶ - x⁴y² - x²y⁴ + y⁶ - x⁴ + 3x²y² - y⁴ - x² - y² + 1\n",
      "SymbolicWedderburn.action(DihedralAction(), g, poly) = x⁶ - x⁴y² - x²y⁴ + y⁶ - x⁴ + 3x²y² - y⁴ - x² - y² + 1\n",
      "SymbolicWedderburn.action(DihedralAction(), g, poly) = x⁶ - x⁴y² - x²y⁴ + y⁶ - x⁴ + 3x²y² - y⁴ - x² - y² + 1\n",
      "SymbolicWedderburn.action(DihedralAction(), g, poly) = x⁶ - x⁴y² - x²y⁴ + y⁶ - x⁴ + 3x²y² - y⁴ - x² - y² + 1\n",
      "SymbolicWedderburn.action(DihedralAction(), g, poly) = x⁶ - x⁴y² - x²y⁴ + y⁶ - x⁴ + 3x²y² - y⁴ - x² - y² + 1\n"
     ]
    }
   ],
   "cell_type": "code",
   "source": [
    "G = DihedralGroup(4)\n",
    "for g in G\n",
    "    @show SymbolicWedderburn.action(DihedralAction(), g, poly)\n",
    "end"
   ],
   "metadata": {},
   "execution_count": 7
  },
  {
   "cell_type": "markdown",
   "source": [
    "We can exploit this symmetry for reducing the problem using the `SymmetricIdeal` certificate as follows:"
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iter: 11 Ap: 9.58e-01 Pobj: -1.6713214e-10 Ad: 9.58e-01 Dobj: -4.5283865e-09 \n",
      "Success: SDP solved\n",
      "Primal objective value: -1.6713214e-10 \n",
      "Dual objective value: -4.5283865e-09 \n",
      "Relative primal infeasibility: 1.16e-14 \n",
      "Relative dual infeasibility: 1.17e-09 \n",
      "Real Relative Gap: -4.36e-09 \n",
      "XZ Relative Gap: 1.21e-09 \n",
      "DIMACS error measures: 2.00e-14 0.00e+00 2.80e-09 0.00e+00 -4.36e-09 1.21e-09\n",
      "CSDP 6.2.0\n",
      "Iter:  0 Ap: 0.00e+00 Pobj:  0.0000000e+00 Ad: 0.00e+00 Dobj:  0.0000000e+00 \n",
      "Iter:  1 Ap: 7.31e-01 Pobj: -5.9540693e+00 Ad: 8.42e-01 Dobj:  3.0671681e-01 \n",
      "Iter:  2 Ap: 8.58e-01 Pobj: -1.6732387e+01 Ad: 8.59e-01 Dobj:  7.8448086e+00 \n",
      "Iter:  3 Ap: 8.56e-01 Pobj: -6.5657892e+00 Ad: 7.36e-01 Dobj:  4.4153053e+00 \n",
      "Iter:  4 Ap: 7.93e-01 Pobj: -3.0377721e+00 Ad: 8.60e-01 Dobj:  8.3888624e-01 \n",
      "Iter:  5 Ap: 8.47e-01 Pobj: -1.8654993e+00 Ad: 8.10e-01 Dobj: -2.9577322e-01 \n",
      "Iter:  6 Ap: 7.39e-01 Pobj: -1.2273103e+00 Ad: 8.05e-01 Dobj: -6.7817689e-01 \n",
      "Iter:  7 Ap: 8.93e-01 Pobj: -1.0495527e+00 Ad: 6.76e-01 Dobj: -8.1824309e-01 \n",
      "Iter:  8 Ap: 6.81e-01 Pobj: -9.8438100e-01 Ad: 7.66e-01 Dobj: -8.8707837e-01 \n",
      "Iter:  9 Ap: 8.56e-01 Pobj: -9.5522629e-01 Ad: 7.27e-01 Dobj: -9.1564040e-01 \n",
      "Iter: 10 Ap: 1.00e+00 Pobj: -9.2940978e-01 Ad: 8.95e-01 Dobj: -9.2673695e-01 \n",
      "Iter: 11 Ap: 8.59e-01 Pobj: -9.2935253e-01 Ad: 9.06e-01 Dobj: -9.3052848e-01 \n",
      "Iter: 12 Ap: 1.00e+00 Pobj: -9.2917505e-01 Ad: 1.00e+00 Dobj: -9.3127516e-01 \n",
      "Iter: 13 Ap: 6.96e-01 Pobj: -9.2976072e-01 Ad: 5.87e-01 Dobj: -9.3164776e-01 \n",
      "Iter: 14 Ap: 1.00e+00 Pobj: -9.3001270e-01 Ad: 3.85e-01 Dobj: -9.3189163e-01 \n",
      "Iter: 15 Ap: 3.46e-01 Pobj: -9.3071764e-01 Ad: 1.00e+00 Dobj: -9.3190533e-01 \n",
      "Iter: 16 Ap: 1.00e+00 Pobj: -9.3058083e-01 Ad: 9.82e-01 Dobj: -9.3222213e-01 \n",
      "Iter: 17 Ap: 3.75e-01 Pobj: -9.3121390e-01 Ad: 1.00e+00 Dobj: -9.3246428e-01 \n",
      "Iter: 18 Ap: 4.04e-01 Pobj: -9.3131190e-01 Ad: 5.84e-01 Dobj: -9.3262908e-01 \n",
      "Iter: 19 Ap: 4.39e-01 Pobj: -9.3163773e-01 Ad: 2.42e-01 Dobj: -9.3271332e-01 \n",
      "Iter: 20 Ap: 3.16e-01 Pobj: -9.3173799e-01 Ad: 4.08e-01 Dobj: -9.3277433e-01 \n",
      "Iter: 21 Ap: 6.55e-01 Pobj: -9.3179189e-01 Ad: 6.39e-01 Dobj: -9.3280691e-01 \n",
      "Iter: 22 Ap: 1.00e+00 Pobj: -9.3180218e-01 Ad: 1.00e+00 Dobj: -9.3281932e-01 \n",
      "Iter: 23 Ap: 7.71e-01 Pobj: -9.3182701e-01 Ad: 4.94e-01 Dobj: -9.3282851e-01 \n",
      "Iter: 24 Ap: 1.00e+00 Pobj: -9.3182027e-01 Ad: 9.25e-01 Dobj: -9.3282972e-01 \n",
      "value(t) = -0.9318227936087737\n",
      "DynamicPolynomials.Polynomial{true, Float64}[xy², x, x³]\n",
      "DynamicPolynomials.Polynomial{true, Float64}[x²y, y, y³]\n",
      "DynamicPolynomials.Polynomial{true, Float64}[1.0, -0.7071067811865472x² - 0.7071067811865475y²]\n",
      "DynamicPolynomials.Polynomial{true, Float64}[xy]\n",
      "DynamicPolynomials.Polynomial{true, Float64}[-0.7071067811865472x² + 0.7071067811865475y²]\n"
     ]
    }
   ],
   "cell_type": "code",
   "source": [
    "import CSDP\n",
    "function solve(G)\n",
    "    solver = CSDP.Optimizer\n",
    "    model = Model(solver)\n",
    "    @variable(model, t)\n",
    "    @objective(model, Max, t)\n",
    "    pattern = Symmetry.Pattern(G, DihedralAction())\n",
    "    con_ref = @constraint(model, poly - t in SOSCone(), symmetry = pattern)\n",
    "    optimize!(model)\n",
    "    @show value(t)\n",
    "\n",
    "\n",
    "    for g in gram_matrix(con_ref).sub_gram_matrices\n",
    "        println(g.basis.polynomials)\n",
    "    end\n",
    "end\n",
    "solve(G)"
   ],
   "metadata": {},
   "execution_count": 8
  },
  {
   "cell_type": "markdown",
   "source": [
    "We notice that we indeed find `-3825/4096` and that symmetry was exploited."
   ],
   "metadata": {}
  },
  {
   "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
}