RCS_nearest-neighbor_MARP_search.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "5f57c8a3",
   "metadata": {},
   "source": [
    "# RCS (Nearest-Neighbor) Benchmark"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "48d34546-e54a-442f-85c5-e9a28d7c38d1",
   "metadata": {},
   "source": [
    "## Getting Started"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7b171b3b-012c-4dad-870a-e8be44c24128",
   "metadata": {},
   "source": [
    "You need [`pyqrack`](https://pypi.org/project/pyqrack/) to run this notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "b921e879-46df-4382-a1c7-804b6f2d0854",
   "metadata": {},
   "outputs": [],
   "source": [
    "# For example, if your Jupyter installation uses pip:\n",
    "# import sys\n",
    "# !{sys.executable} -m pip install pyqrack"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1e55837c-018f-4cb2-a9cb-86f356d7a2d8",
   "metadata": {},
   "source": [
    "In the Python package itself, there should be an executable called `qrack_cl_precompile`. This \"pre-compiles\" the OpenCL \"just-in-time\" (\"JIT\") device program, for your system's accelerators. You might want to find and run this utility first, to avoid the need to \"recompile\" the OpenCL device program every time you first load Qrack into your environment."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f0b34caa-b73e-4059-a27b-ad3e828a623d",
   "metadata": {},
   "source": [
    "## Configuration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "90b257ae",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Random trial qubit width\n",
    "qubit_width = 36\n",
    "# Random trial circuit layer depth\n",
    "layer_depth = 6\n",
    "# Increment for SDRP\n",
    "sdrp_inc = 0.125"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "1f21a89b",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2024-02-21 12:22:28.984079: W tensorflow/compiler/tf2tensorrt/utils/py_utils.cc:38] TF-TRT Warning: Could not find TensorRT\n"
     ]
    }
   ],
   "source": [
    "import math\n",
    "import numpy as np\n",
    "import os\n",
    "import random\n",
    "import time\n",
    "\n",
    "from pyqrack import QrackSimulator, QrackCircuit\n",
    "\n",
    "# For more specific details about all available Qrack environment variables,\n",
    "# See the C++ repository README: https://github.com/unitaryfund/qrack\n",
    "\n",
    "# This is the maximum total number of fully-entangled qubits you expect to fit in general RAM.\n",
    "os.environ['QRACK_MAX_CPU_QB']='32'\n",
    "\n",
    "# Above this threshold, \"QTensorNetwork\" restricts simulations to \"past light cone.\"\n",
    "# At or below the threshold, much more work can be reused.\n",
    "os.environ['QRACK_QTENSORNETWORK_THRESHOLD_QB']='-1'\n",
    "\n",
    "# These below are approximation options. (By default, Qrack simulates in the \"ideal.\")\n",
    "\n",
    "# This is a number between \"0\" (\"ideal\") and \"1\" (\"round to exactly Clifford\") for near-Clifford rounding.\n",
    "# os.environ['QRACK_NONCLIFFORD_ROUNDING_THRESHOLD']='1\n",
    "\n",
    "# This is a number between \"0\" (\"ideal\") and \"1\" (\"destroy all entanglement\") for \"SDRP,\"\n",
    "# \"Schmidt decomposition rounding parameter\". (https://arxiv.org/abs/2304.14969)\n",
    "# os.environ['QRACK_QUNIT_SEPARABILITY_THRESHOLD']='0.6'\n",
    "\n",
    "# This is a number between \"0\" (\"ideal\") and \"1\" (\"combine all binary decision tree branches\")\n",
    "# that sets the allowable \"epsilon\" between \"QBdt\" branches to consider them equal.\n",
    "# os.environ['QRACK_QBDT_SEPARABILITY_THRESHOLD']='0.0001'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fb4ef620",
   "metadata": {},
   "source": [
    "## Run the Benchmark"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "9e742f47",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "36 qubits by 6 circuit layers:\n",
      "Device #0, Loaded binary from: /home/iamu/.qrack/qrack_ocl_dev_NVIDIA_GeForce_RTX_2070_Super.ir\n",
      "0.128010116997757  seconds, 1.252129010262608e-10 out of 1.0 fidelity, SDRP level 1\n",
      "Default platform: NVIDIA CUDA\n",
      "Default device: #0, NVIDIA GeForce RTX 2070 Super\n",
      "OpenCL device #0: NVIDIA GeForce RTX 2070 Super\n",
      "Failed to enqueue kernel, error code: -4\n",
      "(MARP search took  0.4804771119997895  seconds overall.)\n"
     ]
    }
   ],
   "source": [
    "_sqrt1_2 = 1 / np.sqrt(2)\n",
    "_sqrtI = np.sqrt(1j)\n",
    "\n",
    "def u(circ, q, th, ph, lm):\n",
    "    circ.mtrx([np.cos(th / 2), -np.exp(1j * lm) * np.sin(th / 2), np.exp(1j * ph) * np.sin(th / 2), np.exp(1j * (ph + lm)) * np.cos(th / 2)], q)\n",
    "\n",
    "def h(circ, q):\n",
    "    circ.mtrx([_sqrt1_2, _sqrt1_2, _sqrt1_2, -_sqrt1_2], q)\n",
    "\n",
    "def x(circ, q):\n",
    "    circ.mtrx([0, 1, 1, 0], q)\n",
    "\n",
    "def sqrtx(circ, q):\n",
    "    circ.mtrx([(1 + 1j) / 2, (1 - 1j) / 2, (1 - 1j) / 2, (1 + 1j) / 2], q)\n",
    "\n",
    "def y(circ, q):\n",
    "    circ.mtrx([0, -1j, 1j, 0], q)\n",
    "\n",
    "def sqrty(circ, q):\n",
    "    circ.mtrx([(1 + 1j) / 2, -(1 + 1j) / 2, (1 + 1j) / 2, (1 + 1j) / 2], q)\n",
    "\n",
    "def z(circ, q):\n",
    "    circ.mtrx([1, 0, 0, -1], q)\n",
    "\n",
    "def s(circ, q):\n",
    "    circ.mtrx([1, 0, 0, 1j], q)\n",
    "\n",
    "def t(circ, q):\n",
    "    circ.mtrx([1, 0, 0, _sqrtI], q)\n",
    "\n",
    "def rz(circ, q):\n",
    "    circ.mtrx([1, 0, 0, np.exp(1j * random.uniform(0, 2 * np.pi))], q)\n",
    "\n",
    "def cx(circ, q1, q2):\n",
    "    circ.ucmtrx([q1], [0, 1, 1, 0], q2, 1)\n",
    "\n",
    "def cy(circ, q1, q2):\n",
    "    circ.ucmtrx([q1], [0, -1j, 1j, 0], q2, 1)\n",
    "\n",
    "def cz(circ, q1, q2):\n",
    "    circ.ucmtrx([q1], [1, 0, 0, -1], q2, 1)\n",
    "\n",
    "def acx(circ, q1, q2):\n",
    "    circ.ucmtrx([q1], [0, 1, 1, 0], q2, 0)\n",
    "\n",
    "def acy(circ, q1, q2):\n",
    "    circ.ucmtrx([q1], [0, -1j, 1j, 0], q2, 0)\n",
    "\n",
    "def acz(circ, q1, q2):\n",
    "    circ.ucmtrx([q1], [1, 0, 0, -1], q2, 0)\n",
    "\n",
    "def swap(circ, q1, q2):\n",
    "    circ.swap(q1, q2)\n",
    "\n",
    "def nswap(circ, q1, q2):\n",
    "    circ.ucmtrx([q1], [1, 0, 0, -1], q2, 0)\n",
    "    circ.swap(q1, q2)\n",
    "    circ.ucmtrx([q1], [1, 0, 0, -1], q2, 0)\n",
    "\n",
    "def pswap(circ, q1, q2):\n",
    "    circ.ucmtrx([q1], [1, 0, 0, -1], q2, 0)\n",
    "    circ.swap(q1, q2)\n",
    "\n",
    "def mswap(circ, q1, q2):\n",
    "    circ.swap(q1, q2)\n",
    "    circ.ucmtrx([q1], [1, 0, 0, -1], q2, 0)\n",
    "\n",
    "def iswap(circ, q1, q2):\n",
    "    circ.swap(q1, q2)\n",
    "    circ.ucmtrx([q1], [1, 0, 0, -1], q2, 1)\n",
    "    circ.mtrx([1, 0, 0, 1j], q1)\n",
    "    circ.mtrx([1, 0, 0, 1j], q2)\n",
    "\n",
    "def iiswap(circ, q1, q2):\n",
    "    circ.mtrx([1, 0, 0, -1j], q2)\n",
    "    circ.mtrx([1, 0, 0, -1j], q1)\n",
    "    circ.ucmtrx([q1], [1, 0, 0, -1], q2, 1)\n",
    "    circ.swap(q1, q2)\n",
    "\n",
    "# This is a circuit with nearest-neighbor couplers that aims for high Haar-randomness.\n",
    "def generate_random_circuit(n, l):\n",
    "    n_range = range(n)\n",
    "    l_range = range(l)\n",
    "\n",
    "    # This is the coupler order deemed \"hard\" for Sycamore simulation in 2019:\n",
    "    gateSequence = [ 0, 3, 2, 1, 2, 1, 0, 3 ]\n",
    "\n",
    "    # Qubit layout is perfect square greater-than-or-equal to width,\n",
    "    # with one partial remainder row if width is not a perfect square.\n",
    "    row_len = math.ceil(math.sqrt(n))\n",
    "\n",
    "    # Half (phased) swap, half singly-controlled\n",
    "    two_bit_gates = swap, pswap, mswap, nswap, iswap, iiswap, cx, cy, cz, acx, acy, acz\n",
    "\n",
    "    # Discrete:\n",
    "    # one_bit_gates = h, (x, sqrtx), (y, sqrty), (z, s, t)\n",
    "\n",
    "    # Continuous (rz):\n",
    "    one_bit_gates = h, (x, sqrtx), (y, sqrty), rz\n",
    "\n",
    "    circ = QrackCircuit(False)\n",
    "\n",
    "    for _ in l_range:\n",
    "        # Single-qubit gates\n",
    "        for i in n_range:\n",
    "            # Random gate choice from a discrete set:\n",
    "            g = random.choice(one_bit_gates)\n",
    "            if type(g) is tuple:\n",
    "                # The \"gate choice\" segments a Pauli axis, in binary choices\n",
    "                # (We could just pull a 50/50 bit for each segment gate,\n",
    "                # but, here, we're trying to avoid no-op in all cases.)\n",
    "                choice = random.randint(1, (1 << len(g)) - 1)\n",
    "                for j in range(len(g)):\n",
    "                    if (choice >> j) & 1:\n",
    "                        g[j](circ, i)\n",
    "            else:\n",
    "                g(circ, i)\n",
    "            \n",
    "            # One Euler angle axis per gate:\n",
    "            # h(circ, i)\n",
    "            # rz(circ, i)\n",
    "\n",
    "            # Three Euler angle axes per gate:\n",
    "            # u(circ, i, random.uniform(0, 2 * math.pi), random.uniform(0, 2 * math.pi), random.uniform(0, 2 * math.pi))\n",
    "\n",
    "        # 2-qubit couplers\n",
    "        gate = gateSequence.pop(0)\n",
    "        gateSequence.append(gate)\n",
    "        for row in range(1, row_len, 2):\n",
    "            for col in range(row_len):\n",
    "                temp_row = row\n",
    "                temp_col = col\n",
    "                temp_row = temp_row + (1 if (gate & 2) else -1);\n",
    "                temp_col = temp_col + (1 if (gate & 1) else 0)\n",
    "\n",
    "                if (temp_row < 0) or (temp_col < 0) or (temp_row >= row_len) or (temp_col >= row_len):\n",
    "                    # Row and/or column selected were out of range\n",
    "                    continue\n",
    "\n",
    "                b1 = row * row_len + col\n",
    "                b2 = temp_row * row_len + temp_col\n",
    "\n",
    "                if (b1 >= n) or (b2 >= n):\n",
    "                    # Bits selected were out-of-range\n",
    "                    continue\n",
    "\n",
    "                # Swap bits, 50% of the time\n",
    "                if random.randint(0, 1) == 1:\n",
    "                    temp = b1\n",
    "                    b1 = b2\n",
    "                    b2 = temp\n",
    "                \n",
    "                g = random.choice(two_bit_gates)\n",
    "                g(circ, b1, b2)\n",
    "\n",
    "    # You can see the circuit definition and load it from this file name:\n",
    "    circ.out_to_file('marp_circuit.qc')\n",
    "\n",
    "    return circ\n",
    "\n",
    "def bench_qrack(n, l):\n",
    "    circ = generate_random_circuit(n, l)\n",
    "    \n",
    "    sdrp = 1\n",
    "    found_marp = False\n",
    "    result = (0, 0, 1)\n",
    "\n",
    "    while not found_marp:\n",
    "        # This is a number between \"0\" (\"ideal\") and \"1\" (\"destroy all entanglement\") for \"SDRP,\"\n",
    "        # \"Schmidt decomposition rounding parameter\". (https://arxiv.org/abs/2304.14969)\n",
    "        if (sdrp > 1e-6):\n",
    "            os.environ['QRACK_QUNIT_SEPARABILITY_THRESHOLD']=str(sdrp)\n",
    "        else:\n",
    "            del os.environ['QRACK_QUNIT_SEPARABILITY_THRESHOLD']\n",
    "            sdrp = 0\n",
    "            found_marp = True\n",
    "\n",
    "        start = time.perf_counter()\n",
    "        # SDRP is fixed at the point the simulator instance is created.\n",
    "        sim = QrackSimulator(n)\n",
    "        try:\n",
    "            circ.run(sim)\n",
    "            fidelity = sim.get_unitary_fidelity()\n",
    "            sim.m_all()\n",
    "            result = (time.perf_counter() - start, fidelity, sdrp)\n",
    "            print(result[0], \" seconds,\", result[1], \"out of 1.0 fidelity, SDRP level\", result[2])\n",
    "        except:\n",
    "            found_marp = True\n",
    "\n",
    "        sdrp = sdrp - sdrp_inc\n",
    "\n",
    "    return result\n",
    "\n",
    "print(qubit_width, \"qubits by\", layer_depth, \"circuit layers:\")\n",
    "start = time.perf_counter()\n",
    "result = bench_qrack(qubit_width, layer_depth)\n",
    "marp_time = time.perf_counter() - start\n",
    "\n",
    "print(\"(MARP search took \", marp_time, \" seconds overall.)\")"
   ]
  }
 ],
 "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}