Clifford_RZ_Mitiq_RCS.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "composed-midwest",
   "metadata": {},
   "source": [
    "# Qrack Clifford+RZ and Mitiq for High-Width RCS"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "social-scott",
   "metadata": {},
   "source": [
    "You need the `pyqrack` package to run this notebook. [`vm6502q/pyqrack`](https://github.com/vm6502q/pyqrack) is a pure Python wrapper on the [`vm6502q/qrack`](https://github.com/vm6502q/qrack) quantum computer simulation framework core library. The preferred method of installation is from source code, at those GitHub repositories, but a package with default build precompiled binaries is available on [pypi](https://pypi.org/project/pyqrack/0.2.0/)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8a647098-dd6d-4ee0-806f-abc2f5cd37ab",
   "metadata": {},
   "source": [
    "You also need the `mitiq` package. `unitaryfund/mitiq` is a Python package for error mitigation. It relies intermediate representations of circuits like Qiskit's `QuantumCircuit` or Cirq's `Circuit`. In the example below, we use Qiskit's `QuantumCircuit`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "e342bcef",
   "metadata": {},
   "outputs": [],
   "source": [
    "width = 54\n",
    "depth = 20\n",
    "max_magic = 10\n",
    "trials = 64"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "cubic-thousand",
   "metadata": {},
   "outputs": [],
   "source": [
    "# For example, if your Jupyter installation uses pip:\n",
    "# import sys\n",
    "# !{sys.executable} -m pip install pyqrack\n",
    "\n",
    "import os\n",
    "import collections\n",
    "import math\n",
    "import multiprocessing\n",
    "import numpy as np\n",
    "import random\n",
    "\n",
    "from pyqrack import QrackSimulator, Pauli\n",
    "from mitiq import zne\n",
    "from qiskit.circuit.quantumcircuit import QuantumCircuit\n",
    "from qiskit.compiler import transpile\n",
    "from mitiq.zne.scaling.folding import fold_global\n",
    "from mitiq.zne.inference import RichardsonFactory\n",
    "\n",
    "# Change these according to your OpenCL devices and system specifications:\n",
    "os.environ['QRACK_MAX_PAGE_QB']='27'\n",
    "os.environ['QRACK_MAX_ALLOC_MB']='15872'\n",
    "os.environ['QRACK_MAX_PAGING_QB']='30'\n",
    "os.environ['QRACK_MAX_CPU_QB']='-1'\n",
    "# os.environ['QRACK_QUNITMULTI_DEVICES']='1'\n",
    "# os.environ['QRACK_QPAGER_DEVICES']='1'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "59c4e106-4dcb-4c60-938f-a93da5f4a487",
   "metadata": {},
   "outputs": [],
   "source": [
    "def cx(circ, q1, q2):\n",
    "    circ.cx(q1, q2)\n",
    "\n",
    "def cy(circ, q1, q2):\n",
    "    circ.cy(q1, q2)\n",
    "\n",
    "def cz(circ, q1, q2):\n",
    "    circ.cz(q1, q2)\n",
    "\n",
    "def acx(circ, q1, q2):\n",
    "    circ.x(q1)\n",
    "    circ.cx(q1, q2)\n",
    "    circ.x(q1)\n",
    "\n",
    "def acy(circ, q1, q2):\n",
    "    circ.x(q1)\n",
    "    circ.cy(q1, q2)\n",
    "\n",
    "def acz(circ, q1, q2):\n",
    "    circ.x(q1)\n",
    "    circ.cz(q1, q2)\n",
    "    circ.x(q1)\n",
    "\n",
    "def swap(circ, q1, q2):\n",
    "    circ.swap(q1, q2)\n",
    "\n",
    "def nswap(circ, q1, q2):\n",
    "    circ.cz(q1, q2)\n",
    "    circ.swap(q1, q2)\n",
    "    circ.cz(q1, q2)\n",
    "\n",
    "def pswap(circ, q1, q2):\n",
    "    circ.cz(q1, q2)\n",
    "    circ.swap(q1, q2)   \n",
    "\n",
    "def mswap(circ, q1, q2):\n",
    "    circ.swap(q1, q2)\n",
    "    circ.cz(q1, q2)\n",
    "\n",
    "def iswap(circ, q1, q2):\n",
    "    circ.swap(q1, q2)\n",
    "    circ.cz(q1, q2)\n",
    "    circ.s(q1)\n",
    "    circ.s(q2)\n",
    "\n",
    "def iiswap(circ, q1, q2):\n",
    "    circ.sdg(q2)\n",
    "    circ.sdg(q1)\n",
    "    circ.cz(q1, q2)\n",
    "    circ.swap(q1, q2)\n",
    "\n",
    "def random_circuit(width, depth, magic):\n",
    "    circuit = QuantumCircuit(width)\n",
    "\n",
    "    magic_fraction = (3 * (width - 1) * depth / max_magic) if (width == 54) else (3 * width * depth / max_magic)\n",
    "    \n",
    "    # Nearest-neighbor couplers:\n",
    "    gateSequence = [ 0, 3, 2, 1, 2, 1, 0, 3 ]\n",
    "    two_bit_gates = swap, pswap, mswap, nswap, iswap, iiswap, cx, cy, cz, acx, acy, acz\n",
    "\n",
    "    col_len = math.floor(math.sqrt(width))\n",
    "    while (((width // col_len) * col_len) != width):\n",
    "        col_len -= 1\n",
    "    row_len = width // col_len\n",
    "    dead_qubit = 3\n",
    "\n",
    "    for i in range(depth):\n",
    "        # Single bit gates\n",
    "        for j in range(width):\n",
    "            if (width == 54) and (j == dead_qubit):\n",
    "                continue\n",
    "\n",
    "            for _ in range(3):\n",
    "                # We're trying to cover 3 Pauli axes\n",
    "                # with Euler angle axes x-z-x. \n",
    "                circuit.h(j)\n",
    "\n",
    "                # We can trace out a quarter rotations around the Bloch sphere with stabilizer.\n",
    "                rnd = random.randint(0, 3)\n",
    "                if rnd & 1:\n",
    "                    circuit.s(j)\n",
    "                if rnd & 2:\n",
    "                    circuit.z(j)\n",
    "\n",
    "                # For each axis, there is a chance of \"magic.\"\n",
    "                if (magic > 0) and ((magic_fraction * random.random()) < 1):\n",
    "                    angle = random.uniform(0, math.pi / 2)\n",
    "                    circuit.rz(angle, j)\n",
    "                    magic -= 1\n",
    "\n",
    "        # Nearest-neighbor couplers:\n",
    "        ############################\n",
    "        gate = gateSequence.pop(0)\n",
    "        gateSequence.append(gate)\n",
    "        for row in range(1, row_len, 2):\n",
    "            for col in range(col_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",
    "                    continue\n",
    "\n",
    "                b1 = row * row_len + col\n",
    "                b2 = temp_row * row_len + temp_col\n",
    "\n",
    "                if (b1 >= width) or (b2 >= width):\n",
    "                    continue\n",
    "\n",
    "                g = random.choice(two_bit_gates)\n",
    "                g(circuit, b1, b2)\n",
    "    \n",
    "    return circuit"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8c935d73-868e-4fc1-b75d-aff5f90ca566",
   "metadata": {},
   "source": [
    "Run this next cell several times to produce a noisy value, (optionally tuning the width, depth, and magic settings above,) and run the last cell to mitigate with Mitiq."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "ab200420-d7da-4977-89b1-021e707d4faf",
   "metadata": {},
   "outputs": [],
   "source": [
    "def logit(x):\n",
    "    # Theoretically, these limit points are \"infinite,\"\n",
    "    # but precision caps out between 36 and 37:\n",
    "    if x > (1 - 5e-17):\n",
    "        return 37\n",
    "    # For the negative limit, the precisions caps out\n",
    "    # between -37 and -38\n",
    "    elif x < 1e-17:\n",
    "        return -38\n",
    "    return max(-38, min(37, np.log(x / (1 - x))))\n",
    "\n",
    "def expit(x):\n",
    "    # Theoretically, these limit points are \"infinite,\"\n",
    "    # but precision caps out between 36 and 37:\n",
    "    if x >= 37:\n",
    "        return 1.0\n",
    "    # For the negative limit, the precisions caps out\n",
    "    # between -37 and -38\n",
    "    elif x <= -38:\n",
    "        return 0.0\n",
    "    return 1 / (1 + np.exp(-x))\n",
    "\n",
    "def execute(circuit):\n",
    "    \"\"\"Returns the mirror circuit expectation value for unsigned integer overall bit string.\"\"\"\n",
    "\n",
    "    circ = circuit.compose(circuit.inverse())\n",
    "    qsim = QrackSimulator(circ.width(), isSchmidtDecompose=False, isOpenCL=False)\n",
    "    qsim.run_qiskit_circuit(circ, 0)\n",
    "\n",
    "    qubits = list(range(width))\n",
    "    perm = [False] * width\n",
    "\n",
    "    # This is fidelity of the mirror circuit:\n",
    "    p = (qsim.prob_perm_rdm(qubits, perm, r=True) + qsim.prob_perm_rdm(qubits, perm, r=False)) / 2\n",
    "\n",
    "    # So as not to exceed floor at 0.0 and ceiling at 1.0, (assuming 0 < p < 1,)\n",
    "    # we mitigate its logit function value (https://en.wikipedia.org/wiki/Logit)\n",
    "    return logit(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "b0504807-3e1b-4e43-ba10-9f6ac4e14026",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fidelity with ZNE:  0.2555995354052111\n",
      "Width/depth/magic (before circuit mirror):  53 / 20 / 10\n",
      "CPU times: user 327 ms, sys: 117 ms, total: 444 ms\n",
      "Wall time: 18min 45s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "scale_count = 10\n",
    "max_scale = 2\n",
    "factory = RichardsonFactory(scale_factors=[(1 + (max_scale - 1) * x / scale_count) for x in range(0, scale_count)])\n",
    "\n",
    "def trial(procnum):\n",
    "    circ = random_circuit(width, depth, max_magic)\n",
    "    qsim = QrackSimulator(width, isSchmidtDecompose=False, isOpenCL=False)\n",
    "\n",
    "    return zne.execute_with_zne(circ, execute, scale_noise=fold_global, factory=factory)\n",
    "\n",
    "pool = multiprocessing.Pool(processes = multiprocessing.cpu_count())\n",
    "result = pool.map(trial, range(trials))\n",
    "pool.close()\n",
    "mitigated_fidelity = sum([expit(x) for x in result]) / trials\n",
    "\n",
    "print(\"Fidelity with ZNE: \", mitigated_fidelity)\n",
    "print(\"Width/depth/magic (before circuit mirror): \", 53 if width == 54 else width, \"/\", depth, \"/\", max_magic)"
   ]
  }
 ],
 "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
}