Tutorial-Symbolic_Tensor_Rotation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<script async src=\"https://www.googletagmanager.com/gtag/js?id=UA-59152712-8\"></script>\n",
    "<script>\n",
    "  window.dataLayer = window.dataLayer || [];\n",
    "  function gtag(){dataLayer.push(arguments);}\n",
    "  gtag('js', new Date());\n",
    "\n",
    "  gtag('config', 'UA-59152712-8');\n",
    "</script>\n",
    "\n",
    "# Symbolic Tensor (Quaternion) Rotation\n",
    "\n",
    "## Author: Ken Sible\n",
    "\n",
    "## The following module presents an algorithm for symbolic vector or tensor rotation using SymPy. It validates the algorithm against known results and demonstrates common subexpression elimination for optimized symbolic computation.\n",
    "\n",
    "### NRPy+ Source Code for this module:\n",
    "1. [tensor_rotation.py](../edit/tensor_rotation.py); [\\[**tutorial**\\]](Tutorial-Symbolic_Tensor_Rotation.ipynb) The tensor_rotation.py script will perform symbolic tensor rotation using the following function: `rotate(tensor, axis, angle)`. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='toc'></a>\n",
    "\n",
    "# Table of Contents\n",
    "$$\\label{toc}$$\n",
    "\n",
    "0. [Step 0](#prelim): Derivation of quaternion rotation from matrix rotation using [linear algebra](https://en.wikipedia.org/wiki/Linear_algebra) and [ring theory](https://en.wikipedia.org/wiki/Ring_theory)\n",
    "1. [Step 1](#algorithm): Discussion of the tensor rotation algorithm using the [SymPy](https://www.sympy.org) package for symbolic manipulation </header_section>\n",
    "1. [Step 2](#validation): Validation and demonstration of the tensor rotation algorithm, including common subexpression elimination\n",
    "    1. [Step 2.a](#axisrotation): Verification of three-dimensional rotation about a coordinate axis using an equivalent rotation matrix\n",
    "    1. [Step 2.b](#symtensor): Verification of symbolic tensor rotation and the invariance of argument permutation for a [symmetric tensor](https://en.wikipedia.org/wiki/Symmetric_tensor)\n",
    "    1. [Step 2.c](#cse): Demonstration of [common subexpression elimination](https://en.wikipedia.org/wiki/Common_subexpression_elimination) applied to a rotated symbolic matrix\n",
    "1. [Step 3](#latex_pdf_output): Output this notebook to $\\LaTeX$-formatted PDF file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='prelim'></a>\n",
    "\n",
    "# Mathematical Background (Optional): Quaternion Rotation \\[Back to [top](#toc)\\]\n",
    "$$\\label{prelim}$$ \n",
    "\n",
    "Let $\\vec{v}$ denote a vector in $\\mathbb{R}^2$. We recall from linear algebra that $\\vec{v}$ rotated through an angle $\\theta$ about the x-axis, denoted $\\vec{v}'$, has the following matrix formula\n",
    "\n",
    "$$\\vec{v}'=\\begin{bmatrix}\\cos\\theta & -\\sin\\theta \\\\ \\sin\\theta & \\cos\\theta \\end{bmatrix}\\vec{v}.$$\n",
    "\n",
    "Let $z=a+bi\\in\\mathbb{C}$ for some $a,b\\in\\mathbb{R}$. Consider the corresponding (or [isomorphic](https://en.wikipedia.org/wiki/Isomorphism)) vector $\\vec{z}=(a,b)\\in\\mathbb{R}^2$. We observe from the rotation formula that $\\vec{v}'=a(\\cos\\theta,\\sin\\theta)+b(-\\sin\\theta,\\cos\\theta)$ after expanding the matrix product. Let $w=c+di\\in\\mathbb{C}$ for some $c,d\\in\\mathbb{R}$. We recall the definition of the complex product as $zw=(ac-bd)+i(ad+bc)$ where $i^2=-1$. Hence, $z'=(a+bi)(\\cos\\theta+i\\sin\\theta)=(a+bi)e^{i\\theta}$ after comparing the rotated vector $\\vec{v}'$ with the complex product as defined above. Therefore, the following compact formula arises for vector rotation (given the correspondence between $\\mathbb{C}$ and $\\mathbb{R}^2$):\n",
    "\n",
    "$$z'=e^{i\\theta}z.$$\n",
    "\n",
    "However, for vector rotation in three-dimensional space, we curiously require a four-dimensional, non-commutative extension of the complex number system. The following discussion will provide an overview of the connection with vector rotation. For further reading, we suggest the document ([quaternion_rotation](http://graphics.stanford.edu/courses/cs348a-17-winter/Papers/quaternion.pdf)).\n",
    "\n",
    "$$\\mathcal{H}=\\{a+bi+cj+dk:a,b,c,d\\in\\mathbb{R}\\text{ and }i^2=j^2=k^2=ijk=-1\\}$$\n",
    "\n",
    "Consider the special case of rotating a vector $\\vec{v}\\in\\mathbb{R}^3$ through an angle $\\theta$ about a normalized rotation axis $\\vec{n}$ perpendicular to the vector $\\vec{v}$. We typically decompose a quaternion into a scalar and vector component whenever performing vector rotation, such as the decomposition of $q=a+bi+cj+dk$ into $q=(a,\\vec{w})$ where $\\vec{w}=(b,c,d)$. For future convenience, we define the following quaternions for vector rotation: $v=(0,\\vec{v})$, $v'=(0,\\vec{v}')$, and $n=(0,\\vec{n})$.\n",
    "\n",
    "From the fundamental formula of quaternion algebra, $i^2=j^2=k^2=ijk=-1$, we could derive quaternion multiplication, known as the [Hamilton product](https://en.wikipedia.org/wiki/Quaternion#Hamilton_product). Let $q_1,q_2\\in\\mathcal{H}$ where both $q_1$ and $q_2$ have zero scalar component. The Hamilton product of $q_1$ and $q_2$, after some straightforward verification, has the form $q_1q_2=(-\\vec{q}_1\\cdot\\vec{q}_2,\\vec{q}_1\\times\\vec{q}_2)$. Hence, $nv=(0,\\vec{n}\\times\\vec{v})$ since $\\vec{n}$ and $\\vec{v}$ are orthogonal. We observe that the projection of the rotated vector $\\vec{v}'$ onto $\\vec{v}$ is $\\cos\\theta\\,\\vec{v}$ and the projection of $\\vec{v}'$ onto $\\vec{n}\\times\\vec{v}$ is $\\sin\\theta\\,(\\vec{n}\\times\\vec{v})$, and hence $\\vec{v}'=\\cos\\theta\\,\\vec{v}+\\sin\\theta\\,(\\vec{n}\\times\\vec{v})$. We define the quaternion exponential as $e^{n\\theta}=\\cos\\theta+n\\sin\\theta$, analogous to the complex exponential. Finally, we identify the formula for vector rotation in $\\mathbb{R}^3$ by comparing the rotated vector $\\vec{v}'$ with the Hamilton product:\n",
    "\n",
    "$$v'=(\\cos\\theta+n\\sin\\theta)v=e^{n\\theta}v$$\n",
    "\n",
    "The tensor rotation algorithm defined by the function `rotate(tensor, axis, angle)` in `tensor_rotation.py` does general vector and tensor rotation in $\\mathbb{R}^3$ about an arbitrary rotation axis, not necessarily orthogonal to the original vector or tensor. For general three-dimensional rotation, we define the rotation quaternion as $q=e^{n(\\theta/2)}$ and the conjugate of $q$ as $q^*=e^{-n(\\theta/2)}$. The quaternion rotation operator has the following form for general vector and tensor rotation.  \n",
    "\n",
    "$$\\mathcal{L}[v]=qvq^*,\\,\\,\\mathcal{L}[M]=(q(qMq^*)^\\text{T}q^*)^\\text{T}$$\n",
    "\n",
    "We should remark that quaternion-matrix multiplication is defined as column-wise quaternion multiplication [(source)](https://people.dsv.su.se/~miko1432/rb/Rotations%20of%20Tensors%20using%20Quaternions%20v0.3.pdf). Furthermore, we claim that $vq^*=qv$ whenever the rotation axis $\\vec{n}$ and the vector $\\vec{v}$ are perpendicular (straightforward verification), which will recover the rotation formula for the special case."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='algorithm'></a>\n",
    "\n",
    "# Step 1: The Tensor Rotation Algorithm \\[Back to [top](#toc)\\]\n",
    "$$\\label{algorithm}$$\n",
    "\n",
    "### Algorithm Pseudocode\n",
    "```\n",
    "function rotate(tensor, axis, angle):\n",
    "    initialize quaternion q from rotation axis and angle\n",
    "    if tensor is a nested list (i.e. matrix)\n",
    "        convert tensor to a symbolic matrix\n",
    "        if not (size of symbolic matrix is (3, 3))\n",
    "            throw error for invalid matrix size\n",
    "        end if\n",
    "        initialize empty vector M\n",
    "        for column in tensor\n",
    "            append quaternion(column) onto M\n",
    "        end for\n",
    "        M = q * M *(conjugate of q) // rotate each column\n",
    "        replace each column of tensor with the vector part\n",
    "            of the associated column quaternion in M\n",
    "        initialize empty vector M\n",
    "        for row in tensor\n",
    "            append quaternion(row) onto M\n",
    "        end for\n",
    "        M = q * M *(conjugate of q) // rotate each row\n",
    "        replace each row of tensor with the vector part\n",
    "            of the associated row quaternion in M\n",
    "        convert tensor to a nested list\n",
    "        return tensor\n",
    "    else if tensor is a list (i.e. vector)\n",
    "        if not (length of vector is 3):\n",
    "            throw error for invalid vector length\n",
    "        v = q * quaternion(tensor) * (conjugate of q)\n",
    "        replace each element of tensor with the vector part\n",
    "            of the associated vector quaternion v\n",
    "        return tensor\n",
    "    end if\n",
    "    throw error for unsupported tensor type\n",
    "end function\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:25:31.671923Z",
     "iopub.status.busy": "2021-03-07T17:25:31.670621Z",
     "iopub.status.idle": "2021-03-07T17:25:31.988467Z",
     "shell.execute_reply": "2021-03-07T17:25:31.987746Z"
    }
   },
   "outputs": [],
   "source": [
    "\"\"\" Symbolic Tensor (Quaternion) Rotation\n",
    "\n",
    "The following script will perform symbolic tensor rotation using quaternions.\n",
    "\"\"\"\n",
    "\n",
    "from sympy import Quaternion as quat\n",
    "from sympy import Matrix\n",
    "from sympy.functions import transpose\n",
    "\n",
    "# Input:  tensor = 3-vector or (3x3)-matrix\n",
    "#         axis   = rotation axis (normal 3-vector)\n",
    "#         angle  = rotation angle (in radians)\n",
    "# Output: rotated tensor (of original type)\n",
    "def rotate(tensor, axis, angle):\n",
    "    # Quaternion-Matrix Multiplication\n",
    "    def mul(*args):\n",
    "        if isinstance(args[0], list):\n",
    "            q, M = args[1], args[0]\n",
    "            for i, col in enumerate(M):\n",
    "                M[i] = col * q\n",
    "        else:\n",
    "            q, M = args[0], args[1]\n",
    "            for i, col in enumerate(M):\n",
    "                M[i] = q * col\n",
    "        return M\n",
    "    # Rotation Quaternion (Axis, Angle)\n",
    "    q = quat.from_axis_angle(axis, angle)\n",
    "    if isinstance(tensor[0], list):\n",
    "        tensor = Matrix(tensor)\n",
    "        if tensor.shape != (3, 3):\n",
    "            raise Exception('Invalid Matrix Size')\n",
    "        # Rotation Formula: M' = (q.(q.M.q*)^T.q*)^T\n",
    "        M = [quat(0, *tensor[:, i]) for i in range(tensor.shape[1])]\n",
    "        M = mul(q, mul(M, q.conjugate()))\n",
    "        for i in range(tensor.shape[1]):\n",
    "            tensor[:, i] = [M[i].b, M[i].c, M[i].d]\n",
    "        M = [quat(0, *tensor[i, :]) for i in range(tensor.shape[0])]\n",
    "        M = mul(q, mul(M, q.conjugate()))\n",
    "        for i in range(tensor.shape[0]):\n",
    "            tensor[i, :] = [[M[i].b, M[i].c, M[i].d]]\n",
    "        return tensor.tolist()\n",
    "    if isinstance(tensor, list):\n",
    "        if len(tensor) != 3:\n",
    "            raise Exception('Invalid Vector Length')\n",
    "        # Rotation Formula: v' = q.v.q*\n",
    "        v = q * quat(0, *tensor) * q.conjugate()\n",
    "        return [v.b, v.c, v.d]\n",
    "    raise Exception('Unsupported Tensor Type')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='validation'></a>\n",
    "\n",
    "# Step 2: Validation and Demonstration \\[Back to [top](#toc)\\]\n",
    "$$\\label{validation}$$\n",
    "\n",
    "In the following section, we demonstrate the rotation algorithm, specifically for symbolic tensor rotation, and validate that the numeric or symbolic output from the algorithm does agree with a known result, usually obtained from a rotation matrix. The format of each code cell has the structure: generate expected result from a rotation matrix or NRPy+, generate received result from the rotation algorithm, and assert that these are both equivalent."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='axisrotation'></a>\n",
    "\n",
    "## Step 2.a: Verifying Coordinate Axis Rotation \\[Back to [top](#toc)\\]\n",
    "$$\\label{axisrotation}$$\n",
    "\n",
    "We recall that any general three-dimensional rotation can be expressed as the composition of a rotation about each coordinate axis (see [rotation theorem](https://en.wikipedia.org/wiki/Euler%27s_rotation_theorem)). Therefore, we validate our rotation algorithm using only rotations about each coordinate axis rather than about an arbitrary axis in three-dimensional space. Consider the following vector $\\vec{v}$ and matrix $M$ defined below using the SymPy package."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:25:32.001439Z",
     "iopub.status.busy": "2021-03-07T17:25:32.000750Z",
     "iopub.status.idle": "2021-03-07T17:25:32.004048Z",
     "shell.execute_reply": "2021-03-07T17:25:32.004523Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\vec{v}=\\left[\\begin{matrix}1\\\\0\\\\1\\end{matrix}\\right],\\,M=\\left[\\begin{matrix}1 & 2 & 1\\\\0 & 1 & 0\\\\2 & 1 & 2\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sympy.matrices import rot_axis1, rot_axis2, rot_axis3\n",
    "from sympy import latex, pi\n",
    "from IPython.display import Math\n",
    "\n",
    "v, angle = [1, 0, 1], pi/2\n",
    "M = [[1, 2, 1], [0, 1, 0], [2, 1, 2]]\n",
    "Math('\\\\vec{v}=%s,\\\\,M=%s' % (latex(Matrix(v)), latex(Matrix(M))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We further recall that for any rotation matrix $R$ and vector $\\vec{v}$, the rotated vector $\\vec{v}'$ has the formula $\\vec{v}'=R\\vec{v}$ and the rotated matrix $M'$ has the formula $M'=RMR^\\text{-1}$, where $R^{-1}=R^\\text{T}$ since every rotation matrix is an [orthogonal matrix](https://en.wikipedia.org/wiki/Orthogonal_matrix)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:25:32.011228Z",
     "iopub.status.busy": "2021-03-07T17:25:32.010512Z",
     "iopub.status.idle": "2021-03-07T17:25:32.063061Z",
     "shell.execute_reply": "2021-03-07T17:25:32.063570Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\vec{v}'=\\left[\\begin{matrix}1\\\\-1\\\\0\\end{matrix}\\right],\\,M'=\\left[\\begin{matrix}1 & -1 & 2\\\\-2 & 2 & -1\\\\0 & 0 & 1\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# vector rotation about x-axis\n",
    "expected = rot_axis1(-angle) * Matrix(v)\n",
    "received = Matrix(rotate(v, [1, 0, 0], angle))\n",
    "assert expected == received; v_ = received\n",
    "\n",
    "# matrix rotation about x-axis\n",
    "expected = rot_axis1(-angle) * Matrix(M) * transpose(rot_axis1(-angle))\n",
    "received = Matrix(rotate(M, [1, 0, 0], angle))\n",
    "assert expected == received; M_ = received\n",
    "\n",
    "Math('\\\\vec{v}\\'=%s,\\\\,M\\'=%s' % (latex(Matrix(v_)), latex(Matrix(M_))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:25:32.106209Z",
     "iopub.status.busy": "2021-03-07T17:25:32.105308Z",
     "iopub.status.idle": "2021-03-07T17:25:32.108975Z",
     "shell.execute_reply": "2021-03-07T17:25:32.108444Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\vec{v}'=\\left[\\begin{matrix}1\\\\0\\\\-1\\end{matrix}\\right],\\,M'=\\left[\\begin{matrix}2 & 1 & -2\\\\0 & 1 & 0\\\\-1 & -2 & 1\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# vector rotation about y-axis\n",
    "expected = rot_axis2(-angle) * Matrix(v)\n",
    "received = Matrix(rotate(v, [0, 1, 0], angle))\n",
    "assert expected == received; v_ = received\n",
    "\n",
    "# matrix rotation about y-axis\n",
    "expected = rot_axis2(-angle) * Matrix(M) * transpose(rot_axis2(-angle))\n",
    "received = Matrix(rotate(M, [0, 1, 0], angle))\n",
    "assert expected == received; M_ = received\n",
    "\n",
    "Math('\\\\vec{v}\\'=%s,\\\\,M\\'=%s' % (latex(Matrix(v_)), latex(Matrix(M_))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:25:32.120802Z",
     "iopub.status.busy": "2021-03-07T17:25:32.119842Z",
     "iopub.status.idle": "2021-03-07T17:25:32.124276Z",
     "shell.execute_reply": "2021-03-07T17:25:32.124799Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\vec{v}'=\\left[\\begin{matrix}0\\\\1\\\\1\\end{matrix}\\right],\\,M'=\\left[\\begin{matrix}1 & 0 & 0\\\\-2 & 1 & 1\\\\-1 & 2 & 2\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# vector rotation about z-axis\n",
    "expected = rot_axis3(-angle) * Matrix(v)\n",
    "received = Matrix(rotate(v, [0, 0, 1], angle))\n",
    "assert expected == received; v_ = received\n",
    "\n",
    "# matrix rotation about z-axis\n",
    "expected = rot_axis3(-angle) * Matrix(M) * transpose(rot_axis3(-angle))\n",
    "received = Matrix(rotate(M, [0, 0, 1], angle))\n",
    "assert expected == received; M_ = received\n",
    "\n",
    "Math('\\\\vec{v}\\'=%s,\\\\,M\\'=%s' % (latex(Matrix(v_)), latex(Matrix(M_))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='symtensor'></a>\n",
    "\n",
    "## Step 2.b: Verifying Symbolic Tensor Rotation \\[Back to [top](#toc)\\]\n",
    "$$\\label{symtensor}$$\n",
    "\n",
    "The rotation algorithm does support symbolic rotation, as shown below with the second rank, symmetric tensor $h^{\\mu\\nu}$ rotated about the 4-vector $v^\\mu$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:25:32.130688Z",
     "iopub.status.busy": "2021-03-07T17:25:32.129971Z",
     "iopub.status.idle": "2021-03-07T17:25:32.315449Z",
     "shell.execute_reply": "2021-03-07T17:25:32.314789Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle hUU=\\left[\\begin{matrix}hUU_{00} & hUU_{01} & hUU_{02}\\\\hUU_{01} & hUU_{11} & hUU_{12}\\\\hUU_{02} & hUU_{12} & hUU_{22}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sympy import symbols, simplify, cse\n",
    "import indexedexp as ixp # NRPy+: symbolic indexed expressions\n",
    "\n",
    "angle = symbols('theta', real=True)\n",
    "vU    = ixp.declarerank1(\"vU\")\n",
    "hUU   = ixp.declarerank2(\"hUU\", \"sym01\")\n",
    "rotatedhUU, N = rotate(hUU, vU, angle), len(hUU)\n",
    "\n",
    "Math('hUU=%s' % latex(Matrix(hUU)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We demonstrate that a completely symbolic rotation applied to the tensor $h^{\\mu\\nu}$ does preserve the index symmetry ($h^{\\mu\\nu}=h^{\\nu\\mu}$)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:25:32.338764Z",
     "iopub.status.busy": "2021-03-07T17:25:32.323324Z",
     "iopub.status.idle": "2021-03-07T17:25:52.860403Z",
     "shell.execute_reply": "2021-03-07T17:25:52.859907Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Assertion Passed: rotatedhUU[1][0] == rotatedhUU[0][1]\n",
      "Assertion Passed: rotatedhUU[2][0] == rotatedhUU[0][2]\n",
      "Assertion Passed: rotatedhUU[2][1] == rotatedhUU[1][2]\n"
     ]
    }
   ],
   "source": [
    "for i in range(N):\n",
    "    for j in range(N):\n",
    "        if j >= i: continue\n",
    "        assert simplify(rotatedhUU[i][j] - rotatedhUU[j][i]) == 0\n",
    "        print('Assertion Passed: rotatedhUU[{i}][{j}] == rotatedhUU[{j}][{i}]'.format(i=i, j=j))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='cse'></a>\n",
    "\n",
    "## Step 2.c: Common Subexpression Elimination \\[Back to [top](#toc)\\]\n",
    "$$\\label{cse}$$\n",
    "\n",
    "If the rotation algorithm is given any symbolic input, then the resulting expression will support common subexpression elimination."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:25:52.937291Z",
     "iopub.status.busy": "2021-03-07T17:25:52.901641Z",
     "iopub.status.idle": "2021-03-07T17:25:53.027579Z",
     "shell.execute_reply": "2021-03-07T17:25:53.027069Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\text{cse}(hUU)=\\left[\\begin{matrix}x_{1} x_{33} + x_{3} x_{34} - x_{36} x_{5} + x_{37} x_{9} & x_{1} x_{36} + x_{3} x_{37} + x_{33} x_{5} - x_{34} x_{9} & x_{1} x_{34} - x_{3} x_{33} + x_{36} x_{9} + x_{37} x_{5}\\\\x_{1} x_{41} + x_{3} x_{42} - x_{44} x_{5} + x_{45} x_{9} & x_{1} x_{44} + x_{3} x_{45} + x_{41} x_{5} - x_{42} x_{9} & x_{1} x_{42} - x_{3} x_{41} + x_{44} x_{9} + x_{45} x_{5}\\\\x_{1} x_{49} + x_{3} x_{50} - x_{5} x_{52} + x_{53} x_{9} & x_{1} x_{52} + x_{3} x_{53} + x_{49} x_{5} - x_{50} x_{9} & x_{1} x_{50} - x_{3} x_{49} + x_{5} x_{53} + x_{52} x_{9}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "<IPython.core.display.Math object>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cse_rotatedhUU = cse(Matrix(rotatedhUU))[1][0]\n",
    "Math('\\\\text{cse}(hUU)=%s' % latex(Matrix(cse_rotatedhUU)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='latex_pdf_output'></a>\n",
    "\n",
    "# Step 3: Output this notebook to $\\LaTeX$-formatted PDF file \\[Back to [top](#toc)\\]\n",
    "$$\\label{latex_pdf_output}$$\n",
    "\n",
    "The following code cell converts this Jupyter notebook into a proper, clickable $\\LaTeX$-formatted PDF file. After the cell is successfully run, the generated PDF may be found in the root NRPy+ tutorial directory, with filename\n",
    "[Tutorial-Symbolic_Tensor_Rotation.pdf](Tutorial-Symbolic_Tensor_Rotation.pdf). (Note that clicking on this link may not work; you may need to open the PDF file through another means.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-03-07T17:25:53.032545Z",
     "iopub.status.busy": "2021-03-07T17:25:53.031471Z",
     "iopub.status.idle": "2021-03-07T17:25:56.530327Z",
     "shell.execute_reply": "2021-03-07T17:25:56.529096Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Created Tutorial-Symbolic_Tensor_Rotation.tex, and compiled LaTeX file to\n",
      "    PDF file Tutorial-Symbolic_Tensor_Rotation.pdf\n"
     ]
    }
   ],
   "source": [
    "import cmdline_helper as cmd    # NRPy+: Multi-platform Python command-line interface\n",
    "cmd.output_Jupyter_notebook_to_LaTeXed_PDF(\"Tutorial-Symbolic_Tensor_Rotation\")"
   ]
  }
 ],
 "metadata": {
  "file_extension": ".py",
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.1"
  },
  "mimetype": "text/x-python",
  "name": "python",
  "npconvert_exporter": "python",
  "pygments_lexer": "ipython3",
  "version": 3
 },
 "nbformat": 4,
 "nbformat_minor": 4
}