{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "# Cyclic symmetry for Sums of Hermitian Squares"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "We start by defining the Cyclic group."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using GroupsCore\n",
    "import PermutationGroups\n",
    "\n",
    "struct CyclicElem <: GroupElement\n",
    "    n::Int\n",
    "    id::Int\n",
    "end\n",
    "Base.:(==)(a::CyclicElem, b::CyclicElem) = a.n == b.n && a.id == b.id\n",
    "Base.inv(el::CyclicElem) = CyclicElem(el.n, (el.n - el.id) % el.n)\n",
    "\n",
    "function Base.:*(a::CyclicElem, b::CyclicElem)\n",
    "    return CyclicElem(a.n, (a.id + b.id) % a.n)\n",
    "end\n",
    "Base.:^(el::CyclicElem, k::Integer) = CyclicElem(el.n, (el.id * k) % el.n)\n",
    "\n",
    "Base.conj(a::CyclicElem, b::CyclicElem) = inv(b) * a * b\n",
    "Base.:^(a::CyclicElem, b::CyclicElem) = conj(a, b)\n",
    "\n",
    "function PermutationGroups.order(el::CyclicElem)\n",
    "    return div(el.n, gcd(el.n, el.id))\n",
    "end\n",
    "\n",
    "struct CyclicGroup <: Group\n",
    "    n::Int\n",
    "end\n",
    "Base.eltype(::CyclicGroup) = CyclicElem\n",
    "Base.one(c::Union{CyclicGroup, CyclicElem}) = CyclicElem(c.n, 0)\n",
    "PermutationGroups.gens(c::CyclicGroup) = [CyclicElem(c.n, 1)]\n",
    "PermutationGroups.order(::Type{T}, c::CyclicGroup) where {T} = convert(T, c.n)\n",
    "function Base.iterate(c::CyclicGroup, prev::CyclicElem=CyclicElem(c.n, -1))\n",
    "    id = prev.id + 1\n",
    "    if id >= c.n\n",
    "        return nothing\n",
    "    else\n",
    "        next = CyclicElem(c.n, id)\n",
    "        return next, next\n",
    "    end\n",
    "end"
   ],
   "metadata": {},
   "execution_count": 1
  },
  {
   "cell_type": "markdown",
   "source": [
    "Now we define that the cyclic group acts on monomial by permuting variables\n",
    "cyclically. So for instance, `CyclicElem(3, 1)` would transform\n",
    "`x_1^3*x_2*x_3^4` into `x_1^4*x_2^3*x_3`."
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "x₁⁴x₂³x₃",
      "text/latex": "$$ x_{1}^{4}x_{2}^{3}x_{3} $$"
     },
     "metadata": {},
     "execution_count": 2
    }
   ],
   "cell_type": "code",
   "source": [
    "import MultivariatePolynomials as MP\n",
    "\n",
    "using SumOfSquares\n",
    "\n",
    "struct Action{V<:MP.AbstractVariable} <: Symmetry.OnMonomials\n",
    "    variables::Vector{V}\n",
    "end\n",
    "Symmetry.SymbolicWedderburn.coeff_type(::Action) = Float64\n",
    "function Symmetry.SymbolicWedderburn.action(a::Action, el::CyclicElem, mono::MP.AbstractMonomial)\n",
    "    return prod(MP.powers(mono), init=MP.constant_monomial(mono)) do (var, exp)\n",
    "        index = findfirst(isequal(var), a.variables)\n",
    "        new_index = mod1(index + el.id, el.n)\n",
    "        return a.variables[new_index]^exp\n",
    "    end\n",
    "end\n",
    "\n",
    "using DynamicPolynomials\n",
    "@polyvar x[1:3]\n",
    "action = Action(x)\n",
    "g = CyclicElem(3, 1)\n",
    "Symmetry.SymbolicWedderburn.action(action, g, x[1]^3 * x[2] * x[3]^4)"
   ],
   "metadata": {},
   "execution_count": 2
  },
  {
   "cell_type": "markdown",
   "source": [
    "The following polynomial `poly` is invariant under the action of the group `G`."
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "x₃² + x₂x₃ + x₂² + x₁x₃ + x₁x₂ + x₁²",
      "text/latex": "$$ x_{3}^{2} + x_{2}x_{3} + x_{2}^{2} + x_{1}x_{3} + x_{1}x_{2} + x_{1}^{2} $$"
     },
     "metadata": {},
     "execution_count": 3
    }
   ],
   "cell_type": "code",
   "source": [
    "N = 3\n",
    "G = CyclicGroup(N)\n",
    "poly = sum(x[i] * x[mod1(i + 1, N)] for i in 1:N) + sum(x.^2)\n",
    "Symmetry.SymbolicWedderburn.action(action, g, poly)"
   ],
   "metadata": {},
   "execution_count": 3
  },
  {
   "cell_type": "markdown",
   "source": [
    "Let's now find the minimum of `p` by exploiting this symmetry."
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iter: 18 Ap: 9.59e-01 Pobj: -3.5512000e-03 Ad: 9.59e-01 Dobj: -3.5512523e-03 \n",
      "Success: SDP solved\n",
      "Primal objective value: -3.5512000e-03 \n",
      "Dual objective value: -3.5512523e-03 \n",
      "Relative primal infeasibility: 1.17e-13 \n",
      "Relative dual infeasibility: 5.36e-10 \n",
      "Real Relative Gap: -5.20e-08 \n",
      "XZ Relative Gap: 1.22e-09 \n",
      "DIMACS error measures: 1.26e-13 0.00e+00 1.44e-09 0.00e+00 -5.20e-08 1.22e-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: 1.00e+00 Pobj: -6.8633327e+00 Ad: 9.00e-01 Dobj:  1.9835999e+01 \n",
      "Iter:  2 Ap: 1.00e+00 Pobj: -2.8142041e+00 Ad: 9.00e-01 Dobj:  3.4715380e+00 \n",
      "Iter:  3 Ap: 1.00e+00 Pobj: -1.5335395e-02 Ad: 9.00e-01 Dobj:  4.1414777e-01 \n",
      "Iter:  4 Ap: 1.00e+00 Pobj: -2.4302191e-05 Ad: 1.00e+00 Dobj:  4.4865070e-05 \n",
      "Iter:  5 Ap: 1.00e+00 Pobj: -7.8953240e-07 Ad: 1.00e+00 Dobj: -1.2669850e-06 \n",
      "Iter:  6 Ap: 1.00e+00 Pobj: -5.8902884e-08 Ad: 1.00e+00 Dobj: -5.8810369e-07 \n",
      "Iter:  7 Ap: 1.00e+00 Pobj: -8.1684162e-09 Ad: 1.00e+00 Dobj:  1.9520853e-09 \n"
     ]
    },
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "* Solver : CSDP\n\n* Status\n  Result count       : 1\n  Termination status : OPTIMAL\n  Message from the solver:\n  \"Problem solved to optimality.\"\n\n* Candidate solution (result #1)\n  Primal status      : FEASIBLE_POINT\n  Dual status        : FEASIBLE_POINT\n  Objective value    : 2.22045e-16\n  Dual objective value : -2.20229e-10\n\n* Work counters\n  Solve time (sec)   : 1.07503e-03\n"
     },
     "metadata": {},
     "execution_count": 4
    }
   ],
   "cell_type": "code",
   "source": [
    "import CSDP\n",
    "solver = CSDP.Optimizer\n",
    "model = Model(solver)\n",
    "@variable(model, t)\n",
    "@objective(model, Max, t)\n",
    "pattern = Symmetry.Pattern(G, action)\n",
    "con_ref = @constraint(model, poly - t in SOSCone(), symmetry = pattern)\n",
    "optimize!(model)\n",
    "solution_summary(model)"
   ],
   "metadata": {},
   "execution_count": 4
  },
  {
   "cell_type": "markdown",
   "source": [
    "Let's look at the symmetry adapted basis used."
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "DynamicPolynomials.Polynomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}, Float64}[1.0, -0.5773502691896257x₃ - 0.5773502691896257x₂ - 0.5773502691896257x₁]\n",
      "DynamicPolynomials.Polynomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}, Float64}[-0.816496580927726x₃ + 0.4082482904638631x₂ + 0.4082482904638631x₁]\n",
      "DynamicPolynomials.Polynomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}, Float64}[-0.7071067811865477x₂ + 0.7071067811865477x₁]\n"
     ]
    },
    {
     "output_type": "display_data",
     "data": {
      "text/plain": "2×2 SymMatrix{Float64}:\n 1.07382e-17  0.0\n 0.0          2.0"
     },
     "metadata": {}
    },
    {
     "output_type": "display_data",
     "data": {
      "text/plain": "1×1 SymMatrix{Float64}:\n 0.4999999999999981"
     },
     "metadata": {}
    },
    {
     "output_type": "display_data",
     "data": {
      "text/plain": "1×1 SymMatrix{Float64}:\n 0.4999999999999981"
     },
     "metadata": {}
    }
   ],
   "cell_type": "code",
   "source": [
    "for gram in gram_matrix(con_ref).blocks\n",
    "    println(gram.basis.polynomials)\n",
    "    display(gram.Q)\n",
    "end"
   ],
   "metadata": {},
   "execution_count": 5
  },
  {
   "cell_type": "markdown",
   "source": [
    "Let's look into more details at the last two elements of the basis."
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "2-element Vector{DynamicPolynomials.Polynomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}, Float64}}:\n -0.8164965809277261x₃ + 0.4082482904638631x₂ + 0.4082482904638631x₁\n -0.7071067811865475x₂ + 0.7071067811865475x₁"
     },
     "metadata": {},
     "execution_count": 6
    }
   ],
   "cell_type": "code",
   "source": [
    "basis = [(x[1] + x[2] - 2x[3])/√6, (x[1] - x[2])/√2]"
   ],
   "metadata": {},
   "execution_count": 6
  },
  {
   "cell_type": "markdown",
   "source": [
    "This actually constitutes the basis for an invariant subspace corresponding\n",
    "to a group character of degree 2 and multiplicity 1.\n",
    "This means that it decomposes the semidefinite matrix into 2 blocks of size\n",
    "1-by-1 that are equal. Indeed, we see above that `gram.Q` is identically equal for both.\n",
    "As the group is generated by one element `g`, we can just verify it by verifying\n",
    "its invariance under `g`.\n",
    "The image of each element under the basis is:"
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "2-element Vector{DynamicPolynomials.Polynomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}, Float64}}:\n 0.4082482904638631x₃ + 0.4082482904638631x₂ - 0.8164965809277261x₁\n -0.7071067811865475x₃ + 0.7071067811865475x₂"
     },
     "metadata": {},
     "execution_count": 7
    }
   ],
   "cell_type": "code",
   "source": [
    "image = [Symmetry.SymbolicWedderburn.action(action, g, p) for p in basis]"
   ],
   "metadata": {},
   "execution_count": 7
  },
  {
   "cell_type": "markdown",
   "source": [
    "We can see that they are both still in the same 2-dimensional subspace."
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "2-element Vector{DynamicPolynomials.Polynomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}, Float64}}:\n 0.4082482904638631x₃ + 0.4082482904638629x₂ - 0.816496580927726x₁\n -0.7071067811865476x₃ + 0.7071067811865475x₂ + 5.551115123125783e-17x₁"
     },
     "metadata": {},
     "execution_count": 8
    }
   ],
   "cell_type": "code",
   "source": [
    "a = -1/2\n",
    "b = √3/2\n",
    "[a -b; b a] * basis"
   ],
   "metadata": {},
   "execution_count": 8
  },
  {
   "cell_type": "markdown",
   "source": [
    "In fact, these last two basis comes from the real decomposition of a complex one."
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iter:  8 Ap: 1.00e+00 Pobj:  2.2204460e-16 Ad: 9.00e-01 Dobj: -2.2022945e-10 \n",
      "Success: SDP solved\n",
      "Primal objective value: 2.2204460e-16 \n",
      "Dual objective value: -2.2022945e-10 \n",
      "Relative primal infeasibility: 8.86e-16 \n",
      "Relative dual infeasibility: 3.82e-10 \n",
      "Real Relative Gap: -2.20e-10 \n",
      "XZ Relative Gap: 2.29e-09 \n",
      "DIMACS error measures: 1.53e-15 0.00e+00 7.89e-10 0.00e+00 -2.20e-10 2.29e-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: 1.00e+00 Pobj: -8.9892401e+00 Ad: 9.07e-01 Dobj:  4.8627163e+01 \n",
      "Iter:  2 Ap: 1.00e+00 Pobj: -5.4165963e+00 Ad: 9.18e-01 Dobj:  1.0591774e+01 \n",
      "Iter:  3 Ap: 1.00e+00 Pobj: -6.3166907e-01 Ad: 8.98e-01 Dobj:  3.0288603e+00 \n",
      "Iter:  4 Ap: 1.00e+00 Pobj: -4.2564141e-02 Ad: 9.09e-01 Dobj:  3.9201832e-01 \n",
      "Iter:  5 Ap: 1.00e+00 Pobj: -4.3454962e-03 Ad: 9.28e-01 Dobj:  4.0461107e-02 \n",
      "Iter:  6 Ap: 1.00e+00 Pobj: -4.2815474e-04 Ad: 9.34e-01 Dobj:  3.8633036e-03 \n",
      "Iter:  7 Ap: 1.00e+00 Pobj: -4.3920196e-05 Ad: 9.88e-01 Dobj:  1.7392285e-04 \n",
      "Iter:  8 Ap: 1.00e+00 Pobj: -2.3698636e-06 Ad: 9.99e-01 Dobj:  4.5447245e-06 \n",
      "Iter:  9 Ap: 1.00e+00 Pobj: -9.4705116e-08 Ad: 1.00e+00 Dobj: -2.3966141e-06 \n",
      "Iter: 10 Ap: 1.00e+00 Pobj: -4.0174413e-09 Ad: 9.69e-01 Dobj: -9.9538646e-08 \n",
      "DynamicPolynomials.Polynomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}, ComplexF64}[(1.0 - 0.0im), (-0.5773502691896257 - 0.0im)x₃ + (-0.5773502691896257 - 0.0im)x₂ + (-0.5773502691896257 - 0.0im)x₁]\n",
      "DynamicPolynomials.Polynomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}, ComplexF64}[(-0.577350269189626 - 0.0im)x₃ + (0.2886751345948128 + 0.5000000000000001im)x₂ + (0.28867513459481314 - 0.49999999999999994im)x₁]\n",
      "DynamicPolynomials.Polynomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}, ComplexF64}[(-0.577350269189626 - 0.0im)x₃ + (0.28867513459481314 - 0.49999999999999994im)x₂ + (0.2886751345948128 + 0.5000000000000001im)x₁]\n"
     ]
    },
    {
     "output_type": "display_data",
     "data": {
      "text/plain": "2×2 VectorizedHermitianMatrix{ComplexF64, Float64, ComplexF64}:\n 1.67745e-10+0.0im  0.0+0.0im\n         0.0+0.0im  2.0+0.0im"
     },
     "metadata": {}
    },
    {
     "output_type": "display_data",
     "data": {
      "text/plain": "1×1 VectorizedHermitianMatrix{ComplexF64, Float64, ComplexF64}:\n 0.5000000408028554 + 0.0im"
     },
     "metadata": {}
    },
    {
     "output_type": "display_data",
     "data": {
      "text/plain": "1×1 VectorizedHermitianMatrix{ComplexF64, Float64, ComplexF64}:\n 0.4999999591971749 + 0.0im"
     },
     "metadata": {}
    }
   ],
   "cell_type": "code",
   "source": [
    "import CSDP\n",
    "solver = CSDP.Optimizer\n",
    "model = Model(solver)\n",
    "@variable(model, t)\n",
    "@objective(model, Max, t)\n",
    "pattern = Symmetry.Pattern(G, action)\n",
    "cone = SumOfSquares.NonnegPolyInnerCone{MOI.HermitianPositiveSemidefiniteConeTriangle}()\n",
    "con_ref = @constraint(model, poly - t in cone, symmetry = pattern)\n",
    "optimize!(model)\n",
    "solution_summary(model)\n",
    "\n",
    "for gram in gram_matrix(con_ref).blocks\n",
    "    println(gram.basis.polynomials)\n",
    "    display(gram.Q)\n",
    "end"
   ],
   "metadata": {},
   "execution_count": 9
  },
  {
   "cell_type": "markdown",
   "source": [
    "We can see that the real invariant subspace was in fact coming from two complex conjugate complex invariant subspaces:"
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "(0.4082482904638631 - 0.7071067811865475im)x₃ + (0.4082482904638631 + 0.7071067811865475im)x₂ + (-0.8164965809277261 + 0.0im)x₁",
      "text/latex": "$$ (0.4082482904638631 - 0.7071067811865475im)x_{3} + (0.4082482904638631 + 0.7071067811865475im)x_{2} + (-0.8164965809277261 + 0.0im)x_{1} $$"
     },
     "metadata": {},
     "execution_count": 10
    }
   ],
   "cell_type": "code",
   "source": [
    "complex_basis = basis[1] + im * basis[2]\n",
    "image = Symmetry.SymbolicWedderburn.action(action, g, complex_basis)"
   ],
   "metadata": {},
   "execution_count": 10
  },
  {
   "cell_type": "markdown",
   "source": [
    "And there is a direct correspondance between the representation of the real and complex versions:"
   ],
   "metadata": {}
  },
  {
   "outputs": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "(0.4082482904638631 - 0.7071067811865476im)x₃ + (0.4082482904638629 + 0.7071067811865475im)x₂ + (-0.816496580927726 + 5.551115123125783e-17im)x₁",
      "text/latex": "$$ (0.4082482904638631 - 0.7071067811865476im)x_{3} + (0.4082482904638629 + 0.7071067811865475im)x_{2} + (-0.816496580927726 + 5.551115123125783e-17im)x_{1} $$"
     },
     "metadata": {},
     "execution_count": 11
    }
   ],
   "cell_type": "code",
   "source": [
    "(a + b * im) * complex_basis"
   ],
   "metadata": {},
   "execution_count": 11
  },
  {
   "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
}