Add boundary-aware overloads for canonicalize_for_rounds and sliding window decoder

docs/sphinx/api/qec/cpp_api.rst

@@ -82,12 +82,14 @@ Parity Check Matrix Utilities
 .. doxygenfunction:: cudaq::qec::generate_random_pcm(std::size_t, std::size_t, std::size_t, int, std::mt19937_64 &&);
 .. doxygenfunction:: cudaq::qec::generate_timelike_sparse_detector_matrix(std::uint32_t num_syndromes_per_round, std::uint32_t num_rounds, bool include_first_round = false)
 .. doxygenfunction:: cudaq::qec::generate_timelike_sparse_detector_matrix(std::uint32_t num_syndromes_per_round, std::uint32_t num_rounds, std::vector<std::int64_t> first_round_matrix)
-.. 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_pcm_for_rounds(const cudaqx::tensor<uint8_t> &, std::uint32_t, std::uint32_t, std::uint32_t, bool, bool, std::uint32_t);
 .. 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 std::vector<std::vector<std::uint32_t>> &, 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_from_sparse_vec(const std::vector<std::int64_t>& sparse_vec, std::size_t num_rows, std::size_t num_cols)
 .. doxygenfunction:: cudaq::qec::pcm_is_sorted(const cudaqx::tensor<uint8_t> &, std::uint32_t);
+.. doxygenfunction:: cudaq::qec::pcm_is_sorted(const std::vector<std::vector<std::uint32_t>> &, std::uint32_t, std::uint32_t);
 .. doxygenfunction:: cudaq::qec::pcm_to_sparse_vec(const cudaqx::tensor<uint8_t>& pcm)
 .. 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 &&);

docs/sphinx/api/qec/sliding_window_api.rst

@@ -121,6 +121,9 @@
         - `window_size` (int): The number of rounds of syndrome data in each window. (Defaults to 1.)
         - `step_size` (int): The number of rounds to advance the window by each time. (Defaults to 1.)
         - `num_syndromes_per_round` (int): The number of syndromes per round. (Must be provided.)
+        - `num_boundary_syndromes` (int): The number of boundary syndromes, i.e. the number of 
+          detectors in the first and last round of the memory experiment. (Defaults to 0, meaning all 
+          layers have `num_syndromes_per_round` detectors.)
         - `straddle_start_round` (bool): When forming a window, should error
           mechanisms that span the start round and any preceding rounds be included? (Defaults to False.)
         - `straddle_end_round` (bool): When forming a window, should error

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

@@ -63,23 +63,50 @@ struct detector_error_model {
   /// 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. Columns sharing the same detector AND observable signature are
-  /// merged, with their rates composed so the resulting model matches the
-  /// input model. By default, zero-syndrome columns that still flip an
-  /// observable (undetectable logical errors) are retained so the model's
-  /// observable-flip probability is preserved. Set @p
-  /// remove_zero_syndrome_errors to true to drop all columns with no detector
-  /// signature, which is appropriate when the canonicalized DEM is consumed
-  /// only for round-based decoding (where such columns carry no syndrome).
-  ///
+  /// @brief Put the detector_error_matrix into canonical form for round-based
+  /// decoding: topologically order the columns and merge columns that share the
+  /// same detector AND observable signature, composing their rates so the
+  /// canonicalized model matches the input.
+  /// @param num_syndromes_per_round Number of syndromes per round, used to
+  /// order columns by the rounds their detectors span; 0 disables the per-round
+  /// key.
+  /// @param remove_zero_syndrome_errors If false (default), zero-syndrome
+  /// columns that still flip an observable (undetectable logical errors) are
+  /// retained so the observable-flip probability is preserved; if true, all
+  /// columns with no detector signature are dropped (appropriate when the DEM
+  /// is consumed only for round-based decoding).
   /// @note Canonicalization does not preserve cross-column exclusivity
-  /// structure. Each output column is given a fresh unique error id and is
-  /// treated as independent of every other column; any `error_ids` correlation
-  /// present in the input model is discarded.
+  /// structure: each output column is given a fresh unique error id and is
+  /// treated as independent, so any `error_ids` correlation in the input model
+  /// is discarded.
   void canonicalize_for_rounds(uint32_t num_syndromes_per_round,
                                bool remove_zero_syndrome_errors = false);
+
+  /// @brief Boundary-aware overload of canonicalize_for_rounds for
+  /// memory-experiment DEMs whose first and last detector layers (the
+  /// boundaries) are narrower than the interior layers
+  /// @param num_syndromes_per_round Interior-layer width (syndromes per
+  /// interior round).
+  /// @param num_boundary_syndromes Width of the leading and trailing boundary
+  /// layers: the first and last @p num_boundary_syndromes rows form the
+  /// boundary rounds and each intermediate block of @p num_syndromes_per_round
+  /// rows is an interior round.
+  /// @param remove_zero_syndrome_errors Same meaning as in the scalar overload.
+  /// @throws std::invalid_argument if @p num_syndromes_per_round is 0 or
+  /// @p num_boundary_syndromes > @p num_syndromes_per_round.
+  void canonicalize_for_rounds(uint32_t num_syndromes_per_round,
+                               uint32_t num_boundary_syndromes,
+                               bool remove_zero_syndrome_errors);
+
+private:
+  /// Shared implementation of the canonicalize_for_rounds overloads. Given the
+  /// sparse detector matrix and a precomputed topological @p column_order,
+  /// merges columns sharing a full (detector, observable) signature and
+  /// reorders/reduces the matrices accordingly.
+  void canonicalize_for_rounds_impl(
+      const std::vector<std::vector<std::uint32_t>> &row_indices,
+      const std::vector<std::uint32_t> &column_order,
+      bool remove_zero_syndrome_errors);
 };
 
 /// Parse the Stim DEM string @p dem_text into detector/observable flip

libs/qec/include/cudaq/qec/pcm_utils.h

@@ -115,6 +115,16 @@ std::vector<std::uint32_t>
 get_sorted_pcm_column_indices(const cudaqx::tensor<uint8_t> &pcm,
                               std::uint32_t num_syndromes_per_round = 0);
 
+/// @brief Boundary-aware overload of the above: the first and last
+/// @p num_boundary_syndromes rows form the boundary rounds and each interior
+/// round spans @p num_syndromes_per_round rows. throws std::invalid_argument if
+/// @p num_syndromes_per_round is zero or @p num_boundary_syndromes > @p
+/// num_syndromes_per_round.
+std::vector<std::uint32_t> get_sorted_pcm_column_indices(
+    const std::vector<std::vector<std::uint32_t>> &row_indices,
+    std::uint32_t num_syndromes_per_round,
+    std::uint32_t num_boundary_syndromes);
+
 /// @brief Check if a PCM is sorted.
 /// @param pcm The PCM to check.
 /// @param num_syndromes_per_round The number of syndromes per round.
@@ -130,6 +140,18 @@ bool pcm_is_sorted(const cudaqx::tensor<uint8_t> &pcm,
 bool pcm_is_sorted(const std::vector<std::vector<std::uint32_t>> &sparse_pcm,
                    std::uint32_t num_syndromes_per_round = 0);
 
+/// @brief Boundary-aware overload of pcm_is_sorted: the first and last
+/// @p num_boundary_syndromes rows form the boundary rounds and each interior
+/// round spans @p num_syndromes_per_round rows (a [B | K*S | B] layout).
+/// @param sparse_pcm The sparse PCM to check (in the same format as
+/// dense_to_sparse())
+/// @param num_syndromes_per_round The interior-round width.
+/// @param num_boundary_syndromes The width of the first/last boundary rounds.
+/// @return True if the PCM is sorted for this boundary layout, false otherwise.
+bool pcm_is_sorted(const std::vector<std::vector<std::uint32_t>> &sparse_pcm,
+                   std::uint32_t num_syndromes_per_round,
+                   std::uint32_t num_boundary_syndromes);
+
 /// @brief Reorder the columns of a PCM according to the given column order.
 /// Note: this may return a subset of the columns in the original PCM if the
 /// \p column_order does not contain all of the columns in the original PCM.
@@ -172,6 +194,62 @@ simplify_pcm(const cudaqx::tensor<uint8_t> &pcm,
              const std::vector<double> &weights,
              std::uint32_t num_syndromes_per_round = 0);
 
+namespace details {
+
+/// @brief Internal helper (not part of the public API). Maps between detector
+/// rounds and detector rows for a (possibly non-uniform) round layout. The rows
+/// are laid out as [B | S | ... | S | B], where B == num_boundary_syndromes is
+/// the width of the first/last boundary rounds and S == num_syndromes_per_round
+/// is the interior width. B == 0 or B == S is the uniform layout (every round
+/// has width S).
+struct round_layout {
+  std::size_t S = 0;          ///< interior round width
+  std::size_t B = 0;          ///< boundary round width
+  std::size_t num_rows = 0;   ///< total number of detector rows
+  std::size_t num_rounds = 0; ///< number of rounds (detector layers)
+  bool boundary = false;      ///< true iff B is a distinct boundary width
+
+  round_layout() = default;
+  round_layout(std::size_t num_syndromes_per_round,
+               std::size_t num_boundary_syndromes, std::size_t total_rows)
+      : S(num_syndromes_per_round), B(num_boundary_syndromes),
+        num_rows(total_rows) {
+    boundary = (B != 0 && B != S);
+    num_rounds = boundary ? (num_rows - 2 * B) / S + 2 : num_rows / S;
+  }
+
+  /// Global row index at which round @p r begins; round_start(num_rounds) is
+  /// the trailing sentinel (== num_rows).
+  std::size_t round_start(std::size_t r) const {
+    if (!boundary)
+      return r * S;
+    if (r == 0)
+      return 0;
+    if (r >= num_rounds)
+      return num_rows;
+    return B + (r - 1) * S;
+  }
+
+  /// Number of detector rows in round @p r (B for the first/last round, else
+  /// S).
+  std::size_t round_width(std::size_t r) const {
+    return round_start(r + 1) - round_start(r);
+  }
+
+  /// The round that global detector row @p row belongs to.
+  std::size_t row_to_round(std::size_t row) const {
+    if (!boundary)
+      return row / S;
+    if (row < B)
+      return 0;
+    if (row >= num_rows - B)
+      return num_rounds - 1;
+    return 1 + (row - B) / S;
+  }
+};
+
+} // namespace details
+
 /// @brief Get a sub-PCM for a range of rounds. It is recommended (but not
 /// required) that you call sort_pcm_columns() before calling this function.
 /// @param pcm The PCM to get a sub-PCM for.
@@ -182,20 +260,31 @@ simplify_pcm(const cudaqx::tensor<uint8_t> &pcm,
 /// start_round (defaults to false)
 /// @param straddle_end_round Whether to include columns that straddle the
 /// end_round (defaults to false)
+/// @param num_boundary_syndromes Width of the narrower first/last boundary
+/// rounds for a non-uniform [B | K*S | B] detector layout (0 == uniform).
 /// @return A tuple with the new PCM with the columns in the range [start_round,
 /// end_round], the first column included, and the last column included.
 std::tuple<cudaqx::tensor<uint8_t>, std::uint32_t, std::uint32_t>
 get_pcm_for_rounds(const cudaqx::tensor<uint8_t> &pcm,
                    std::uint32_t num_syndromes_per_round,
                    std::uint32_t start_round, std::uint32_t end_round,
                    bool straddle_start_round = false,
-                   bool straddle_end_round = false);
+                   bool straddle_end_round = false,
+                   std::uint32_t num_boundary_syndromes = 0);
 
 /// @brief Same semantics as the overload taking a dense tensor \p pcm, but
 /// reads from \p pcm as ``sparse_binary_matrix`` so the full dense PCM is not
 /// required (only the returned sub-matrix is dense). Parameter meanings match
 /// the dense overload.
 ///
+/// @param pcm The PCM to get a sub-PCM for.
+/// @param num_syndromes_per_round The number of syndromes per round.
+/// @param start_round The start round (0-based).
+/// @param end_round The end round (0-based).
+/// @param straddle_start_round Whether to include columns that straddle the
+/// start_round (defaults to false)
+/// @param straddle_end_round Whether to include columns that straddle the
+/// end_round (defaults to false)
 /// @param pcm_is_canonical If true, the caller asserts \p pcm has
 /// sorted-unique per-group indices (i.e. is the output of
 /// `sparse_binary_matrix::canonicalize`
@@ -215,13 +304,16 @@ get_pcm_for_rounds(const cudaqx::tensor<uint8_t> &pcm,
 /// wrong round assignments. If unsure, leave the flag `false`.
 /// @return A tuple with the new PCM with the columns in the range [start_round,
 /// end_round], the first column included, and the last column included.
+/// @param num_boundary_syndromes Width of the narrower first/last boundary
+/// rounds for a non-uniform [B | K*S | B] detector layout (0 == uniform).
 std::tuple<cudaqx::tensor<uint8_t>, std::uint32_t, std::uint32_t>
 get_pcm_for_rounds(const sparse_binary_matrix &pcm,
                    std::uint32_t num_syndromes_per_round,
                    std::uint32_t start_round, std::uint32_t end_round,
                    bool straddle_start_round = false,
                    bool straddle_end_round = false,
-                   bool pcm_is_canonical = false);
+                   bool pcm_is_canonical = false,
+                   std::uint32_t num_boundary_syndromes = 0);
 
 /// @brief Generate a random PCM with the given parameters.
 ///