Merge branch 'main' of github.com:quicophy/mdopt

docs/source/toric.ipynb

@@ -4,20 +4,20 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "# Decoding Quantum Toric Codes"
+    "# Decoding Toric Code"
    ]
   },
   {
    "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "In this experiment, we’ll use ``mdopt`` to compute the threshold of the toric code. Hereafter, we assume an independent noise model as well as perfect syndrome measurements."
+    "In this experiment, we’ll use ``mdopt`` to decode toric code. Hereafter, we assume an independent noise model as well as perfect syndrome measurements. In this example, we will mostly follow the procedure described in a similar [tutorial](https://pymatching.readthedocs.io/en/stable/toric-code-example.html) from PyMatching with the main difference being the use of ``mdopt`` for decoding."
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": 2,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -38,66 +38,93 @@
     "    XOR_RIGHT,\n",
     ")\n",
     "from examples.decoding.decoding import (\n",
-    "    css_code_checks,\n",
-    "    css_code_logicals,\n",
-    "    css_code_logicals_sites,\n",
-    "    css_code_constraint_sites,\n",
-    ")\n",
-    "from examples.decoding.decoding import (\n",
     "    apply_constraints,\n",
     "    apply_bitflip_bias,\n",
     ")\n",
     "from examples.decoding.decoding import (\n",
     "    pauli_to_mps,\n",
-    "    decode_shor,\n",
+    "    toric_code_x_checks,\n",
+    "    toric_code_x_logicals,\n",
     ")"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 27,
+   "execution_count": 6,
    "metadata": {},
    "outputs": [],
    "source": [
-    "def repetition_code(n):\n",
-    "    \"\"\"\n",
-    "    Parity check matrix of a repetition code with length n.\n",
-    "    \"\"\"\n",
-    "    row_ind, col_ind = zip(*((i, j) for i in range(n) for j in (i, (i + 1) % n)))\n",
-    "    data = np.ones(2 * n, dtype=np.uint8)\n",
-    "    return csc_matrix((data, (row_ind, col_ind)))\n",
+    "lattice_size = 4\n",
+    "num_logicals = 2\n",
+    "num_sites = 2 * (lattice_size ** 2) + num_logicals\n",
     "\n",
     "\n",
-    "def toric_code_x_checks(L):\n",
-    "    \"\"\"\n",
-    "    Sparse check matrix for the X stabilisers of a toric code with\n",
-    "    lattice size L, constructed as the hypergraph product of\n",
-    "    two repetition codes.\n",
-    "    \"\"\"\n",
-    "    Hr = repetition_code(L)\n",
-    "    H = hstack(\n",
-    "        [kron(Hr, eye(Hr.shape[1])), kron(eye(Hr.shape[0]), Hr.T)], dtype=np.uint8\n",
-    "    )\n",
-    "    H.data = H.data % 2\n",
-    "    H.eliminate_zeros()\n",
-    "    checks = csc_matrix(H).toarray()\n",
-    "    return [list(np.nonzero(check)[0]) for check in checks]\n",
-    "\n",
+    "checks = toric_code_x_checks(lattice_size)\n",
+    "logicals = toric_code_x_logicals(lattice_size)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "[[0, 4, 16, 19],\n",
+       " [1, 5, 16, 17],\n",
+       " [2, 6, 17, 18],\n",
+       " [3, 7, 18, 19],\n",
+       " [4, 8, 20, 23],\n",
+       " [5, 9, 20, 21],\n",
+       " [6, 10, 21, 22],\n",
+       " [7, 11, 22, 23],\n",
+       " [8, 12, 24, 27],\n",
+       " [9, 13, 24, 25],\n",
+       " [10, 14, 25, 26],\n",
+       " [11, 15, 26, 27],\n",
+       " [0, 12, 28, 31],\n",
+       " [1, 13, 28, 29],\n",
+       " [2, 14, 29, 30],\n",
+       " [3, 15, 30, 31]]"
+      ]
+     },
+     "execution_count": 7,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "checks"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 8,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "[[0, 1, 2, 3], [16, 20, 24, 28]]"
+      ]
+     },
+     "execution_count": 8,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "logicals"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 27,
+   "metadata": {},
+   "outputs": [],
+   "source": [
     "\n",
-    "def toric_code_x_logicals(L):\n",
-    "    \"\"\"\n",
-    "    Sparse binary matrix with each row corresponding to an X logical operator\n",
-    "    of a toric code with lattice size L. Constructed from the\n",
-    "    homology groups of the repetition codes using the Kunneth\n",
-    "    theorem.\n",
-    "    \"\"\"\n",
-    "    H1 = csc_matrix(([1], ([0], [0])), shape=(1, L), dtype=np.uint8)\n",
-    "    H0 = csc_matrix(np.ones((1, L), dtype=np.uint8))\n",
-    "    x_logicals = block_diag([kron(H1, H0), kron(H0, H1)])\n",
-    "    x_logicals.data = x_logicals.data % 2\n",
-    "    x_logicals.eliminate_zeros()\n",
-    "    logicals = csc_matrix(x_logicals).toarray()\n",
-    "    return [list(np.nonzero(logical)[0]) for logical in logicals]\n",
     "\n",
     "\n",
     "def toric_code_constraint_sites(L):\n",
@@ -870,7 +897,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.8"
+   "version": "3.10.13"
   },
   "orig_nbformat": 4,
   "vscode": {

examples/decoding/classical.py

@@ -149,6 +149,9 @@
             success = 0
 
         failures.append(1 - success)
+        logging.info(
+            f"Finished experiment {l} for NUM_BITS={NUM_BITS}, CHI_MAX={CHI_MAX_CONTRACTOR}, PROB_ERROR={PROB_ERROR}"
+        )
 
     failures_statistics[(NUM_BITS, CHI_MAX_CONTRACTOR, PROB_ERROR)] = failures
     failures_key = (

examples/decoding/decoding.py

@@ -14,6 +14,7 @@
 from matrex import msro
 from opt_einsum import contract
 from more_itertools import powerset
+from scipy.sparse import eye, csc_matrix, hstack, kron, block_diag
 from qecstruct import (
     Rng,
     BinarySymmetricChannel,
@@ -935,3 +936,121 @@ def decode_shor(
         return "Z", logical
     if np.argmax(logical) == 3:
         return "Y", logical
+
+
+def repetition_code(num_sites: int) -> csc_matrix:
+    """
+    Generate the parity check matrix for a repetition code of length `num_sites`.
+
+    A repetition code is a simple error-correcting code in which the same data bit
+    is repeated multiple times to ensure error detection and correction. The parity
+    check matrix created by this function represents each data bit and its
+    immediate neighbor, allowing for the detection of single-bit errors in the
+    code.
+    This code is taken from PyMatching library.
+
+    Parameters
+    ----------
+    num_sites : int
+        The number of bits in the repetition code. This is also the length of the
+        code.
+
+    Returns
+    -------
+    scipy.sparse.csc_matrix
+        The parity check matrix of the repetition code. This matrix is returned
+        as a sparse matrix in compressed sparse column format, where each row
+        represents a parity check involving two bits (the current and the next,
+        with wrap-around).
+
+    Examples
+    --------
+    >>> from scipy.sparse import csc_matrix
+    >>> num_sites = 3
+    >>> matrix = repetition_code(num_sites)
+    >>> matrix.toarray()
+    array([[1, 1, 0],
+           [0, 1, 1],
+           [1, 0, 1]], dtype=uint8)
+    """
+
+    row_ind, col_ind = zip(
+        *((i, j) for i in range(num_sites) for j in (i, (i + 1) % num_sites))
+    )
+    data = np.ones(2 * num_sites, dtype=np.uint8)
+    return csc_matrix((data, (row_ind, col_ind)))
+
+
+def toric_code_x_checks(lattice_size: int) -> List[List[int]]:
+    """
+    Generate a sparse check matrix for the X stabilizers of a toric code.
+
+    The toric code is defined on a 2D lattice and this function generates the X stabilizers
+    based on the hypergraph product of two repetition codes, which is key to detecting
+    bit-flip errors in a quantum error correction setting.
+    This code is taken from PyMatching library.
+
+    Parameters
+    ----------
+    L : int
+        The lattice size, which determines the dimensions of the resulting check matrix.
+        This is also used to define the size of each repetition code.
+
+    Returns
+    -------
+    list of lists of int
+        A list of lists, where each inner list contains the indices of non-zero entries
+        (qubits involved) in each row of the check matrix, representing an X stabilizer.
+
+    Notes
+    -----
+    The function constructs the parity check matrix by creating two repetition codes and
+    taking the Kronecker product with identity matrices to align the checks correctly across
+    the 2D lattice, treating each row as an X stabilizer of the toric code.
+    """
+
+    Hr = repetition_code(lattice_size)
+    H = hstack(
+        [kron(Hr, eye(Hr.shape[1])), kron(eye(Hr.shape[0]), Hr.T)], dtype=np.uint8
+    )
+    H.data = H.data % 2
+    H.eliminate_zeros()
+    checks = csc_matrix(H).toarray()
+    return [list(np.nonzero(check)[0]) for check in checks]
+
+
+def toric_code_x_logicals(lattice_size: int) -> List[List[int]]:
+    """
+    Generate a sparse binary matrix representing the X logical operators for a toric code.
+
+    This function constructs logical operators based on the homology groups of the repetition
+    codes, using the Kunneth theorem. Each row of the resulting matrix represents an X logical
+    operator, which spans non-trivial loops around the torus in the toric code.
+    This code is taken from PyMatching library.
+
+    Parameters
+    ----------
+    L : int
+        The lattice size, which determines the dimensions of each repetition code used in the
+        construction of the logical operators.
+
+    Returns
+    -------
+    list of lists of int
+        A list of lists where each inner list contains the indices of non-zero entries
+        in each row of the logical operator matrix.
+
+    Notes
+    -----
+    The X logical operators are constructed by taking the Kronecker product of specific
+    sparse matrices representing minimal non-trivial cycles of the repetition codes,
+    reflecting the topological nature of the toric code.
+    """
+
+    H1 = csc_matrix(([1], ([0], [0])), shape=(1, lattice_size), dtype=np.uint8)
+    H0 = csc_matrix(np.ones((1, lattice_size), dtype=np.uint8))
+    x_logicals = block_diag([kron(H1, H0), kron(H0, H1)])
+    x_logicals.data = x_logicals.data % 2
+    x_logicals.eliminate_zeros()
+    logicals = csc_matrix(x_logicals).toarray()
+    return [list(np.nonzero(logical)[0]) for logical in logicals]