{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import gpt as g"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Lecture 2: Linear algebra with lattices"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We start by creating a $4^4$ double precision grid."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "grid = g.grid([4, 4, 4, 4], g.double)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we create a parallel pseudorandom number generator and two color vectors from a complex normal distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "GPT :       6.818323 s : Initializing gpt.random(seed for rng,vectorized_ranlux24_389_64) took 0.00201893 s\n"
     ]
    }
   ],
   "source": [
    "rng = g.random(\"seed for rng\")\n",
    "v1 = g.vcolor(grid)\n",
    "v2 = g.vcolor(grid)\n",
    "rng.cnormal([v1,v2]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We inspect a one-dimensional slice of $v_1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.58527248-0.85360387j, -1.580872  -1.50785674j,\n",
       "        -1.04090035+0.40186692j],\n",
       "       [-0.23550862-0.90891595j,  0.73888871+0.21903038j,\n",
       "         1.10074976+1.52119033j],\n",
       "       [ 0.79255516+0.64832885j, -0.22812566+0.54092838j,\n",
       "         1.17416062-0.72316771j],\n",
       "       [-0.81076797+0.6141707j ,  0.6581139 -0.02166552j,\n",
       "         0.17446881+0.19891321j]])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v1[0,0,0,:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Expressions\n",
    "Let us now take a linear combination of $\\frac12 v_1 + 3 v_2$.  Note that there is a difference between abstract expressions and evaluating them to lattice fields."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "GPT :      28.000815 s : This is an expression:  + ((0.5+0j))*lattice(ot_vector_color(3),double) + ((3+0j))*lattice(ot_vector_color(3),double)\n",
      "GPT :      28.006040 s : And this is the coordinate (0,0,0,0) of the evaluated lattice object: tensor([ 0.09215363+4.58520759j  0.53058119-4.0223991j  -0.79159505+2.5800456j ],ot_vector_color(3))\n"
     ]
    }
   ],
   "source": [
    "expr = 1./2. * v1 + 3. * v2\n",
    "\n",
    "g.message(\"This is an expression:\", expr)\n",
    "\n",
    "result = g.eval(expr)\n",
    "\n",
    "g.message(\"And this is the coordinate (0,0,0,0) of the evaluated lattice object:\", result[0,0,0,0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since the \"=\" operator in Python cannot be overloaded and always assigns the Python object on the right side to the symbol to left, GPT uses the operator \"@=\" to assign the result of an expression to a lattice field."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "GPT :      39.886252 s : Difference between g.eval(...) and @= operation: 0.0\n"
     ]
    }
   ],
   "source": [
    "result2 = g.lattice(v1) # creates a lattice of same type than v1 without initializing the values\n",
    "result2 @= expr\n",
    "\n",
    "g.message(\"Difference between g.eval(...) and @= operation:\", g.norm2(result2 - result))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is also a short-hand notation for g.eval(...), which is just g(...)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "GPT :     138.149730 s : Difference between g.eval(...) and g(...): 0.0\n"
     ]
    }
   ],
   "source": [
    "result3 = g(expr)\n",
    "\n",
    "g.message(f\"Difference between g.eval(...) and g(...):\", g.norm2(result3 - result))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also work on lists of lattice objects at the same time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "GPT :     440.245933 s : Tests: 4.406174692996534e-29 5.068846293762706e-29\n"
     ]
    }
   ],
   "source": [
    "vv_a = [v1, v2]\n",
    "vv_b = [v2, v1]\n",
    "\n",
    "res_list = g(2*g.expr(vv_a) - 3*g.expr(vv_b))\n",
    "g.message(\n",
    "    \"Tests:\",\n",
    "    g.norm2(res_list[0] - 2*v1 + 3*v2),\n",
    "    g.norm2(res_list[1] - 2*v2 + 3*v1)\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where we needed to wrap the list of lattice object in an explicit g.expr statement.\n",
    "\n",
    "It is also useful to know of the increment operator:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "result += v1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, let us construct site-local inner and outer products from our vectors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "GPT :       0.479349 s : Origin of inner_product field: (-0.14170995506280581+0.7288503091897771j)\n",
      "GPT :       0.482812 s : Origin of outer_product field: tensor([[-1.50115922+0.86831104j  0.67227476-1.01352379j -0.6240419 +0.54129306j]\n",
      "                       :  [-2.72189862+2.44771251j  0.94667551-2.38631285j -1.05290498+1.38997313j]\n",
      "                       :  [ 0.53787765+1.79054559j -0.89617918-0.95709307j  0.41277376+0.7891515j ]],ot_matrix_su_n_fundamental_group(3))\n"
     ]
    }
   ],
   "source": [
    "local_inner_product = g.complex(grid)\n",
    "local_inner_product @= g.adj(v1) * v2\n",
    "\n",
    "local_outer_product = g.mcolor(grid)\n",
    "local_outer_product @= v1 * g.adj(v2)\n",
    "\n",
    "g.message(\"Origin of inner_product field:\", local_inner_product[0,0,0,0])\n",
    "\n",
    "g.message(\"Origin of outer_product field:\", local_outer_product[0,0,0,0])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "GPT :       0.496987 s : v1 . v2 (-44.827357722843956-52.851620669240525j)\n",
      "GPT :       0.499137 s : sum(local_inner_product) (-44.82735772284396-52.851620669240525j)\n",
      "GPT :       0.500471 s : v1 . v1 (1592.0538142528967-2.14659970599242e-15j)\n",
      "GPT :       0.501776 s : norm2(v1) 1592.0538142528967\n"
     ]
    }
   ],
   "source": [
    "v12 = g.inner_product(v1, v2)\n",
    "v11 = g.inner_product(v1, v1)\n",
    "n1 = g.norm2(v1)\n",
    "\n",
    "g.message(\"v1 . v2\", v12)\n",
    "g.message(\"sum(local_inner_product)\", g.sum(local_inner_product))\n",
    "\n",
    "g.message(\"v1 . v1\", v11)\n",
    "g.message(\"norm2(v1)\", n1)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We move on to SU(3) matrices now.  g.mcolor is short-hand for the QCD gauge group, but GPT has general SU(N) and U(1) implemented.  We initialize a SU(3) matrix with a random near-unit element and then compute its site-wise determinant and inverse and check the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "GPT :       0.526540 s : Origin of V: tensor([[ 0.99992903-0.00286107j -0.00411134-0.01041353j  0.00277436+0.00083957j]\n",
      "                       :  [ 0.00410881-0.01043508j  0.99991181+0.00142932j -0.00571668+0.00398393j]\n",
      "                       :  [-0.00270793+0.00080031j  0.00572503+0.00402551j  0.9999705 +0.00143171j]],ot_matrix_su_n_fundamental_group(3))\n"
     ]
    }
   ],
   "source": [
    "V = g.mcolor(grid)\n",
    "rng.normal_element(V, scale=0.01)\n",
    "\n",
    "g.message(\"Origin of V:\", V[0,0,0,0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "GPT :       0.562489 s : Slice of local determinant: [[1.-3.14017533e-20j]\n",
      "                       :  [1.-6.78225009e-19j]\n",
      "                       :  [1.-2.12880792e-20j]\n",
      "                       :  [1.+1.37204908e-18j]]\n"
     ]
    }
   ],
   "source": [
    "det_V = g.matrix.det(V)\n",
    "\n",
    "g.message(\"Slice of local determinant:\", det_V[0,0,0,:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "GPT :       0.575219 s : Difference between matrix inverse and adjoint for unitary matrix: 1.849405542690569e-29\n"
     ]
    }
   ],
   "source": [
    "inv_V = g.matrix.inv(V)\n",
    "\n",
    "g.message(\"Difference between matrix inverse and adjoint for unitary matrix:\", g.norm2(inv_V - g.adj(V)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also check that the logarithm of the matrix is anti-Hermitian."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.1647695871037737e-29"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "logV = g.matrix.log(V)\n",
    "g.norm2(logV + g.adj(logV))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are also component-wise operations, see, e.g., the cosine:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "tensor([[0.54036423+2.40739705e-03j 1.00004577-4.28141474e-05j\n",
       "  0.9999965 -2.32927134e-06j]\n",
       " [1.000046  +4.28763992e-05j 0.54037706-1.20266577e-03j\n",
       "  0.9999916 +2.27747790e-05j]\n",
       " [0.99999665+2.16718170e-06j 0.99999171-2.30460831e-05j\n",
       "  0.54032768-1.20471735e-03j]],ot_matrix_su_n_fundamental_group(3))"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "g.component.cos(V)[0,0,0,0]"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (default)",
   "language": "python",
   "name": "python3-default"
  },
  "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.13.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}