{
 "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",
    "import MultivariateBases as MB\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": [
      "-------------------------------------------------------------\n",
      "           Clarabel.jl v0.10.0  -  Clever Acronym              \n",
      "                   (c) Paul Goulart                          \n",
      "                University of Oxford, 2022                   \n",
      "-------------------------------------------------------------\n",
      "\n",
      "problem:\n",
      "  variables     = 5\n",
      "  constraints   = 14\n",
      "  nnz(P)        = 0\n",
      "  nnz(A)        = 28\n",
      "  cones (total) = 3\n",
      "    : Zero        = 1,  numel = 10\n",
      "    : Nonnegative = 1,  numel = 1\n",
      "    : SecondOrder = 1,  numel = 3\n",
      "\n",
      "settings:\n",
      "  linear algebra: direct / qdldl, precision: Float64\n",
      "  max iter = 200, time limit = Inf,  max step = 0.990\n",
      "  tol_feas = 1.0e-08, tol_gap_abs = 1.0e-08, tol_gap_rel = 1.0e-08,\n",
      "  static reg : on, ϵ1 = 1.0e-08, ϵ2 = 4.9e-32\n",
      "  dynamic reg: on, ϵ = 1.0e-13, δ = 2.0e-07\n",
      "  iter refine: on, reltol = 1.0e-13, abstol = 1.0e-12, \n",
      "               max iter = 10, stop ratio = 5.0\n",
      "  equilibrate: on, min_scale = 1.0e-04, max_scale = 1.0e+04\n",
      "               max iter = 10\n",
      "\n",
      "iter    pcost        dcost       gap       pres      dres      k/t        μ       step      \n",
      "---------------------------------------------------------------------------------------------\n",
      "  0   3.0203e-16   2.6765e-16  3.44e-17  2.27e-01  2.62e-01  1.00e+00  1.87e+00   ------   \n",
      "  1   8.8771e-03  -3.1810e-02  4.07e-02  1.63e-02  1.86e-02  1.89e-02  1.59e-01  9.34e-01  \n",
      "  2   7.3065e-05  -2.2596e-04  2.99e-04  1.67e-04  1.91e-04  3.05e-04  1.76e-03  9.89e-01  \n",
      "  3   7.3003e-07  -2.2578e-06  2.99e-06  1.67e-06  1.91e-06  3.05e-06  1.76e-05  9.90e-01  \n",
      "  4   7.3003e-09  -2.2577e-08  2.99e-08  1.67e-08  1.91e-08  3.05e-08  1.76e-07  9.90e-01  \n",
      "  5   7.3005e-11  -2.2577e-10  2.99e-10  1.67e-10  1.91e-10  3.05e-10  1.76e-09  9.90e-01  \n",
      "---------------------------------------------------------------------------------------------\n",
      "Terminated with status = solved\n",
      "solve time =  538ms\n"
     ]
    },
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "solution_summary(; result = 1, verbose = false)\n├ solver_name          : Clarabel\n├ Termination\n│ ├ termination_status : OPTIMAL\n│ ├ result_count       : 1\n│ └ raw_status         : SOLVED\n├ Solution (result = 1)\n│ ├ primal_status        : FEASIBLE_POINT\n│ ├ dual_status          : FEASIBLE_POINT\n│ ├ objective_value      : -7.30048e-11\n│ └ dual_objective_value : 2.25773e-10\n└ Work counters\n  ├ solve_time (sec)   : 5.37622e-01\n  └ barrier_iterations : 5"
     },
     "metadata": {},
     "execution_count": 4
    }
   ],
   "cell_type": "code",
   "source": [
    "import Clarabel\n",
    "solver = Clarabel.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": [
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "BlockDiagonalGramMatrix with 2 blocks:\n[1] Block with row/column basis:\n     Simple basis:\n  FixedBasis([1.0·1 + 0.0·x₃ + 0.0·x₂ + 0.0·x₁, 0.0·1 - 0.5773502691896257·x₃ - 0.5773502691896257·x₂ - 0.5773502691896257·x₁])\n    And entries in a 2×2 SymMatrix{Float64}:\n     7.300482440797396e-11  0.0\n     0.0                    2.0\n[2] Block with row/column basis:\n     Semisimple basis with 2 simple sub-bases:\n  FixedBasis([0.0·1 - 0.816496580927726·x₃ + 0.4082482904638631·x₂ + 0.4082482904638631·x₁])\n  FixedBasis([0.0·1 + 0.0·x₃ - 0.7071067811865477·x₂ + 0.7071067811865477·x₁])\n    And entries in a 1×1 SymMatrix{Float64}:\n     0.4999999999999991"
     },
     "metadata": {},
     "execution_count": 5
    }
   ],
   "cell_type": "code",
   "source": [
    "gram_matrix(con_ref)"
   ],
   "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": [
      "-------------------------------------------------------------\n",
      "           Clarabel.jl v0.10.0  -  Clever Acronym              \n",
      "                   (c) Paul Goulart                          \n",
      "                University of Oxford, 2022                   \n",
      "-------------------------------------------------------------\n",
      "\n",
      "chordal decomposition:\n",
      "  compact format = on, dual completion = on\n",
      "  merge method = clique_graph\n",
      "  PSD cones initial             = 3\n",
      "  PSD cones decomposable        = 1\n",
      "  PSD cones after decomposition = 4\n",
      "  PSD cones after merges        = 4\n",
      "\n",
      "problem:\n",
      "  variables     = 10\n",
      "  constraints   = 28\n",
      "  nnz(P)        = 0\n",
      "  nnz(A)        = 41\n",
      "  cones (total) = 5\n",
      "    : Zero        = 1,  numel = 10\n",
      "    : PSDTriangle = 4,  numel = (6,6,3,3)\n",
      "\n",
      "settings:\n",
      "  linear algebra: direct / qdldl, precision: Float64\n",
      "  max iter = 200, time limit = Inf,  max step = 0.990\n",
      "  tol_feas = 1.0e-08, tol_gap_abs = 1.0e-08, tol_gap_rel = 1.0e-08,\n",
      "  static reg : on, ϵ1 = 1.0e-08, ϵ2 = 4.9e-32\n",
      "  dynamic reg: on, ϵ = 1.0e-13, δ = 2.0e-07\n",
      "  iter refine: on, reltol = 1.0e-13, abstol = 1.0e-12, \n",
      "               max iter = 10, stop ratio = 5.0\n",
      "  equilibrate: on, min_scale = 1.0e-04, max_scale = 1.0e+04\n",
      "               max iter = 10\n",
      "\n",
      "iter    pcost        dcost       gap       pres      dres      k/t        μ       step      \n",
      "---------------------------------------------------------------------------------------------\n",
      "  0   0.0000e+00  -0.0000e+00  0.00e+00  3.72e-01  5.55e-01  1.00e+00  1.72e+00   ------   \n",
      "  1   1.7404e-01   1.6823e-01  5.82e-03  5.08e-02  8.11e-02  9.32e-02  3.09e-01  9.04e-01  \n",
      "  2   8.1131e-04   2.5968e-03  1.79e-03  6.25e-04  1.00e-03  2.96e-03  4.26e-03  9.89e-01  \n",
      "  3   8.0996e-06   2.5932e-05  1.78e-05  6.24e-06  1.00e-05  2.95e-05  4.26e-05  9.90e-01  \n",
      "  4   8.0995e-08   2.5932e-07  1.78e-07  6.24e-08  1.00e-07  2.95e-07  4.26e-07  9.90e-01  \n",
      "  5   8.0995e-10   2.5932e-09  1.78e-09  6.24e-10  1.00e-09  2.95e-09  4.26e-09  9.90e-01  \n",
      "---------------------------------------------------------------------------------------------\n",
      "Terminated with status = solved\n",
      "solve time =  1.65s\n"
     ]
    },
    {
     "output_type": "execute_result",
     "data": {
      "text/plain": "BlockDiagonalGramMatrix with 3 blocks:\n[1] Block with row/column basis:\n     Simple basis:\n  FixedBasis([(1.0 - 0.0im)·1 + (0.0 + 0.0im)·x₃ + (0.0 + 0.0im)·x₂ + (0.0 + 0.0im)·x₁, (0.0 + 0.0im)·1 + (-0.5773502691896257 - 0.0im)·x₃ + (-0.5773502691896257 - 0.0im)·x₂ + (-0.5773502691896257 - 0.0im)·x₁])\n    And entries in a 2×2 VectorizedHermitianMatrix{ComplexF64, Float64, ComplexF64}:\n     8.099519427648694e-10 + 0.0im                 0.0 + 0.0im\n                       0.0 + 0.0im  2.0000000000000013 + 0.0im\n[2] Block with row/column basis:\n     Simple basis:\n  FixedBasis([(0.0 + 0.0im)·1 + (-0.577350269189626 - 0.0im)·x₃ + (0.28867513459481314 - 0.49999999999999994im)·x₂ + (0.2886751345948128 + 0.5000000000000001im)·x₁])\n    And entries in a 1×1 VectorizedHermitianMatrix{ComplexF64, Float64, ComplexF64}:\n     0.4999999999999958 + 0.0im\n[3] Block with row/column basis:\n     Simple basis:\n  FixedBasis([(0.0 + 0.0im)·1 + (-0.577350269189626 - 0.0im)·x₃ + (0.2886751345948128 + 0.5000000000000001im)·x₂ + (0.28867513459481314 - 0.49999999999999994im)·x₁])\n    And entries in a 1×1 VectorizedHermitianMatrix{ComplexF64, Float64, ComplexF64}:\n     0.5000000000000036 + 0.0im"
     },
     "metadata": {},
     "execution_count": 9
    }
   ],
   "cell_type": "code",
   "source": [
    "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",
    "gram_matrix(con_ref)"
   ],
   "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.11.5"
  },
  "kernelspec": {
   "name": "julia-1.11",
   "display_name": "Julia 1.11.5",
   "language": "julia"
  }
 },
 "nbformat": 4
}