{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Decomposition algorithm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Array{DynamicPolynomials.PolyVar{true},1}:\n",
       " x1\n",
       " x2\n",
       " x3"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using TensorDec\n",
    "using LinearAlgebra\n",
    "X = @ring x1 x2 x3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We want to find a sparse representation of the following series known up to degree 3:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4.0dx1*dx3 + 15.0dx2 + 6.0dx3 + 20.0dx1*dx2 + 6.0dx3^3 + 43.0dx2^2dx3 - 26.0dx1^3 + 129.0dx2^3 + 30.0dx1^2dx2 + 15.0dx2*dx3 + 20.0dx1*dx2*dx3 + 6.0dx1^2 + 6.0dx3^2 + 4.0dx1 + 43.0dx2^2 + 6.0dx1^2dx3 + 4.0dx1*dx3^2 + 72.0dx1*dx2^2 + 15.0dx2*dx3^2 + 6.0"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sigma = dual(6.0 + 4.0*x1 + 15.0*x2 + 6.0*x3 + 6.0*x1^2 + 20.0*x1*x2 + 4.0*x1*x3 + 43.0*x2^2 + 15.0*x2*x3 + 6.0*x3^2 - 26.0*x1^3 + 30.0*x1^2*x2 + 6.0*x1^2*x3 + 72.0*x1*x2^2 + 20.0*x1*x2*x3 + 4.0*x1*x3^2 + 129.0*x2^3 + 43.0*x2^2*x3 + 15.0*x2*x3^2 + 6.0*x3^3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10-element Array{DynamicPolynomials.Monomial{true},1}:\n",
       " 1   \n",
       " x1  \n",
       " x2  \n",
       " x3  \n",
       " x1² \n",
       " x1x2\n",
       " x1x3\n",
       " x2² \n",
       " x2x3\n",
       " x3² "
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "L1 = monoms(X,1)\n",
    "L2 = monoms(X,2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4×10 Array{Float64,2}:\n",
       "  6.0   4.0  15.0   6.0    6.0  20.0   4.0   43.0  15.0   6.0\n",
       "  4.0   6.0  20.0   4.0  -26.0  30.0   6.0   72.0  20.0   4.0\n",
       " 15.0  20.0  43.0  15.0   30.0  72.0  20.0  129.0  43.0  15.0\n",
       "  6.0   4.0  15.0   6.0    6.0  20.0   4.0   43.0  15.0   6.0"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "H = hankel(sigma,L1,L2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The rank of $H_{\\sigma}$ will give us an idea on the dimension of $\\mathcal{A}_\\sigma$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rank(H)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We check that $\\{1, x_1, x_2\\}$ is a basis of $\\mathcal{A}_\\sigma$: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Array{DynamicPolynomials.Monomial{true},1}:\n",
       " 1 \n",
       " x1\n",
       " x2"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "B0 = L1[1:3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Array{Float64,2}:\n",
       "  6.0   4.0  15.0\n",
       "  4.0   6.0  20.0\n",
       " 15.0  20.0  43.0"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "H0 = hankel(sigma, B0, B0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rank(H0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us compute the shifted (truncated) Hankel operators."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Array{Float64,2}:\n",
       "  4.0    6.0  20.0\n",
       "  6.0  -26.0  30.0\n",
       " 20.0   30.0  72.0"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "H1 = hankel(sigma, B0, B0*x1)\n",
    "H2 = hankel(sigma, B0, B0*x2)\n",
    "H3 = hankel(sigma, B0, B0*x3);\n",
    "H  = [H1,H2,H3]\n",
    "H[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Array{Float64,2}:\n",
       "  1.11022e-16   9.14286  -0.571429\n",
       "  1.0           3.85714   1.57143 \n",
       " -1.11022e-16  -4.28571   1.14286 "
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M = [ H0^(-1)*H[i] for i in 1:3 ]\n",
    "M[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The eigenvalues and eigenvectors of $M_{x_1}$ are"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We deduce the operators of multiplication by the variables in the basis $B_0$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Eigen{Float64,Float64,Array{Float64,2},Array{Float64,1}}\n",
       "eigenvalues:\n",
       "3-element Array{Float64,1}:\n",
       " -0.9999999999999991\n",
       "  4.000000000000002 \n",
       "  2.000000000000002 \n",
       "eigenvectors:\n",
       "3×3 Array{Float64,2}:\n",
       "  0.963087  -0.811107  -0.762001\n",
       " -0.120386  -0.324443  -0.127   \n",
       " -0.240772   0.486664   0.635001"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v, E = eigen(M[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The matrices $M_{x_i}$ are diagonal in this basis:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Array{Float64,2}:\n",
       " -1.0          -6.99441e-15  -3.66374e-15\n",
       "  4.21885e-15   4.0          -4.44089e-15\n",
       " -4.66294e-15  -3.9968e-15    2.0        "
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "D = [E^(-1)*M[i]*E for i in 1:3]\n",
    "D[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Array{Float64,2}:\n",
       "  1.0          -4.44089e-16  -1.44329e-15\n",
       "  8.88178e-16   2.0           2.66454e-15\n",
       " -3.55271e-15   2.66454e-15   3.0        "
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "D[2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Array{Float64,2}:\n",
       "  1.0           3.33067e-16  1.11022e-16\n",
       " -9.4369e-16    1.0          6.66134e-16\n",
       "  5.55112e-16  -6.66134e-16  1.0        "
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "D[3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looking at the corresponding terms on the diagonal, we get the coordinates of the points $\\Xi$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Array{Float64,2}:\n",
       " -1.0  4.0  2.0\n",
       "  1.0  2.0  3.0\n",
       "  1.0  1.0  1.0"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Xi = [ D[i][j,j] for i in 1:3, j in 1:3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We normalize the eigenvectors by $v_i \\over v_i(\\xi_i)$ and get the interpolation polynomials at the points $\\xi_i$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Array{DynamicPolynomials.Polynomial{true,Float64},1}:\n",
       " -0.14285714285714324x1 - 0.2857142857142862x2 + 1.142857142857143 \n",
       " 0.28571428571428614x1 - 0.4285714285714279x2 + 0.7142857142857121 \n",
       " -0.14285714285714332x1 + 0.7142857142857134x2 - 0.8571428571428543"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Dg = E'*vcat(fill(1.,1,3), Xi[1:2,:])\n",
    "E = E*Dg^(-1)\n",
    "U = E'*B0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We deduce the weights $w_i=\\sigma(u_i)$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×1 Array{Float64,2}:\n",
       "  1.999999999999992 \n",
       " -1.0000000000000018\n",
       "  5.000000000000002 "
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "w = hankel(sigma, U, [L1[1]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the command `decompose`, we can get directly the same decomposition: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "([-1.0, 5.0, 2.0], [4.0 2.0 -1.0; 2.0 3.0 1.0; 1.0 1.0 1.0])"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "w, Xi = decompose(sigma)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3×3 Array{Float64,2}:\n",
       " 4.0  2.0  -1.0\n",
       " 2.0  3.0   1.0\n",
       " 1.0  1.0   1.0"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Xi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3-element Array{Float64,1}:\n",
       " -1.0000000000000129\n",
       "  5.000000000000011 \n",
       "  1.9999999999999998"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "w"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The series decomposes as $2 \\mathfrak{e}_{(-1,1,1)} + 5 \\mathfrak{e}_{(2,3,1)} - \\mathfrak{e}_{(4,2,1)}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Julia 1.0.0",
   "language": "julia",
   "name": "julia-1.0"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.0.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}