Add Parity Check Matrix generation and utility functions (NVIDIA#147)

This PR introduces a CUDA-Q QEC Detector Error Model type (similar to
Stim's Detector Error Model but different).

A Detector Error Model can be generated from CUDA-Q circuits and then
fed to CUDA-Q QEC decoders as inputs so that the decoders know how to
decode the syndromes that come from those CUDA-Q circuits.

The main implementation and interfaces are done in C++. They are also
exposed to Python via pybind bindings.

In addition to the core functionality listed above, various PCM utility
functions (like column sorting, canonicalizing, combining, etc.) that
are used to create the DEM are also exposed to the user.

---------

Signed-off-by: Ben Howe <[email protected]>

Author: bmhowe23

docs/sphinx/api/qec/cpp_api.rst

@@ -22,6 +22,15 @@ Code
 .. doxygenclass:: cudaq::qec::surface_code::surface_code
     :members:
 
+Detector Error Model
+====================
+
+.. doxygenstruct:: cudaq::qec::detector_error_model
+    :members:
+
+.. doxygenfunction:: cudaq::qec::dem_from_memory_circuit(const code &, operation, std::size_t, cudaq::noise_model &)
+.. doxygenfunction:: cudaq::qec::x_dem_from_memory_circuit(const code &, operation, std::size_t, cudaq::noise_model &)
+.. doxygenfunction:: cudaq::qec::z_dem_from_memory_circuit(const code &, operation, std::size_t, cudaq::noise_model &)
 
 Decoder Interfaces
 ==================
@@ -42,6 +51,21 @@ NVIDIA QLDPC Decoder
 
 .. include:: nv_qldpc_decoder_api.rst
 
+Parity Check Matrix Utilities
+=============================
+
+.. doxygenfunction:: cudaq::qec::dense_to_sparse(const cudaqx::tensor<uint8_t> &)
+.. doxygenfunction:: cudaq::qec::generate_random_pcm(std::size_t, std::size_t, std::size_t, int, std::mt19937_64 &&);
+.. doxygenfunction:: cudaq::qec::get_pcm_for_rounds(const cudaqx::tensor<uint8_t> &, std::uint32_t, std::uint32_t, std::uint32_t, bool, bool);
+.. doxygenfunction:: cudaq::qec::get_sorted_pcm_column_indices(const std::vector<std::vector<std::uint32_t>> &, std::uint32_t);
+.. doxygenfunction:: cudaq::qec::get_sorted_pcm_column_indices(const cudaqx::tensor<uint8_t> &, std::uint32_t);
+.. doxygenfunction:: cudaq::qec::pcm_extend_to_n_rounds(const cudaqx::tensor<uint8_t> &, std::size_t, std::uint32_t);
+.. doxygenfunction:: cudaq::qec::pcm_is_sorted(const cudaqx::tensor<uint8_t> &, std::uint32_t);
+.. doxygenfunction:: cudaq::qec::reorder_pcm_columns(const cudaqx::tensor<uint8_t> &, const std::vector<std::uint32_t> &, uint32_t, uint32_t);
+.. doxygenfunction:: cudaq::qec::shuffle_pcm_columns(const cudaqx::tensor<uint8_t> &, std::mt19937_64 &&);
+.. doxygenfunction:: cudaq::qec::simplify_pcm(const cudaqx::tensor<uint8_t> &, const std::vector<double> &, std::uint32_t);
+.. doxygenfunction:: cudaq::qec::sort_pcm_columns(const cudaqx::tensor<uint8_t> &, std::uint32_t);
+
 Common
 =============
 

docs/sphinx/api/qec/python_api.rst

@@ -10,6 +10,16 @@ Code
 .. autoclass:: cudaq_qec.Code
     :members:
 
+Detector Error Model
+====================
+
+.. autoclass:: cudaq_qec.DetectorErrorModel
+    :members:
+
+.. autofunction:: cudaq_qec.dem_from_memory_circuit
+.. autofunction:: cudaq_qec.x_dem_from_memory_circuit
+.. autofunction:: cudaq_qec.z_dem_from_memory_circuit
+
 Decoder Interfaces
 ==================
 
@@ -35,3 +45,16 @@ Common
 .. autofunction:: cudaq_qec.sample_memory_circuit
 
 .. autofunction:: cudaq_qec.sample_code_capacity
+
+Parity Check Matrix Utilities
+-------------
+
+.. autofunction:: cudaq_qec.generate_random_pcm
+.. autofunction:: cudaq_qec.get_pcm_for_rounds
+.. autofunction:: cudaq_qec.get_sorted_pcm_column_indices
+.. autofunction:: cudaq_qec.pcm_extend_to_n_rounds
+.. autofunction:: cudaq_qec.pcm_is_sorted
+.. autofunction:: cudaq_qec.reorder_pcm_columns
+.. autofunction:: cudaq_qec.shuffle_pcm_columns
+.. autofunction:: cudaq_qec.simplify_pcm
+.. autofunction:: cudaq_qec.sort_pcm_columns

docs/sphinx/examples/qec/python/circuit_level_noise.py

@@ -27,7 +27,10 @@
 nShots = 3
 nRounds = 4
 
-# error probabily
+# Uncomment for repeatability
+# cudaq.set_random_seed(13)
+
+# error probability
 p = 0.01
 noise = cudaq.NoiseModel()
 noise.add_all_qubit_channel("x", cudaq.Depolarization2(p), 1)
@@ -37,14 +40,18 @@
 # our expected measurement in this state is 0
 expected_value = 0
 
+# For large runs, set verbose to False to suppress output
+verbose = nShots <= 10
+
 # sample the steane memory circuit with noise on each cx gate
 # reading out the syndromes after each stabilizer round (xor'd against the previous)
 # and readout out the data qubits at the end of the experiment
 syndromes, data = qec.sample_memory_circuit(steane, statePrep, nShots, nRounds,
                                             noise)
-print("From sample function:\n")
-print("syndromes:\n", syndromes)
-print("data:\n", data)
+if verbose:
+    print("From sample function:\n")
+    print("syndromes:\n", syndromes)
+    print("data:\n", data)
 
 # Get a decoder
 decoder = qec.get_decoder("single_error_lut", H)
@@ -54,7 +61,8 @@
 logical_measurements = (Lz @ data.transpose()) % 2
 # only one logical qubit, so do not need the second axis
 logical_measurements = logical_measurements.flatten()
-print("LMz:\n", logical_measurements)
+if verbose:
+    print("LMz:\n", logical_measurements)
 
 # organize data by shot and round if desired
 syndromes = syndromes.reshape((nShots, nRounds, syndromes.shape[1]))
@@ -63,33 +71,39 @@
 # through the stabilizer rounds
 pauli_frame = np.array([0, 0], dtype=np.uint8)
 for shot in range(0, nShots):
-    print("shot:", shot)
+    if verbose:
+        print("shot:", shot)
     for syndrome in syndromes[shot]:
-        print("syndrome:", syndrome)
+        if verbose:
+            print("syndrome:", syndrome)
         # Decode the syndrome
         results = decoder.decode(syndrome)
         convergence = results.converged
         result = results.result
         data_prediction = np.array(result, dtype=np.uint8)
 
         # see if the decoded result anti-commutes with the observables
-        print("decode result:", data_prediction)
+        if verbose:
+            print("decode result:", data_prediction)
         decoded_observables = (observables @ data_prediction) % 2
-        print("decoded_observables:", decoded_observables)
+        if verbose:
+            print("decoded_observables:", decoded_observables)
 
         # update pauli frame
         pauli_frame = (pauli_frame + decoded_observables) % 2
-        print("pauli frame:", pauli_frame)
+        if verbose:
+            print("pauli frame:", pauli_frame)
 
     # after pauli frame has tracked corrections through the rounds
     # apply the pauli frame correction to the measurement, and see
     # if this matches the state we intended to prepare
     # We prepared |0>, so we check if logical measurement Mz + Pf_X = 0
     corrected_mz = (logical_measurements[shot] + pauli_frame[0]) % 2
-    print("Expected value:", expected_value)
-    print("Corrected value:", corrected_mz)
+    if verbose:
+        print("Expected value:", expected_value)
+        print("Corrected value:", corrected_mz)
     if (corrected_mz != expected_value):
         nLogicalErrors += 1
 
 # Count how many shots the decoder failed to correct the errors
-print("Number of logical errors:", nLogicalErrors)
+print(f"Number of logical errors out of {nShots} shots: {nLogicalErrors}")

libs/core/include/cuda-qx/core/heterogeneous_map.h

@@ -1,5 +1,5 @@
 /****************************************************************-*- C++ -*-****
- * Copyright (c) 2024 NVIDIA Corporation & Affiliates.                         *
+ * Copyright (c) 2024 - 2025 NVIDIA Corporation & Affiliates.                  *
  * All rights reserved.                                                        *
  *                                                                             *
  * This source code and the accompanying materials are made available under    *
@@ -183,6 +183,10 @@ class heterogeneous_map {
   /// @return The number of key-value pairs in the map
   std::size_t size() const { return items.size(); }
 
+  /// @brief Check if the map is empty
+  /// @return true if the map is empty, false otherwise
+  bool empty() const { return items.empty(); }
+
   /// @brief Check if the map contains a key
   /// @param key The key to check
   /// @return true if the key exists, false otherwise

libs/core/include/cuda-qx/core/kwargs_utils.h

@@ -67,6 +67,9 @@ inline heterogeneous_map hetMapFromKwargs(const py::kwargs &kwargs) {
       } else {
         throw std::runtime_error("Unsupported array data type in kwargs.");
       }
+    } else if (py::isinstance<py::dict>(value)) {
+      // Check if it is a kwargs dict
+      result.insert(key, hetMapFromKwargs(value.cast<py::dict>()));
     } else {
       throw std::runtime_error(
           "Invalid python type for mapping kwargs to a heterogeneous_map.");
@@ -77,14 +80,25 @@ inline heterogeneous_map hetMapFromKwargs(const py::kwargs &kwargs) {
 }
 
 template <typename T>
-tensor<T> toTensor(const py::array_t<T> &H) {
+tensor<T> toTensor(const py::array_t<T> &H, bool perform_pcm_checks = false) {
   py::buffer_info buf = H.request();
 
   if (buf.ndim >= 1 && buf.strides[0] == buf.itemsize) {
     throw std::runtime_error("toTensor: data must be in row-major order, but "
                              "column-major order was detected.");
   }
 
+  if (perform_pcm_checks) {
+    if (buf.itemsize != sizeof(uint8_t)) {
+      throw std::runtime_error(
+          "Parity check matrix must be an array of uint8_t.");
+    }
+
+    if (buf.ndim != 2) {
+      throw std::runtime_error("Parity check matrix must be 2-dimensional.");
+    }
+  }
+
   // Create a vector of the array dimensions
   std::vector<std::size_t> shape;
   for (py::ssize_t d : buf.shape) {
@@ -96,4 +110,12 @@ tensor<T> toTensor(const py::array_t<T> &H) {
   tensor_H.borrow(static_cast<T *>(buf.ptr), shape);
   return tensor_H;
 }
+
+/// @brief Convert a py::array_t<uint8_t> to a tensor<uint8_t>. This is the same
+/// as toTensor, but with additional checks.
+template <typename T>
+tensor<T> pcmToTensor(const py::array_t<T> &H) {
+  return toTensor(H, /*perform_pcm_checks=*/true);
+}
+
 } // namespace cudaqx

libs/core/include/cuda-qx/core/tensor.h

@@ -193,7 +193,11 @@ class tensor {
       throw std::runtime_error("Dot product requires rank-2 tensors");
     }
     if (shape()[1] != other.shape()[0]) {
-      throw std::runtime_error("Invalid matrix dimensions for dot product");
+      throw std::runtime_error(
+          "Invalid matrix dimensions for dot product: matrix1 dimensions: " +
+          std::to_string(shape()[0]) + "x" + std::to_string(shape()[1]) +
+          ", matrix2 dimensions: " + std::to_string(other.shape()[0]) + "x" +
+          std::to_string(other.shape()[1]));
     }
 
     std::vector<std::size_t> result_shape = {shape()[0], other.shape()[1]};

libs/qec/include/cudaq/qec/detector_error_model.h

@@ -0,0 +1,68 @@
+/****************************************************************-*- C++ -*-****
+ * Copyright (c) 2025 NVIDIA Corporation & Affiliates.                         *
+ * All rights reserved.                                                        *
+ *                                                                             *
+ * This source code and the accompanying materials are made available under    *
+ * the terms of the Apache License 2.0 which accompanies this distribution.    *
+ ******************************************************************************/
+#pragma once
+
+#include "cuda-qx/core/tensor.h"
+#include <optional>
+
+namespace cudaq::qec {
+
+/// A detector error model (DEM) for a quantum error correction circuit. A
+/// DEM can be created from a QEC circuit and a noise model. It contains
+/// information about which errors flip which detectors. This is used by the
+/// decoder to help make predictions about observables flips.
+///
+/// Shared size parameters among the matrix types.
+/// - \p detector_error_matrix: num_detectors x num_error_mechanisms [d, e]
+/// - \p error_rates: num_error_mechanisms
+/// - \p observables_flips_matrix: num_observables x num_error_mechanisms [k, e]
+///
+/// @note The C++ API for this class may change in the future. The Python API is
+/// more likely to be backwards compatible.
+struct detector_error_model {
+  /// The detector error matrix is a specific kind of circuit-level parity-check
+  /// matrix where each row represents a detector, and each column represents
+  /// an error mechanism. The entries of this matrix are H[i,j] = 1 if detector
+  /// i is triggered by error mechanism j, and 0 otherwise.
+  cudaqx::tensor<uint8_t> detector_error_matrix;
+
+  /// The list of weights has length equal to the number of columns of
+  /// \p detector_error_matrix, which assigns a likelihood to each error
+  /// mechanism.
+  std::vector<double> error_rates;
+
+  /// The observables flips matrix is a specific kind of circuit-level parity-
+  /// check matrix where each row represents a Pauli observable, and each
+  /// column represents an error mechanism. The entries of this matrix are
+  /// O[i,j] = 1 if Pauli observable i is flipped by error mechanism j, and 0
+  /// otherwise.
+  cudaqx::tensor<uint8_t> observables_flips_matrix;
+
+  /// Error mechanism ID. From a probability perspective, each error mechanism
+  /// ID is independent of all other error mechanism ID. For all errors with
+  /// the *same* ID, only one of them can happen. That is - the errors
+  /// containing the same ID are correlated with each other.
+  std::optional<std::vector<std::size_t>> error_ids;
+
+  /// Return the number of rows in the detector_error_matrix.
+  std::size_t num_detectors() const;
+
+  /// Return the number of columns in the detector_error_matrix, error_rates,
+  /// and observables_flips_matrix.
+  std::size_t num_error_mechanisms() const;
+
+  /// Return the number of rows in the observables_flips_matrix.
+  std::size_t num_observables() const;
+
+  /// Put the detector_error_matrix into canonical form, where the rows and
+  /// columns are ordered in a way that is amenable to the round-based decoding
+  /// process.
+  void canonicalize_for_rounds(uint32_t num_syndromes_per_round);
+};
+
+} // namespace cudaq::qec

libs/qec/include/cudaq/qec/experiments.h

@@ -1,5 +1,5 @@
 /****************************************************************-*- C++ -*-****
- * Copyright (c) 2024 NVIDIA Corporation & Affiliates.                         *
+ * Copyright (c) 2024 - 2025 NVIDIA Corporation & Affiliates.                  *
  * All rights reserved.                                                        *
  *                                                                             *
  * This source code and the accompanying materials are made available under    *
@@ -8,6 +8,7 @@
 #pragma once
 
 #include "cudaq/qec/code.h"
+#include "cudaq/qec/detector_error_model.h"
 
 namespace cudaq::qec {
 
@@ -100,4 +101,36 @@ sample_memory_circuit(const code &code, std::size_t numShots,
 std::tuple<cudaqx::tensor<uint8_t>, cudaqx::tensor<uint8_t>>
 sample_memory_circuit(const code &code, std::size_t numShots,
                       std::size_t numRounds, cudaq::noise_model &noise);
+
+/// @brief Given a memory circuit setup, generate a DEM
+/// @param code QEC Code to sample
+/// @param statePrep Initial state preparation operation
+/// @param numRounds Number of stabilizer measurement rounds
+/// @param noise Noise model to apply
+/// @return Detector error model
+cudaq::qec::detector_error_model
+dem_from_memory_circuit(const code &code, operation statePrep,
+                        std::size_t numRounds, cudaq::noise_model &noise);
+
+/// @brief Given a memory circuit setup, generate a DEM for X stabilizers.
+/// @param code QEC Code to sample
+/// @param statePrep Initial state preparation operation
+/// @param numRounds Number of stabilizer measurement rounds
+/// @param noise Noise model to apply
+/// @return Detector error model
+detector_error_model x_dem_from_memory_circuit(const code &code,
+                                               operation statePrep,
+                                               std::size_t numRounds,
+                                               cudaq::noise_model &noise);
+
+/// @brief Given a memory circuit setup, generate a DEM for Z stabilizers.
+/// @param code QEC Code to sample
+/// @param statePrep Initial state preparation operation
+/// @param numRounds Number of stabilizer measurement rounds
+/// @param noise Noise model to apply
+/// @return Detector error model
+detector_error_model z_dem_from_memory_circuit(const code &code,
+                                               operation statePrep,
+                                               std::size_t numRounds,
+                                               cudaq::noise_model &noise);
 } // namespace cudaq::qec