examples/Example___coefficients_recursive.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "BgD7IaenqG4B"
      },
      "source": [
        "# Example for `coefficients_recursive.py`"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "We  implemented *recursive* formulas for the characteristic polynomial coefficients from paper [**\"On computing the determinant, other characteristic polynomial coefficients, and inverse in Clifford algebras of arbitrary dimension\"**](https://arxiv.org/abs/2005.04015) (*Theorem 4*)  in `coefficients_recursive.py`. The implementation is valid for `clifford` (symbolic) and `galgebra` (numeric) python packages."
      ],
      "metadata": {
        "id": "7uAxfQrMfL3O"
      }
    },
    {
      "cell_type": "code",
      "execution_count": 15,
      "metadata": {
        "cellView": "form",
        "id": "v1-Xts03y-no",
        "outputId": "559cbb9c-2b8a-4f42-cb9c-f4cce244fdc8",
        "colab": {
          "base_uri": "https://localhost:8080/"
        }
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Requirement already satisfied: clifford in /usr/local/lib/python3.7/dist-packages (1.4.0)\n",
            "Requirement already satisfied: numpy>=1.17 in /usr/local/lib/python3.7/dist-packages (from clifford) (1.21.5)\n",
            "Requirement already satisfied: numba>0.46 in /usr/local/lib/python3.7/dist-packages (from clifford) (0.51.2)\n",
            "Requirement already satisfied: sparse in /usr/local/lib/python3.7/dist-packages (from clifford) (0.13.0)\n",
            "Requirement already satisfied: h5py in /usr/local/lib/python3.7/dist-packages (from clifford) (3.1.0)\n",
            "Requirement already satisfied: scipy in /usr/local/lib/python3.7/dist-packages (from clifford) (1.4.1)\n",
            "Requirement already satisfied: llvmlite<0.35,>=0.34.0.dev0 in /usr/local/lib/python3.7/dist-packages (from numba>0.46->clifford) (0.34.0)\n",
            "Requirement already satisfied: setuptools in /usr/local/lib/python3.7/dist-packages (from numba>0.46->clifford) (57.4.0)\n",
            "Requirement already satisfied: cached-property in /usr/local/lib/python3.7/dist-packages (from h5py->clifford) (1.5.2)\n",
            "Requirement already satisfied: galgebra in /usr/local/lib/python3.7/dist-packages (0.5.0)\n",
            "Requirement already satisfied: sympy in /usr/local/lib/python3.7/dist-packages (from galgebra) (1.7.1)\n",
            "Requirement already satisfied: mpmath>=0.19 in /usr/local/lib/python3.7/dist-packages (from sympy->galgebra) (1.2.1)\n"
          ]
        }
      ],
      "source": [
        "#@title Install and import `clifford` and `galgebra` packages\n",
        "!pip install clifford\n",
        "!pip install galgebra\n",
        "import clifford as cl\n",
        "import sympy\n",
        "from galgebra.ga import Ga"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 16,
      "metadata": {
        "cellView": "form",
        "id": "5CHYUDKZrfff"
      },
      "outputs": [],
      "source": [
        "#@title Import latex printer for notebook\n",
        "from sympy import init_printing\n",
        "def custom_latex_printer(exp,**options):\n",
        "    from google.colab.output._publish import javascript\n",
        "    url = \"https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.3/latest.js?config=default\"\n",
        "    javascript(url=url)\n",
        "    return sympy.printing.latex(exp,**options)\n",
        "init_printing(use_latex=\"mathjax\",latex_printer=custom_latex_printer)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 17,
      "metadata": {
        "cellView": "form",
        "id": "ZlluRKWdfrq8",
        "outputId": "062ecc3e-7299-457b-9b99-e87b003d1035",
        "colab": {
          "base_uri": "https://localhost:8080/"
        }
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "--2022-03-12 13:56:39--  https://raw.githubusercontent.com/kamranuz/clifford_det/main/coefficients_recursive.py\n",
            "Resolving raw.githubusercontent.com (raw.githubusercontent.com)... 185.199.108.133, 185.199.109.133, 185.199.110.133, ...\n",
            "Connecting to raw.githubusercontent.com (raw.githubusercontent.com)|185.199.108.133|:443... connected.\n",
            "HTTP request sent, awaiting response... 200 OK\n",
            "Length: 974 [text/plain]\n",
            "Saving to: ‘coefficients_recursive.py.1’\n",
            "\n",
            "\r          coefficie   0%[                    ]       0  --.-KB/s               \rcoefficients_recurs 100%[===================>]     974  --.-KB/s    in 0s      \n",
            "\n",
            "2022-03-12 13:56:39 (41.2 MB/s) - ‘coefficients_recursive.py.1’ saved [974/974]\n",
            "\n"
          ]
        }
      ],
      "source": [
        "#@title Import recursive formulas for the characteristic polynomial coefficients\n",
        "!wget https://raw.githubusercontent.com/kamranuz/clifford_det/main/coefficients_recursive.py\n",
        "\n",
        "import coefficients_recursive as recur"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "eNLxNO50o4dk"
      },
      "source": [
        "## `clifford`: numerical Geometric Algebra package for python"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 18,
      "metadata": {
        "cellView": "form",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "5Rs7T_sFP5Lq",
        "outputId": "df942896-aefe-46ec-a585-ef90a95eda4e"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "U = -(2^e1) - (4^e2) + (1^e3) - (1^e4) - (1^e5) + (3^e12) - (1^e13) + (3^e14) - (1^e15) + (1^e23) - (4^e24) - (3^e25) - (2^e34) - (2^e35) - (4^e45) - (3^e123) + (5^e124) - (4^e125) + (5^e134) + (1^e135) - (3^e145) + (1^e234) + (2^e245) + (2^e345) + (2^e1234) + (2^e1235) - (2^e1245) + (2^e1345) - (2^e2345) - (2^e12345)\n",
            "CPU times: user 4.27 ms, sys: 0 ns, total: 4.27 ms\n",
            "Wall time: 4.86 ms\n"
          ]
        }
      ],
      "source": [
        "%%time\n",
        "\n",
        "#@title Create $U\\in Cl_{5,0}$ with random integer coefficients\n",
        "l, b = cl.Cl(5,0)\n",
        "locals().update(b)\n",
        "\n",
        "# create random MV with integer coefficients\n",
        "import random\n",
        "u = sum(random.randint(-5,5)*b[i] for i in b)\n",
        "\n",
        "# Or one can create random MV with rational coefficients, but for the rational coefficients the calculation error is too high\n",
        "# u = l.randomMV(min=-5, max=5)\n",
        "\n",
        "# print random MV\n",
        "print('U =',u)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Gdt5m_NihX1a"
      },
      "source": [
        "Using ***recursive*** formulas for the characteristic polynomial coefficients $C_{(i)}$ in the case $n=5$  from `coefficients_recur.py`, we compute the coefficients $C_{(i)}$ and auxiliary terms $U_{(i)}$, $\\quad i=1,2,\\ldots,8$:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 19,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "nZ9ivUM-P5Ls",
        "outputId": "7451b38c-9b79-471c-84aa-6c5a04fc1627"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "CPU times: user 338 ms, sys: 12.3 ms, total: 350 ms\n",
            "Wall time: 413 ms\n"
          ]
        }
      ],
      "source": [
        "%%time\n",
        "\n",
        "# get the generator of the characteristic polynomial coefficients and the corresponding auxiliary elements for n=5, element u\n",
        "# using recursive algorithm\n",
        "coefficients_recur = recur.get_coefficients(u)\n",
        "auxiliary_recur = recur.get_auxiliary(u)\n",
        "\n",
        "# cast coefficients generator to list of coefficients\n",
        "coefficients_recur = list(coefficients_recur)\n",
        "auxiliary_recur = list(auxiliary_recur)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "rCGP1ewVnAkk"
      },
      "source": [
        "Compare the results of ***explicit*** and ***recursive*** formulas in the case $n=5$"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 20,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "I5AR6mICWEfz",
        "outputId": "5e096afe-ba22-463b-f98f-401370ff515c"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Characteristic polynomial coefficients:  [0.0, -468.0, -2576.0, -74662.0, -127008.0, -7273108.0, -4124176.0, -174147201.0]\n"
          ]
        }
      ],
      "source": [
        "print('Characteristic polynomial coefficients: ', coefficients_recur)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "QyhyenR4hJef"
      },
      "source": [
        "Compute the inverse of the element $U$"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 21,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "ty9YwG21P5Lv",
        "outputId": "20ad01fe-4086-40ad-ccbb-bafb9836466f"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "CPU times: user 155 µs, sys: 14 µs, total: 169 µs\n",
            "Wall time: 177 µs\n"
          ]
        }
      ],
      "source": [
        "%%time\n",
        "adj_u = auxiliary_recur[6]-coefficients_recur[6]\n",
        "det_u = coefficients_recur[7]\n",
        "inv_u = adj_u/det_u"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "jDWXxCsqhE26"
      },
      "source": [
        "Check computed inverse of the element $U$  ($U \\cdot U^{-1}$  should be equal to $1$):"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 22,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 0
        },
        "id": "tXCFBnpqWcKS",
        "outputId": "8bab133d-719e-4740-9c67-7e751c8ff646"
      },
      "outputs": [
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "'U * U_inv = 1.0'"
            ],
            "application/vnd.google.colaboratory.intrinsic+json": {
              "type": "string"
            }
          },
          "metadata": {},
          "execution_count": 22
        }
      ],
      "source": [
        "f'U * U_inv = {u * inv_u}'"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "f7-bPrIDpGZ6"
      },
      "source": [
        "## `GAlgebra`: Symbolic Geometric Algebra/Calculus package for SymPy"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 23,
      "metadata": {
        "id": "Xy5VVGzhqRpz",
        "cellView": "form"
      },
      "outputs": [],
      "source": [
        "#@title Create $V\\in Cl_{5,0}$ \n",
        "\n",
        "e12345 = sympy.symbols('1 2 3 4 5', real=True)\n",
        "eta=[1, 1, 1, 1, 1]\n",
        "o3d = Ga('e1 e2 e3 e4 e5', g=eta, coords=e12345)\n",
        "grad = o3d.grad\n",
        "v = o3d.mv('v', 'mv')"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "Using ***explicit*** formulas for the characteristic polynomial coefficients in the case $n=5$ from `coefficients_explicit_n5.py`, we compute the coefficients $C_{(1)}$ and $C_{(2)}$"
      ],
      "metadata": {
        "id": "U-jD43Rsi9qT"
      }
    },
    {
      "cell_type": "code",
      "execution_count": 24,
      "metadata": {
        "id": "lXOqkpXzs-Qp",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "64133545-f50c-4697-e152-cdb31aad3c85"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "CPU times: user 7 µs, sys: 0 ns, total: 7 µs\n",
            "Wall time: 11 µs\n"
          ]
        }
      ],
      "source": [
        "%%time\n",
        "\n",
        "# get the generator of characteristic polynomial coefficients for n=5, element u\n",
        "# using explisit formulas\n",
        "coefs_generator = recur.get_coefficients(v)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 25,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "QHZIUJxbsbP3",
        "outputId": "935c66b9-2265-40d2-af88-0c02cb15c163"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "CPU times: user 361 µs, sys: 0 ns, total: 361 µs\n",
            "Wall time: 368 µs\n"
          ]
        }
      ],
      "source": [
        "%%time\n",
        "c1 = next(coefs_generator)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 26,
      "metadata": {
        "id": "1RHYi6DXuxul",
        "outputId": "53670599-4955-4a74-ce1b-6bd4a1ac247d",
        "colab": {
          "base_uri": "https://localhost:8080/"
        }
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "CPU times: user 21.1 s, sys: 90.3 ms, total: 21.2 s\n",
            "Wall time: 27.5 s\n"
          ]
        }
      ],
      "source": [
        "%%time\n",
        "c2 = next(coefs_generator)"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "print('C_1:')\n",
        "c1,"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 56
        },
        "id": "TCrJzDgkhvBV",
        "outputId": "1ea81871-fd4b-4174-b6e3-0e431c999f6c"
      },
      "execution_count": 27,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "C_1:\n"
          ]
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<IPython.core.display.HTML object>"
            ],
            "text/html": [
              "<script src='https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.3/latest.js?config=default'></script>"
            ]
          },
          "metadata": {}
        },
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "(8.0⋅v,)"
            ],
            "text/latex": "$\\displaystyle \\left( 8.0 v,\\right)$"
          },
          "metadata": {},
          "execution_count": 27
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "print('C_2:')\n",
        "c2,"
      ],
      "metadata": {
        "id": "Brpc0vqHhr9W",
        "outputId": "cbd828e3-b7ec-4b25-9af6-017b9ff86dfc",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 85
        }
      },
      "execution_count": 28,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "C_2:\n"
          ]
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<IPython.core.display.HTML object>"
            ],
            "text/html": [
              "<script src='https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.3/latest.js?config=default'></script>"
            ]
          },
          "metadata": {}
        },
        {
          "output_type": "execute_result",
          "data": {
            "text/plain": [
              "⎛        2         2          2           2            2             2        \n",
              "⎝- 28.0⋅v  + 4.0⋅v¹  - 4.0⋅v¹²  - 4.0⋅v¹²³  + 4.0⋅v¹²³⁴  + 4.0⋅v¹²³⁴⁵  + 4.0⋅v\n",
              "\n",
              "    2           2            2           2          2           2            2\n",
              "¹²³⁵  - 4.0⋅v¹²⁴  + 4.0⋅v¹²⁴⁵  - 4.0⋅v¹²⁵  - 4.0⋅v¹³  - 4.0⋅v¹³⁴  + 4.0⋅v¹³⁴⁵ \n",
              "\n",
              "           2          2           2          2         2          2           \n",
              " - 4.0⋅v¹³⁵  - 4.0⋅v¹⁴  - 4.0⋅v¹⁴⁵  - 4.0⋅v¹⁵  + 4.0⋅v²  - 4.0⋅v²³  - 4.0⋅v²³⁴\n",
              "\n",
              "2            2           2          2           2          2         2        \n",
              "  + 4.0⋅v²³⁴⁵  - 4.0⋅v²³⁵  - 4.0⋅v²⁴  - 4.0⋅v²⁴⁵  - 4.0⋅v²⁵  + 4.0⋅v³  - 4.0⋅v\n",
              "\n",
              "  2           2          2         2          2         2 ⎞\n",
              "³⁴  - 4.0⋅v³⁴⁵  - 4.0⋅v³⁵  + 4.0⋅v⁴  - 4.0⋅v⁴⁵  + 4.0⋅v⁵ ,⎠"
            ],
            "text/latex": "$\\displaystyle \\left( - 28.0 v^{2} + 4.0 \\left(v^{1}\\right)^{2} - 4.0 \\left(v^{12}\\right)^{2} - 4.0 \\left(v^{123}\\right)^{2} + 4.0 \\left(v^{1234}\\right)^{2} + 4.0 \\left(v^{12345}\\right)^{2} + 4.0 \\left(v^{1235}\\right)^{2} - 4.0 \\left(v^{124}\\right)^{2} + 4.0 \\left(v^{1245}\\right)^{2} - 4.0 \\left(v^{125}\\right)^{2} - 4.0 \\left(v^{13}\\right)^{2} - 4.0 \\left(v^{134}\\right)^{2} + 4.0 \\left(v^{1345}\\right)^{2} - 4.0 \\left(v^{135}\\right)^{2} - 4.0 \\left(v^{14}\\right)^{2} - 4.0 \\left(v^{145}\\right)^{2} - 4.0 \\left(v^{15}\\right)^{2} + 4.0 \\left(v^{2}\\right)^{2} - 4.0 \\left(v^{23}\\right)^{2} - 4.0 \\left(v^{234}\\right)^{2} + 4.0 \\left(v^{2345}\\right)^{2} - 4.0 \\left(v^{235}\\right)^{2} - 4.0 \\left(v^{24}\\right)^{2} - 4.0 \\left(v^{245}\\right)^{2} - 4.0 \\left(v^{25}\\right)^{2} + 4.0 \\left(v^{3}\\right)^{2} - 4.0 \\left(v^{34}\\right)^{2} - 4.0 \\left(v^{345}\\right)^{2} - 4.0 \\left(v^{35}\\right)^{2} + 4.0 \\left(v^{4}\\right)^{2} - 4.0 \\left(v^{45}\\right)^{2} + 4.0 \\left(v^{5}\\right)^{2},\\right)$"
          },
          "metadata": {},
          "execution_count": 28
        }
      ]
    }
  ],
  "metadata": {
    "colab": {
      "collapsed_sections": [
        "eNLxNO50o4dk",
        "f7-bPrIDpGZ6"
      ],
      "name": "Example | coefficients_recursive.ipynb",
      "provenance": []
    },
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.9.5"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}