learnLyncs.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "cc84030b",
   "metadata": {},
   "source": [
    "# Learn Lyncs-API"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "1cede082",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import sys\n",
    "import numpy as np  # type: ignore\n",
    "import lyncs_io as lio  # type: ignore\n",
    "from typing import List"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7da040a0",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\"><b>Warning:</b>If you get the warning <code><FONT COLOR=\"#db5048\">ModuleNotFoundError</code><code>: No module named 'lyncs_io'</code>, you need to change the kernel via Kernel > Change kernel > lyncs.</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "52bfdb32",
   "metadata": {},
   "outputs": [],
   "source": [
    "import plaquette as RPlaq"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e9c73de4",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "\n",
    "<a href=\"https://github.com/Lyncs-API\">Lyncs-API</a> is a Python API for Lattice QCD applications () created by Simone Bacchio, Christodoulos Stylianou and Alexandros Angeli.\n",
    "\n",
    "One such API is <a href=\"https://github.com/Lyncs-API/lyncs.io\">lyncs.io</a>, a suite of I/O functions which allow for quick and simple interfacing with several common file formats, including lattice-oriented formats such as `.lime` and `.openqcd`."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "84b8825b",
   "metadata": {},
   "source": [
    "## Accessing the data\n",
    "\n",
    "Let us use some example gaugefield files to showcase `lyncs`. We use $8\\times24^3$ SU(3) $N_f = 2 + 1$ gaugefields in $3+1d$. These otherwise have the same parameters as FASTSUM's Generation 2 ensembles \\[[1](https://arxiv.org/abs/1412.6411), [2](https://arxiv.org/abs/2007.04188)\\]."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "c35ee30f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# the location, name and ID number of the gaugefields\n",
    "gfDir = 'confs'\n",
    "gfName = 'Gen2_8x24n'\n",
    "gfIDs = [7,8,9]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4b3b0167",
   "metadata": {},
   "source": [
    "Using the `head()` function, we can read the header data stored in the file and access the information in the form of a `dict` object.\n",
    "\n",
    "Here we see that our example gaugefield files are arranged in the shape $N_t \\times N_s^3 \\times N_d \\times N_c^2$. The `dtype` key shows that the datatype is `'<c16'` or (little-endian) double-precision complex. `'_offset'` is the number of bytes in the header and `'plaq'` is the average value of the spatial and temporal plaquette."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "277510c3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "confs/Gen2_8x24n7: {'shape': (8, 24, 24, 24, 4, 3, 3), 'dtype': '<c16', '_offset': 24, 'plaq': 1.6265985010264397}\n",
      "confs/Gen2_8x24n8: {'shape': (8, 24, 24, 24, 4, 3, 3), 'dtype': '<c16', '_offset': 24, 'plaq': 1.6235420123884416}\n",
      "confs/Gen2_8x24n9: {'shape': (8, 24, 24, 24, 4, 3, 3), 'dtype': '<c16', '_offset': 24, 'plaq': 1.6244340856720185}\n"
     ]
    }
   ],
   "source": [
    "# We can probe the header data of the gaugefield files\n",
    "# Loop over each ID\n",
    "for iid in gfIDs:\n",
    "    gfFile = os.path.join(gfDir,f'{gfName}{iid}')\n",
    "    # Read and print header\n",
    "    print(f\"{gfFile}:\", lio.head(gfFile, format='openqcd'))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "917c999e",
   "metadata": {},
   "source": [
    "The `load()` function returns a `numpy.ndarray` containing the gaugefields. Let's load some example data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "0e29b38d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the gaugefields\n",
    "# Make a list of gaugefield data\n",
    "gfData: List[np.ndarray] = []\n",
    "# Loop over each ID\n",
    "for iid in gfIDs:\n",
    "    gfFile = os.path.join(gfDir,f'{gfName}{iid}')\n",
    "    # Load and append\n",
    "    # Here we have specified the full path and the format\n",
    "    # Can figure it out based on extension, but format is clearer\n",
    "    gfData.append(lio.load(gfFile, format='openqcd'))\n",
    "# Convert to array for better indexing\n",
    "gfAr = np.asarray(gfData)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2c43ae3a",
   "metadata": {},
   "source": [
    "`lyncs.io` can also be used to convert from one format to another. For example, we can simply convert from `openqcd` format to `lime` using the `save()` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "4d854b8e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save the gaugefield array gfAr to a lime file\n",
    "lio.save(gfAr,'Gen2_8x24_gfAr.lime')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "99981925",
   "metadata": {},
   "source": [
    "Since our new file has a standard extension, `lyncs.io` can infer the format from the filename. The `head()` function now accesses the `lime` record associated with our new data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "74c58510",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'_lyncs_io': '0.2.3',\n",
       " 'created': '2023-11-02 14:37:29',\n",
       " 'type': \"<class 'numpy.ndarray'>\",\n",
       " 'shape': (3, 8, 24, 24, 24, 4, 3, 3),\n",
       " 'dtype': dtype('>c16'),\n",
       " 'fortran_order': False,\n",
       " 'descr': '<c16',\n",
       " 'nbytes': 191102976,\n",
       " '_offset': 856}"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Read the header of our new lime file\n",
    "lio.head('Gen2_8x24_gfAr.lime')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0cedf639",
   "metadata": {},
   "source": [
    "If we want to access all of the records in the file, we can use `lime.read_records()`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "d7b83558",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[{'_fp': <_io.BufferedReader name='Gen2_8x24_gfAr.lime'>,\n",
       "  'magic_number': 1164413355,\n",
       "  'version': 1,\n",
       "  'msg_bits': 32768,\n",
       "  'nbytes': 181,\n",
       "  'lime_type': 'xlf-info',\n",
       "  'offset': 144,\n",
       "  'begin': True,\n",
       "  'end': False,\n",
       "  'data': b\"_lyncs_io = 0.2.3\\ncreated = 2023-11-02 14:37:29\\ntype = <class 'numpy.ndarray'>\\nshape = (3, 8, 24, 24, 24, 4, 3, 3)\\ndtype = >c16\\nfortran_order = False\\ndescr = <c16\\nnbytes = 191102976\"},\n",
       " {'_fp': <_io.BufferedReader name='Gen2_8x24_gfAr.lime'>,\n",
       "  'magic_number': 1164413355,\n",
       "  'version': 1,\n",
       "  'msg_bits': 0,\n",
       "  'nbytes': 237,\n",
       "  'lime_type': 'lyncs-io-info',\n",
       "  'offset': 472,\n",
       "  'begin': False,\n",
       "  'end': False,\n",
       "  'data': b\"\\x80\\x04\\x95\\xe2\\x00\\x00\\x00\\x00\\x00\\x00\\x00}\\x94(\\x8c\\t_lyncs_io\\x94\\x8c\\x050.2.3\\x94\\x8c\\x07created\\x94\\x8c\\x132023-11-02 14:37:29\\x94\\x8c\\x04type\\x94\\x8c\\x17<class 'numpy.ndarray'>\\x94\\x8c\\x05shape\\x94(K\\x03K\\x08K\\x18K\\x18K\\x18K\\x04K\\x03K\\x03t\\x94\\x8c\\x05dtype\\x94\\x8c\\x05numpy\\x94\\x8c\\x05dtype\\x94\\x93\\x94\\x8c\\x03c16\\x94\\x89\\x88\\x87\\x94R\\x94(K\\x03\\x8c\\x01>\\x94NNNJ\\xff\\xff\\xff\\xffJ\\xff\\xff\\xff\\xffK\\x00t\\x94b\\x8c\\rfortran_order\\x94\\x89\\x8c\\x05descr\\x94\\x8c\\x04<c16\\x94\\x8c\\x06nbytes\\x94J\\x00\\x00d\\x0bu.\"},\n",
       " {'_fp': <_io.BufferedReader name='Gen2_8x24_gfAr.lime'>,\n",
       "  'magic_number': 1164413355,\n",
       "  'version': 1,\n",
       "  'msg_bits': 16384,\n",
       "  'nbytes': 191102976,\n",
       "  'lime_type': 'ildg-binary-data',\n",
       "  'offset': 856,\n",
       "  'begin': False,\n",
       "  'end': True,\n",
       "  'data': None}]"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lio.lime.read_records('Gen2_8x24_gfAr.lime')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f4da143a",
   "metadata": {},
   "source": [
    "If we want to read this lime file back into a `numpy.ndarray` we can use the `load()` function again"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "e03f5f74",
   "metadata": {},
   "outputs": [],
   "source": [
    "lime = lio.load('Gen2_8x24_gfAr.lime')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "edf8e158",
   "metadata": {},
   "source": [
    "## Manipulating the data\n",
    "\n",
    "As an exercise, let's calculate the values for the whole, spatial and temporal plaquettes:\n",
    "<img src=\"latticePlaqDiag.png\" alt=\"plaquette diagram\" width=\"400\" height=\"400\"/>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "1e707f42",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(3, 8, 24, 24, 24, 4, 3, 3)\n"
     ]
    }
   ],
   "source": [
    "# Here show the shape of the data\n",
    "# it is [ID, NT, NX, NY, NZ, mu, colour, colour]\n",
    "print(gfAr.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "013e5718",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "calculated average plaquette 1.626598501026478 in 2.4e+01 seconds\n"
     ]
    }
   ],
   "source": [
    "# Use Ryan's plaquette code to calculate whole of lattice plaquette\n",
    "# the sum of plaquettes, the number of plaquttes, the average plaquette and the time taken\n",
    "sumTrP, nP, ave, time = RPlaq.plaquette(gfAr[0, ...])\n",
    "print(f'calculated average plaquette {ave} in {time:.2} seconds')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "af563690",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "calculated average spatial plaquette 1.0088805585220895 in 1.2e+01 seconds\n"
     ]
    }
   ],
   "source": [
    "# Use Ryan's plaquette code to calculate spatial plaquette\n",
    "ssumTrP, snP, save, stime = RPlaq.plaquette(gfAr[0, ...], muStart=1, muEnd=4, nuEnd=4)\n",
    "print(f'calculated average spatial plaquette {save} in {stime:.2} seconds')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "34861e02",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "calculated average temporal plaquette 2.2443164435307277 in 1.2e+01 seconds\n"
     ]
    }
   ],
   "source": [
    "# Use Ryan's plaquette code to calculate temporal plaquette\n",
    "tsumTrP, tnP, tave, ttime = RPlaq.plaquette(gfAr[0, ...], muStart=0, muEnd=1, nuEnd=4)\n",
    "print(f'calculated average temporal plaquette {tave} in {ttime:.2} seconds')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "da9ef2b6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "average of temporal and spatial plaquettes 1.6265985010264086\n"
     ]
    }
   ],
   "source": [
    "# Just check that this agrees with the whole of lattice plaquette (visually)\n",
    "print(f'average of temporal and spatial plaquettes {(tave+save)/2.0}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "84c85b15",
   "metadata": {},
   "source": [
    "# C\n",
    "\n",
    "Now that we can load and manipulate the data, we might want to increase performance by offloading our calculations to a script written in a compiled language such as C."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "5fa360b2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Prepare a single gaugefield for writing\n",
    "COrder = gfAr[0, ...]\n",
    "# Now save\n",
    "COrder.tofile('Gen2_8x24_gfAr0.C')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "88153b6f",
   "metadata": {},
   "source": [
    "There is a C program `readC.c` supplied. This program will read `Gen2_8x24_gfAr0.C` in and calculate whole, spatial and temporal plaquettes. We do that now"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "f0c8edfc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# The command to compile\n",
    "ccmd = 'gcc -O3 readC.c'\n",
    "os.system(ccmd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "ff553fed",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Type:\t\tsumReTrP:\tnP:\tAvg:\n",
      "whole\t\t1079332.688553\t663552\t1.626599\n",
      "spatial\t\t334722.356184\t331776\t1.008881\n",
      "temporal\t744610.332369\t331776\t2.244316\n",
      "Total execution time: 0.108303s\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# The command to run\n",
    "rcmd = './a.out'\n",
    "os.system(rcmd)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "501132b5",
   "metadata": {},
   "source": [
    "# Fortran\n",
    "We might like to read the gaugefield into Fortran too. While there exists code that read openqcd format into Fortran, these do not undo the checkerboarding or `flattening` of the space-time dimensions. Here we seek to have the data in the same $N_t \\times N_s^3 \\times N_d \\times N_c^2$ format as in the `np.ndarray`.\n",
    "\n",
    "Unlike C which uses row-major ordering, Fortran uses column-major ordering for multidimensional arrays in linear memory. We must specify this reordering when calling `reshape` by passing the argument `order='F'`.\n",
    "\n",
    "<div class=\"alert alert-block alert-info\"><b>Info:</b> Python natively uses neither row nor column ordering. Instead, the allocations are made directly onto the heap (which is not necessarily contiguous) rather than the stack. However, the <code>numpy</code> package is based in C and thus follows the row-major ordering scheme.</div>\n",
    "\n",
    "This solution is based upon [this stack overflow](https://stackoverflow.com/a/49179272) answer. This solution assumes that your Python and Fortran code use the same (little) endianess."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "dacdd362",
   "metadata": {},
   "outputs": [],
   "source": [
    "# First reorder a single gaugefield into 'fortran' order\n",
    "fortOrder = gfAr[0, ...].reshape(gfAr[0, ...].shape, order='F')\n",
    "# we consider a single gaugefield only as Fortran allows only rank 7 arrays \n",
    "# and saving all configurations at once would be rank 8\n",
    "# This could be avoided by some sort of derived type\n",
    "# Now we save\n",
    "fortOrder.T.tofile('Gen2_8x24_gfAr0.fort')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b31c28f9",
   "metadata": {},
   "source": [
    "There is a Fortran program `readFortran.f90` supplied. This program will read `Gen2_8x24_gfAr0.fort` in and calculate whole, spatial and temporal plaquettes. We do that now"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "id": "273febd2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# The command to compile\n",
    "ccmd = 'gfortran -O3 -g -fbacktrace readFortran.f90'\n",
    "os.system(ccmd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "id": "12612168",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " type, sumTrp, nP, ave, time (seconds)\n",
      " whole   1079332.6885531216           663552   1.6265985010264781       0.14878800511360168     \n",
      " spatial   334722.35618422477           331776   1.0088805585220895        6.2325984239578247E-002\n",
      " temporal   744610.33236885071           331776   2.2443164435307277        6.2561005353927612E-002\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# The command to run\n",
    "rcmd = './a.out'\n",
    "os.system(rcmd)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "50d00afc",
   "metadata": {},
   "source": [
    "# F2PY\n",
    "Here we use the <code>F2PY</code> available through <code>NumPy</code> with documentation available [here](https://numpy.org/doc/stable/f2py/).\n",
    "\n",
    "We have written a small module in fortran, which has the essentially the same plaquette code as in `readFortran.f90`. This is in the file `fortPlaq.f90`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "id": "d9defbb3",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/ryan/bin/miniconda3/envs/lyncs/lib/python3.10/site-packages/numpy/f2py/f2py2e.py:686: VisibleDeprecationWarning: distutils has been deprecated since NumPy 1.26.Use the Meson backend instead, or generate wrapperswithout -c and use a custom build script\n",
      "  builder = build_backend(\n",
      "/home/ryan/bin/miniconda3/envs/lyncs/lib/python3.10/site-packages/setuptools/_distutils/cmd.py:66: SetuptoolsDeprecationWarning: setup.py install is deprecated.\n",
      "!!\n",
      "\n",
      "        ********************************************************************************\n",
      "        Please avoid running ``setup.py`` directly.\n",
      "        Instead, use pypa/build, pypa/installer or other\n",
      "        standards-based tools.\n",
      "\n",
      "        See https://blog.ganssle.io/articles/2021/10/setup-py-deprecated.html for details.\n",
      "        ********************************************************************************\n",
      "\n",
      "!!\n",
      "  self.initialize_options()\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "running build\n",
      "running config_cc\n",
      "INFO: unifing config_cc, config, build_clib, build_ext, build commands --compiler options\n",
      "running config_fc\n",
      "INFO: unifing config_fc, config, build_clib, build_ext, build commands --fcompiler options\n",
      "running build_src\n",
      "INFO: build_src\n",
      "INFO: building extension \"fortPlaq\" sources\n",
      "INFO: f2py options: []\n",
      "INFO: f2py:> /tmp/tmp0jz_s1i7/src.linux-x86_64-3.10/fortPlaqmodule.c\n",
      "creating /tmp/tmp0jz_s1i7/src.linux-x86_64-3.10\n",
      "Reading fortran codes...\n",
      "\tReading file 'fortPlaq.f90' (format:free)\n",
      "Post-processing...\n",
      "\tBlock: fortPlaq\n",
      "\t\t\tBlock: types\n",
      "In: :fortPlaq:fortPlaq.f90:types\n",
      "get_parameters: got \"eval() arg 1 must be a string, bytes or code object\" on 8\n",
      "In: :fortPlaq:fortPlaq.f90:types\n",
      "get_parameters: got \"eval() arg 1 must be a string, bytes or code object\" on 8\n",
      "\t\t\tBlock: plaq\n",
      "\t\t\t\tBlock: multiplymatmat\n",
      "In: :fortPlaq:fortPlaq.f90:plaq:multiplymatmat\n",
      "get_parameters: got \"eval() arg 1 must be a string, bytes or code object\" on 8\n",
      "In: :fortPlaq:fortPlaq.f90:plaq:multiplymatmat\n",
      "get_parameters: got \"eval() arg 1 must be a string, bytes or code object\" on 8\n",
      "\t\t\t\tBlock: multiplymatdagmatdag\n",
      "In: :fortPlaq:fortPlaq.f90:plaq:multiplymatdagmatdag\n",
      "get_parameters: got \"eval() arg 1 must be a string, bytes or code object\" on 8\n",
      "In: :fortPlaq:fortPlaq.f90:plaq:multiplymatdagmatdag\n",
      "get_parameters: got \"eval() arg 1 must be a string, bytes or code object\" on 8\n",
      "\t\t\t\tBlock: realtracemultmatmat\n",
      "In: :fortPlaq:fortPlaq.f90:plaq:realtracemultmatmat\n",
      "get_parameters: got \"eval() arg 1 must be a string, bytes or code object\" on 8\n",
      "In: :fortPlaq:fortPlaq.f90:plaq:realtracemultmatmat\n",
      "get_parameters: got \"eval() arg 1 must be a string, bytes or code object\" on 8\n",
      "\t\t\t\tBlock: plaquette\n",
      "In: :fortPlaq:fortPlaq.f90:plaq:plaquette\n",
      "get_parameters: got \"eval() arg 1 must be a string, bytes or code object\" on 8\n",
      "In: :fortPlaq:fortPlaq.f90:plaq:plaquette\n",
      "get_parameters: got \"eval() arg 1 must be a string, bytes or code object\" on 8\n",
      "Applying post-processing hooks...\n",
      "  character_backward_compatibility_hook\n",
      "Post-processing (stage 2)...\n",
      "\tBlock: fortPlaq\n",
      "\t\tBlock: unknown_interface\n",
      "\t\t\tBlock: types\n",
      "\t\t\tBlock: plaq\n",
      "\t\t\t\tBlock: multiplymatmat\n",
      "\t\t\t\tBlock: multiplymatdagmatdag\n",
      "\t\t\t\tBlock: realtracemultmatmat\n",
      "\t\t\t\tBlock: plaquette\n",
      "Building modules...\n",
      "    Building module \"fortPlaq\"...\n",
      "\t\tConstructing F90 module support for \"types\"...\n",
      "\t\t  Variables: dc dp r_dc r_dp\n",
      "\t\tConstructing F90 module support for \"plaq\"...\n",
      "            Constructing wrapper function \"plaq.multiplymatmat\"...\n",
      "              mm = multiplymatmat(left,right)\n",
      "            Constructing wrapper function \"plaq.multiplymatdagmatdag\"...\n",
      "              mm = multiplymatdagmatdag(left,right)\n",
      "            Constructing wrapper function \"plaq.realtracemultmatmat\"...\n",
      "              trmm = realtracemultmatmat(left,right)\n",
      "\t\tCreating wrapper for Fortran subroutine \"plaquette\"(\"plaquette\")...\n",
      "            Constructing wrapper function \"plaq.plaquette\"...\n",
      "              sumtrp,np,time = plaquette(data,mustart,muend,nuend)\n",
      "    Wrote C/API module \"fortPlaq\" to file \"/tmp/tmp0jz_s1i7/src.linux-x86_64-3.10/fortPlaqmodule.c\"\n",
      "    Fortran 90 wrappers are saved to \"/tmp/tmp0jz_s1i7/src.linux-x86_64-3.10/fortPlaq-f2pywrappers2.f90\"\n",
      "INFO:   adding '/tmp/tmp0jz_s1i7/src.linux-x86_64-3.10/fortranobject.c' to sources.\n",
      "INFO:   adding '/tmp/tmp0jz_s1i7/src.linux-x86_64-3.10' to include_dirs.\n",
      "copying /home/ryan/bin/miniconda3/envs/lyncs/lib/python3.10/site-packages/numpy/f2py/src/fortranobject.c -> /tmp/tmp0jz_s1i7/src.linux-x86_64-3.10\n",
      "copying /home/ryan/bin/miniconda3/envs/lyncs/lib/python3.10/site-packages/numpy/f2py/src/fortranobject.h -> /tmp/tmp0jz_s1i7/src.linux-x86_64-3.10\n",
      "INFO:   adding '/tmp/tmp0jz_s1i7/src.linux-x86_64-3.10/fortPlaq-f2pywrappers2.f90' to sources.\n",
      "INFO: build_src: building npy-pkg config files\n",
      "running build_ext\n",
      "INFO: customize UnixCCompiler\n",
      "INFO: customize UnixCCompiler using build_ext\n",
      "INFO: get_default_fcompiler: matching types: '['arm', 'gnu95', 'intel', 'lahey', 'pg', 'nv', 'absoft', 'nag', 'vast', 'compaq', 'intele', 'intelem', 'gnu', 'g95', 'pathf95', 'nagfor', 'fujitsu']'\n",
      "INFO: customize ArmFlangCompiler\n",
      "WARN: Could not locate executable armflang\n",
      "INFO: customize Gnu95FCompiler\n",
      "INFO: Found executable /home/ryan/bin/miniconda3/envs/lyncs/bin/gfortran\n",
      "INFO: customize Gnu95FCompiler\n",
      "INFO: customize Gnu95FCompiler using build_ext\n",
      "INFO: building 'fortPlaq' extension\n",
      "INFO: compiling C sources\n",
      "INFO: C compiler: gcc -pthread -B /home/ryan/bin/miniconda3/envs/lyncs/compiler_compat -Wno-unused-result -Wsign-compare -DNDEBUG -fwrapv -O2 -Wall -fPIC -O2 -isystem /home/ryan/bin/miniconda3/envs/lyncs/include -fPIC -O2 -isystem /home/ryan/bin/miniconda3/envs/lyncs/include -fPIC\n",
      "\n",
      "creating /tmp/tmp0jz_s1i7/tmp\n",
      "creating /tmp/tmp0jz_s1i7/tmp/tmp0jz_s1i7\n",
      "creating /tmp/tmp0jz_s1i7/tmp/tmp0jz_s1i7/src.linux-x86_64-3.10\n",
      "INFO: compile options: '-DNPY_DISABLE_OPTIMIZATION=1 -I/tmp/tmp0jz_s1i7/src.linux-x86_64-3.10 -I/home/ryan/bin/miniconda3/envs/lyncs/lib/python3.10/site-packages/numpy/core/include -I/home/ryan/bin/miniconda3/envs/lyncs/include/python3.10 -c'\n",
      "INFO: gcc: /tmp/tmp0jz_s1i7/src.linux-x86_64-3.10/fortPlaqmodule.c\n",
      "INFO: gcc: /tmp/tmp0jz_s1i7/src.linux-x86_64-3.10/fortranobject.c\n",
      "INFO: compiling Fortran 90 module sources\n",
      "INFO: Fortran f77 compiler: /home/ryan/bin/miniconda3/envs/lyncs/bin/gfortran -Wall -g -ffixed-form -fno-second-underscore -fPIC -O3 -funroll-loops\n",
      "Fortran f90 compiler: /home/ryan/bin/miniconda3/envs/lyncs/bin/gfortran -Wall -g -fno-second-underscore -fPIC -O3 -funroll-loops\n",
      "Fortran fix compiler: /home/ryan/bin/miniconda3/envs/lyncs/bin/gfortran -Wall -g -ffixed-form -fno-second-underscore -Wall -g -fno-second-underscore -fPIC -O3 -funroll-loops\n",
      "INFO: compile options: '-I/tmp/tmp0jz_s1i7/src.linux-x86_64-3.10 -I/home/ryan/bin/miniconda3/envs/lyncs/lib/python3.10/site-packages/numpy/core/include -I/home/ryan/bin/miniconda3/envs/lyncs/include/python3.10 -c'\n",
      "extra options: '-J/tmp/tmp0jz_s1i7/ -I/tmp/tmp0jz_s1i7/'\n",
      "INFO: gfortran:f90: fortPlaq.f90\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "fortPlaq.f90:12:26:\n",
      "\n",
      "   12 |         real(dc) :: r_dc = (1.0_dp, 1.0_dp)\n",
      "      |                          1\n",
      "Warning: Non-zero imaginary part discarded in conversion from 'COMPLEX(8)' to 'REAL(8)' at (1) [-Wconversion]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "INFO: compiling Fortran sources\n",
      "INFO: Fortran f77 compiler: /home/ryan/bin/miniconda3/envs/lyncs/bin/gfortran -Wall -g -ffixed-form -fno-second-underscore -fPIC -O3 -funroll-loops\n",
      "Fortran f90 compiler: /home/ryan/bin/miniconda3/envs/lyncs/bin/gfortran -Wall -g -fno-second-underscore -fPIC -O3 -funroll-loops\n",
      "Fortran fix compiler: /home/ryan/bin/miniconda3/envs/lyncs/bin/gfortran -Wall -g -ffixed-form -fno-second-underscore -Wall -g -fno-second-underscore -fPIC -O3 -funroll-loops\n",
      "INFO: compile options: '-I/tmp/tmp0jz_s1i7/src.linux-x86_64-3.10 -I/home/ryan/bin/miniconda3/envs/lyncs/lib/python3.10/site-packages/numpy/core/include -I/home/ryan/bin/miniconda3/envs/lyncs/include/python3.10 -c'\n",
      "extra options: '-J/tmp/tmp0jz_s1i7/ -I/tmp/tmp0jz_s1i7/'\n",
      "INFO: gfortran:f90: /tmp/tmp0jz_s1i7/src.linux-x86_64-3.10/fortPlaq-f2pywrappers2.f90\n",
      "INFO: /home/ryan/bin/miniconda3/envs/lyncs/bin/gfortran -Wall -g -Wall -g -shared /tmp/tmp0jz_s1i7/tmp/tmp0jz_s1i7/src.linux-x86_64-3.10/fortPlaqmodule.o /tmp/tmp0jz_s1i7/tmp/tmp0jz_s1i7/src.linux-x86_64-3.10/fortranobject.o /tmp/tmp0jz_s1i7/fortPlaq.o /tmp/tmp0jz_s1i7/tmp/tmp0jz_s1i7/src.linux-x86_64-3.10/fortPlaq-f2pywrappers2.o -L/home/ryan/bin/miniconda3/envs/lyncs/bin/../lib/gcc/x86_64-conda-linux-gnu/13.2.0/../../../../x86_64-conda-linux-gnu/lib/../lib -L/home/ryan/bin/miniconda3/envs/lyncs/bin/../lib/gcc/x86_64-conda-linux-gnu/13.2.0/../../../../x86_64-conda-linux-gnu/lib/../lib -lgfortran -o ./fortPlaq.cpython-310-x86_64-linux-gnu.so\n",
      "Removing build directory /tmp/tmp0jz_s1i7\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# First we have to 'compile' the module. Luckily we already have gfortran from earlier\n",
    "ccmd = 'f2py -c fortPlaq.f90 -m fortPlaq'\n",
    "os.system(ccmd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "id": "f006072f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "calculated average plaquette 1.626598501026478 in 0.11 seconds\n",
      "calculated average spatial plaquette 1.0088805585220895 in 0.055 seconds\n",
      "calculated average temporal plaquette 2.2443164435307277 in 0.059 seconds\n"
     ]
    }
   ],
   "source": [
    "# Now import the module\n",
    "import fortPlaq                               \n",
    "fPlaq = fortPlaq.plaq.plaquette\n",
    "\n",
    "# separately calculate plaq\n",
    "sP, nP, time = fPlaq(gfAr[0, ...], mustart=1, muend=4, nuend=4)\n",
    "print(f'calculated average plaquette {sP / nP} in {time:.2} seconds')\n",
    "# calculate splaq\n",
    "sP, nP, time = fPlaq(gfAr[0, ...], mustart=2, muend=4, nuend=4)\n",
    "print(f'calculated average spatial plaquette {sP / nP} in {time:.2} seconds')\n",
    "# calculate tplaq\n",
    "sP, nP, time = fPlaq(gfAr[0, ...], mustart=1, muend=1, nuend=4)\n",
    "print(f'calculated average temporal plaquette {sP / nP} in {time:.2} seconds')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7de37e21",
   "metadata": {},
   "source": [
    "### "
   ]
  }
 ],
 "metadata": {
  "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.10.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}