From 7dc8c56446ba0eb3b79e43e39ee314b9e77580fc Mon Sep 17 00:00:00 2001 From: Eliot Heinrich Date: Tue, 7 Jul 2026 07:05:42 -0700 Subject: [PATCH 1/2] Added boundary-aware overloads of canonicalize_for_rounds and sliding window decoder Signed-off-by: Eliot Heinrich --- docs/sphinx/api/qec/sliding_window_api.rst | 3 + .../include/cudaq/qec/detector_error_model.h | 55 ++- libs/qec/include/cudaq/qec/pcm_utils.h | 10 + libs/qec/lib/decoders/sliding_window.cpp | 123 +++++- libs/qec/lib/decoders/sliding_window.h | 40 +- libs/qec/lib/detector_error_model.cpp | 25 ++ libs/qec/lib/experiments.cpp | 21 +- libs/qec/lib/pcm_utils.cpp | 89 +++-- libs/qec/python/bindings/py_decoder.cpp | 26 +- libs/qec/python/tests/test_sliding_window.py | 66 ++++ .../backend-specific/stim/test_qec_stim.cpp | 358 +++++++++++++++++- libs/qec/unittests/test_decoders.cpp | 26 ++ libs/qec/unittests/test_qec.cpp | 49 +++ 13 files changed, 817 insertions(+), 74 deletions(-) diff --git a/docs/sphinx/api/qec/sliding_window_api.rst b/docs/sphinx/api/qec/sliding_window_api.rst index 2790abd19..03ade7afd 100644 --- a/docs/sphinx/api/qec/sliding_window_api.rst +++ b/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 diff --git a/libs/qec/include/cudaq/qec/detector_error_model.h b/libs/qec/include/cudaq/qec/detector_error_model.h index 64729f3d0..a902da40f 100644 --- a/libs/qec/include/cudaq/qec/detector_error_model.h +++ b/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> &row_indices, + const std::vector &column_order, + bool remove_zero_syndrome_errors); }; /// Parse the Stim DEM string @p dem_text into detector/observable flip diff --git a/libs/qec/include/cudaq/qec/pcm_utils.h b/libs/qec/include/cudaq/qec/pcm_utils.h index f4ee9df90..f4e925449 100644 --- a/libs/qec/include/cudaq/qec/pcm_utils.h +++ b/libs/qec/include/cudaq/qec/pcm_utils.h @@ -115,6 +115,16 @@ std::vector get_sorted_pcm_column_indices(const cudaqx::tensor &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 get_sorted_pcm_column_indices( + const std::vector> &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. diff --git a/libs/qec/lib/decoders/sliding_window.cpp b/libs/qec/lib/decoders/sliding_window.cpp index 3ab604d4d..c6c7796ab 100644 --- a/libs/qec/lib/decoders/sliding_window.cpp +++ b/libs/qec/lib/decoders/sliding_window.cpp @@ -178,16 +178,99 @@ void sliding_window::update_rw_next_read_index() { rw_next_read_index -= num_syndromes_per_window; } +sparse_binary_matrix +sliding_window::prepare_pcm(const cudaq::qec::sparse_binary_matrix &H, + std::size_t num_syndromes_per_round, + std::size_t num_boundary_syndromes) { + const std::size_t S = num_syndromes_per_round; + const std::size_t B = num_boundary_syndromes; + + // Uniform layout (no boundary padding requested) + if (B == 0 || B == S) + return H.canonicalize().to_csc(); + + if (B > S) + throw std::invalid_argument( + "sliding_window: num_boundary_syndromes must be <= " + "num_syndromes_per_round"); + + const auto num_rows = static_cast(H.num_rows()); + // The boundary layout is [B | K*S | B]; verify it is consistent. + if (num_rows < 2 * B || (num_rows - 2 * B) % S != 0) + throw std::invalid_argument( + "sliding_window: number of PCM rows is inconsistent with the given " + "num_syndromes_per_round and num_boundary_syndromes"); + + const std::size_t pad = S - B; + // Shift rows past the initial boundary by (S - B): this zero-pads the initial + // boundary and leaves matching zero rows at the tail for the final boundary. + auto nested = H.to_nested_csc(); + for (auto &col : nested) + for (auto &r : col) + if (r >= static_cast(B)) + r += static_cast(pad); + + const std::size_t padded_num_rows = num_rows + 2 * pad; + return sparse_binary_matrix::from_nested_csc( + static_cast(padded_num_rows), + H.num_cols(), nested) + .canonicalize() + .to_csc(); +} + +std::vector +sliding_window::pad_syndrome(const std::vector &syndrome) const { + const std::size_t S = num_syndromes_per_round; + const std::size_t B = num_boundary_syndromes; + // Unpadded interior portion: syndrome[B, size - B) has length K * S. + const std::size_t num_interior = syndrome.size() - 2 * B; + + std::vector out(this->syndrome_size, 0.0); + + // Initial boundary round + for (std::size_t i = 0; i < B; ++i) + out[i] = syndrome[i]; + + // Interior rounds + for (std::size_t i = 0; i < num_interior; ++i) + out[S + i] = syndrome[B + i]; + + // Final boundary round + const std::size_t final_round_start = this->syndrome_size - S; + const std::size_t final_boundary_src = syndrome.size() - B; + for (std::size_t i = 0; i < B; ++i) + out[final_round_start + i] = syndrome[final_boundary_src + i]; + + return out; +} + +std::vector +sliding_window::pad_round(const std::vector &round) const { + // A boundary layer is [real | zeros]: the real values come first, matching + // where prepare_pcm placed the boundary rows within each padded round. + std::vector out(num_syndromes_per_round, 0.0); + std::copy(round.begin(), round.end(), out.begin()); + return out; +} + sliding_window::sliding_window(const cudaq::qec::sparse_binary_matrix &H, const cudaqx::heterogeneous_map ¶ms) // Canonical CSC is the steady-state contract for decode_window's column - // slices and for validate_inputs's per-column .front()/.back() reads. - : decoder(H.canonicalize().to_csc()) { + // slices and for validate_inputs's per-column .front()/.back() reads. When + // a boundary layout is supplied (num_boundary_syndromes set), + // prepare_pcm zero-pads the narrower boundary layers up to the interior + // width so that this->H is uniform and all the round arithmetic below stays + // uniform. + : decoder( + prepare_pcm(H, params.get("num_syndromes_per_round", 0), + params.get("num_boundary_syndromes", 0))) { // Fetch parameters from the params map. window_size = params.get("window_size", window_size); step_size = params.get("step_size", step_size); num_syndromes_per_round = params.get("num_syndromes_per_round", num_syndromes_per_round); + num_boundary_syndromes = + params.get("num_boundary_syndromes", num_boundary_syndromes); straddle_start_round = params.get("straddle_start_round", straddle_start_round); straddle_end_round = @@ -199,15 +282,25 @@ sliding_window::sliding_window(const cudaq::qec::sparse_binary_matrix &H, inner_decoder_params = params.get( "inner_decoder_params", inner_decoder_params); - // Guard the H.num_rows() / num_syndromes_per_round below. + // Guard the this->H.num_rows() / num_syndromes_per_round below. if (num_syndromes_per_round == 0) throw std::invalid_argument("sliding_window constructor: " "num_syndromes_per_round must be non-zero"); - num_rounds = H.num_rows() / num_syndromes_per_round; + // this->H is the (possibly padded) uniform matrix, so every round has + // num_syndromes_per_round rows. + num_rounds = this->H.num_rows() / num_syndromes_per_round; num_windows = (num_rounds - window_size) / step_size + 1; num_syndromes_per_window = num_syndromes_per_round * window_size; + // Syndromes handed to decode() are in the original (unpadded) layout; the + // padded layout adds (num_syndromes_per_round - num_boundary_syndromes) + // rows at each of the two boundaries. + unpadded_syndrome_size = this->syndrome_size; + if (boundary_padding_active()) + unpadded_syndrome_size -= + 2 * (num_syndromes_per_round - num_boundary_syndromes); + validate_inputs(); // this->H is canonical CSC (ctor init list), so skip the per-call @@ -246,6 +339,15 @@ sliding_window::sliding_window(const cudaq::qec::sparse_binary_matrix &H, } decoder_result sliding_window::decode(const std::vector &syndrome) { + // Boundary-layout inputs are padded to the uniform layout and re-dispatched: + // a full unpadded block via pad_syndrome, a single boundary layer via + // pad_round. Interior layers and already-padded inputs fall through. + if (boundary_padding_active()) { + if (syndrome.size() == unpadded_syndrome_size) + return decode(pad_syndrome(syndrome)); + if (syndrome.size() == num_boundary_syndromes) + return decode(pad_round(syndrome)); + } if (syndrome.size() == this->syndrome_size) { auto t0 = std::chrono::high_resolution_clock::now(); CUDA_QEC_DBG("Decoding whole block"); @@ -314,6 +416,19 @@ std::vector sliding_window::decode_batch( CUDA_QEC_DBG("Returning empty decoder_result (no syndrome)"); return {}; } + // Boundary-layout inputs are padded up to the uniform layout, then + // re-dispatched: a full unpadded block via pad_syndrome, a single boundary + // detector layer via pad_round. Interior layers pass straight through. + if (boundary_padding_active()) { + const std::size_t sz = syndromes[0].size(); + if (sz == unpadded_syndrome_size || sz == num_boundary_syndromes) { + std::vector> padded(syndromes.size()); + for (std::size_t s = 0; s < syndromes.size(); ++s) + padded[s] = (sz == unpadded_syndrome_size) ? pad_syndrome(syndromes[s]) + : pad_round(syndromes[s]); + return decode_batch(padded); + } + } if (syndromes[0].size() == this->syndrome_size) { CUDA_QEC_DBG("Decoding whole block"); // Decode the whole thing, iterating over windows manually. diff --git a/libs/qec/lib/decoders/sliding_window.h b/libs/qec/lib/decoders/sliding_window.h index 826a991d8..7e7a23709 100644 --- a/libs/qec/lib/decoders/sliding_window.h +++ b/libs/qec/lib/decoders/sliding_window.h @@ -27,8 +27,10 @@ class sliding_window : public decoder { std::size_t window_size = 1; /// The number of rounds to advance the window by each time. std::size_t step_size = 1; - /// The number of syndromes per round. + /// The number of syndromes per round (the interior-layer width). std::size_t num_syndromes_per_round = 0; + /// The number of boundary measurements. + std::size_t num_boundary_syndromes = 0; /// When forming a window, should error mechanisms that span the start round /// and any preceding rounds be included? bool straddle_start_round = false; @@ -46,6 +48,8 @@ class sliding_window : public decoder { std::size_t num_windows = 0; std::size_t num_rounds = 0; std::size_t num_syndromes_per_window = 0; + /// Number of detector rows in an unpadded (boundary-layout) syndrome. + std::size_t unpadded_syndrome_size = 0; std::size_t num_rounds_since_last_decode = 0; std::vector> inner_decoders; std::vector first_columns; @@ -108,15 +112,43 @@ class sliding_window : public decoder { /// @brief Decode a single window (internal helper) void decode_window(); + /// @brief Whether the decoder is zero-padding narrower boundary layers up to + /// the interior width (i.e. a boundary-layout was supplied). + bool boundary_padding_active() const { + return num_boundary_syndromes != 0 && + num_boundary_syndromes != num_syndromes_per_round; + } + + /// @brief Pad an all-rounds syndrome from the unpadded boundary + /// layout up to the uniform padded layout by inserting zeros for the + /// boundary-layer padding rows. + std::vector pad_syndrome(const std::vector &syndrome) const; + + /// @brief Pad a single boundary detector layer up to the interior width by + /// appending zeros. + std::vector pad_round(const std::vector &round) const; + + /// @brief Zero-pad the two boundary layers of @p H up to the interior width + /// so all rounds are uniform. The first/last @p num_boundary_syndromes + /// rows are the narrower boundary layers and interior layers have + /// @p num_syndromes_per_round rows. Returns the canonical CSC form expected + /// by the base decoder. If no boundary padding is needed (@p + /// num_boundary_syndromes is 0 or equals @p num_syndromes_per_round), this + /// is just H.canonicalize().to_csc(). + static sparse_binary_matrix + prepare_pcm(const cudaq::qec::sparse_binary_matrix &H, + std::size_t num_syndromes_per_round, + std::size_t num_boundary_syndromes); + public: /// @brief Constructor /// @param H The full parity check matrix for all rounds /// @param params A heterogeneous map containing required parameters: /// - window_size: Size of each decoding window (in rounds) /// - step_size: Step size between consecutive windows (in rounds) - /// - num_rounds: Total number of rounds - /// - num_syndromes_per_round: Number of syndromes per round - /// - inner_decoder: Name of the inner decoder to use + /// - num_syndromes_per_round: Number of syndromes per (interior) round + /// - num_boundary_syndromes: Boundary-layer width (0 if uniform) + /// - inner_decoder_name: Name of the inner decoder to use /// - inner_decoder_params: Parameters for the inner decoder (optional) sliding_window(const cudaq::qec::sparse_binary_matrix &H, const cudaqx::heterogeneous_map ¶ms); diff --git a/libs/qec/lib/detector_error_model.cpp b/libs/qec/lib/detector_error_model.cpp index 79a89b798..dc63d449a 100644 --- a/libs/qec/lib/detector_error_model.cpp +++ b/libs/qec/lib/detector_error_model.cpp @@ -154,6 +154,31 @@ void detector_error_model::canonicalize_for_rounds( auto row_indices = dense_to_sparse(detector_error_matrix); auto column_order = get_sorted_pcm_column_indices(row_indices, num_syndromes_per_round); + canonicalize_for_rounds_impl(row_indices, column_order, + remove_zero_syndrome_errors); +} + +void detector_error_model::canonicalize_for_rounds( + uint32_t num_syndromes_per_round, uint32_t num_boundary_syndromes, + bool remove_zero_syndrome_errors) { + // A boundary wider than the interior would misassign rounds silently. + if (num_boundary_syndromes > num_syndromes_per_round) + throw std::invalid_argument( + "canonicalize_for_rounds: num_boundary_syndromes (" + + std::to_string(num_boundary_syndromes) + + ") must be <= num_syndromes_per_round (" + + std::to_string(num_syndromes_per_round) + ")"); + auto row_indices = dense_to_sparse(detector_error_matrix); + auto column_order = get_sorted_pcm_column_indices( + row_indices, num_syndromes_per_round, num_boundary_syndromes); + canonicalize_for_rounds_impl(row_indices, column_order, + remove_zero_syndrome_errors); +} + +void detector_error_model::canonicalize_for_rounds_impl( + const std::vector> &row_indices, + const std::vector &column_order, + bool remove_zero_syndrome_errors) { const std::size_t num_obs = this->num_observables(); const auto num_cols = column_order.size(); const bool has_error_ids = diff --git a/libs/qec/lib/experiments.cpp b/libs/qec/lib/experiments.cpp index 9dbb9d0d1..f1240e3dc 100644 --- a/libs/qec/lib/experiments.cpp +++ b/libs/qec/lib/experiments.cpp @@ -510,11 +510,22 @@ dem_from_memory_circuit(const code &code, operation statePrep, } dem.detector_error_matrix = std::move(selectedRows); - const std::size_t numReturnSynPerRound = - (keep_z_stabilizers ? numZStabs : 0) + - (keep_x_stabilizers ? numXStabs : 0); - dem.canonicalize_for_rounds(numReturnSynPerRound, - /*remove_zero_syndrome_errors=*/true); + if (keep_x_stabilizers && keep_z_stabilizers) { + // Full both-basis DEM: the boundary layers carry only the fixed basis and + // are narrower than the interior, so canonicalize is boundary-aware. + // Boundary width = fixed-basis count (numZStabs for a Z-basis prep, else + // numXStabs). + const uint32_t numBoundary = is_z_prep ? numZStabs : numXStabs; + dem.canonicalize_for_rounds(numXStabs + numZStabs, numBoundary, + /*remove_zero_syndrome_errors=*/true); + } else { + // A single-basis DEM is uniform (every layer has the same width) + const std::size_t numReturnSynPerRound = + (keep_z_stabilizers ? numZStabs : 0) + + (keep_x_stabilizers ? numXStabs : 0); + dem.canonicalize_for_rounds(numReturnSynPerRound, + /*remove_zero_syndrome_errors=*/true); + } return dem; } diff --git a/libs/qec/lib/pcm_utils.cpp b/libs/qec/lib/pcm_utils.cpp index 8a5dc90c6..db244ef77 100644 --- a/libs/qec/lib/pcm_utils.cpp +++ b/libs/qec/lib/pcm_utils.cpp @@ -187,28 +187,27 @@ void select_pcm_columns_for_round_range( namespace cudaq::qec { -/// @brief Return a vector of column indices that would sort the PCM columns -/// in topological order. -/// @param row_indices For each column, a vector of row indices that have a -/// non-zero value in that column. -/// @param num_syndromes_per_round The number of syndromes per round. (Defaults -/// to 0, which means that no secondary per-round sorting will occur.) -/// @details This function tries to make a matrix that is close to a block -/// diagonal matrix from its input. If \p num_syndromes_per_round is > 0, then -/// the columns are first sorted by rounds numbers in which the checks are -/// performed. Columns are then sorted by the index of the first non-zero entry -/// in the column, and if those match, then they are sorted by the index of the -/// last non-zero entry in the column. This ping pong continues for the indices -/// of the second non-zero element and the second-to-last non-zero element, and -/// so forth. -std::vector get_sorted_pcm_column_indices( - const std::vector> &row_indices, - std::uint32_t num_syndromes_per_round) { +namespace { +/// Topological column sort shared by the scalar and boundary-aware +/// get_sorted_pcm_column_indices overloads. @p round_of maps a detector row to +/// its round; when @p use_rounds is false the round key is skipped. +/// +/// This function tries to make a matrix that is close to a block diagonal +/// matrix from its input. If @p use_rounds is true, then the columns are first +/// sorted by the round numbers in which the checks are performed. Columns are +/// then sorted by the index of the first non-zero entry in the column, and if +/// those match, then they are sorted by the index of the last non-zero entry +/// in the column. This ping pong continues for the indices of the second +/// non-zero element and the second-to-last non-zero element, and so forth. +template +std::vector sorted_pcm_column_indices_impl( + const std::vector> &row_indices, bool use_rounds, + RoundOf round_of) { std::vector column_order(row_indices.size()); std::iota(column_order.begin(), column_order.end(), 0); std::sort(column_order.begin(), column_order.end(), - [&row_indices, num_syndromes_per_round](const std::uint32_t &a, - const std::uint32_t &b) { + [&row_indices, use_rounds, &round_of](const std::uint32_t &a, + const std::uint32_t &b) { const auto &a_vec = row_indices[a]; const auto &b_vec = row_indices[b]; @@ -227,14 +226,12 @@ std::vector get_sorted_pcm_column_indices( auto b_it_head = b_vec.begin(); auto b_it_tail = b_vec.end() - 1; - // First sort by the span of rounds that the errors appear in. We - // can only do this sorting if we know how many syndromes per - // round. - if (num_syndromes_per_round > 0) { - auto a_first_round = *a_it_head / num_syndromes_per_round; - auto a_last_round = *a_it_tail / num_syndromes_per_round; - auto b_first_round = *b_it_head / num_syndromes_per_round; - auto b_last_round = *b_it_tail / num_syndromes_per_round; + // First sort by the span of rounds that the errors appear in. + if (use_rounds) { + auto a_first_round = round_of(*a_it_head); + auto a_last_round = round_of(*a_it_tail); + auto b_first_round = round_of(*b_it_head); + auto b_last_round = round_of(*b_it_tail); if (a_first_round != b_first_round) return a_first_round < b_first_round; if (a_last_round != b_last_round) @@ -291,6 +288,44 @@ std::vector get_sorted_pcm_column_indices( return column_order; } +} // namespace + +std::vector get_sorted_pcm_column_indices( + const std::vector> &row_indices, + std::uint32_t num_syndromes_per_round) { + return sorted_pcm_column_indices_impl( + row_indices, /*use_rounds=*/num_syndromes_per_round > 0, + // Only invoked when use_rounds is true, so the divisor is non-zero. + [num_syndromes_per_round](std::uint32_t r) { + return r / num_syndromes_per_round; + }); +} + +std::vector get_sorted_pcm_column_indices( + const std::vector> &row_indices, + std::uint32_t num_syndromes_per_round, + std::uint32_t num_boundary_syndromes) { + if (num_syndromes_per_round == 0) + throw std::invalid_argument( + "get_sorted_pcm_column_indices: num_syndromes_per_round must be " + "non-zero when num_boundary_syndromes is specified"); + if (num_boundary_syndromes > num_syndromes_per_round) + throw std::invalid_argument( + "get_sorted_pcm_column_indices: num_boundary_syndromes (" + + std::to_string(num_boundary_syndromes) + + ") must be <= num_syndromes_per_round (" + + std::to_string(num_syndromes_per_round) + ")"); + // Detector rows map to rounds as: the first num_boundary_syndromes rows + // are round 0, then each num_syndromes_per_round block is an interior round. + return sorted_pcm_column_indices_impl( + row_indices, /*use_rounds=*/true, + [num_syndromes_per_round, + num_boundary_syndromes](std::uint32_t r) -> std::uint32_t { + if (r < num_boundary_syndromes) + return 0; + return 1 + (r - num_boundary_syndromes) / num_syndromes_per_round; + }); +} bool pcm_is_sorted(const std::vector> &sparse_pcm, std::uint32_t num_syndromes_per_round) { diff --git a/libs/qec/python/bindings/py_decoder.cpp b/libs/qec/python/bindings/py_decoder.cpp index 69af10ad9..76340a891 100644 --- a/libs/qec/python/bindings/py_decoder.cpp +++ b/libs/qec/python/bindings/py_decoder.cpp @@ -793,7 +793,8 @@ void bindDecoder(nb::module_ &mod) { The number of observables in the detector error model )pbdoc") .def("canonicalize_for_rounds", - &detector_error_model::canonicalize_for_rounds, + static_cast( + &detector_error_model::canonicalize_for_rounds), R"pbdoc( Canonicalize the detector error model for a given number of rounds. @@ -812,7 +813,28 @@ void bindDecoder(nb::module_ &mod) { correlation in the input model is discarded. )pbdoc", nb::arg("num_syndromes_per_round"), - nb::arg("remove_zero_syndrome_errors") = false); + nb::arg("remove_zero_syndrome_errors") = false) + .def( + "canonicalize_for_rounds", + static_cast( + &detector_error_model::canonicalize_for_rounds), + R"pbdoc( + Boundary-aware canonicalization for memory-experiment DEMs whose + first and last detector layers (the boundaries) are narrower than + the interior layers. The first ``num_boundary_syndromes`` detector + rows form the initial round, each subsequent block of + ``num_syndromes_per_round`` rows is an interior round, and the + trailing ``num_boundary_syndromes`` rows form the final round. + This makes the round-based column ordering respect the true rounds + even when the boundary width differs from the interior width. + + ``remove_zero_syndrome_errors`` behaves as in the two-argument + overload. Raises ``ValueError`` if ``num_syndromes_per_round`` is + zero or ``num_boundary_syndromes`` exceeds + ``num_syndromes_per_round``. + )pbdoc", + nb::arg("num_syndromes_per_round"), nb::arg("num_boundary_syndromes"), + nb::arg("remove_zero_syndrome_errors") = false); qecmod.def("dem_from_stim_text", &dem_from_stim_text, R"pbdoc( diff --git a/libs/qec/python/tests/test_sliding_window.py b/libs/qec/python/tests/test_sliding_window.py index a34c937c9..321eeb106 100644 --- a/libs/qec/python/tests/test_sliding_window.py +++ b/libs/qec/python/tests/test_sliding_window.py @@ -94,6 +94,72 @@ def test_sliding_window_1(decoder_name, batched, num_rounds, num_windows): assert num_mismatches == 0 +@pytest.mark.parametrize("code_name", ['steane', 'repetition', 'surface_code']) +@pytest.mark.parametrize("num_windows", [1, 2, 3]) +def test_sliding_window_boundary_layout(code_name, num_windows): + # A full both-basis DEM has a non-uniform detector layout: the first and + # last detector layers (the boundaries) carry only the fixed basis, so they + # are narrower than the interior layers. The sliding window handles this via + # num_boundary_syndromes, and must agree with a full decoder. The + # repetition code has no X stabilizers, so it exercises the degenerate + # uniform (boundary width == interior width) path. + cudaq.set_random_seed(13) + code = qec.get_code(code_name, distance=3) + numX = code.get_num_x_stabilizers() + numZ = code.get_num_z_stabilizers() + interior = numX + numZ + num_boundary = numZ # prep0 => Z is the fixed (boundary) basis + + p = 0.001 + noise = cudaq.NoiseModel() + noise.add_all_qubit_channel("x", cudaq.Depolarization2(p), 1) + + # Choose num_rounds so the padded layout gives the requested num_windows + # with window_size = padded_rounds - num_windows + 1, step_size = 1. + num_rounds = 6 + dem = qec.dem_from_memory_circuit(code, qec.operation.prep0, num_rounds, + noise) + H = np.asarray(dem.detector_error_matrix) + rows, cols = H.shape + # Padded (uniform) layout has one extra (interior-boundary) block at each + # boundary, i.e. num_rounds + 1 layers. + padded_rounds = (rows + 2 * (interior - num_boundary)) // interior + assert padded_rounds == num_rounds + 1 + window_size = padded_rounds - num_windows + 1 + + full_decoder = qec.get_decoder('single_error_lut', + dem.detector_error_matrix) + sw = qec.get_decoder("sliding_window", + dem.detector_error_matrix, + window_size=window_size, + step_size=1, + num_syndromes_per_round=interior, + num_boundary_syndromes=num_boundary, + straddle_start_round=False, + straddle_end_round=True, + error_rate_vec=np.array(dem.error_rates), + inner_decoder_name="single_error_lut", + inner_decoder_params={'dummy_param': 1}) + + # Inject one error mechanism per shot (a DEM column) in the unpadded + # boundary layout, and confirm the observable prediction matches the full + # decoder. + np.random.seed(13) + nShots = 300 + O = np.asarray(dem.observables_flips_matrix) + num_mismatches = 0 + for _ in range(nShots): + col = np.random.randint(0, cols) + syndrome = H[:, col].astype(np.uint8) + r_full = np.asarray(full_decoder.decode(syndrome).result) > 0.5 + r_sw = np.asarray(sw.decode(syndrome).result) > 0.5 + of_full = (O @ r_full.astype(np.uint8)) % 2 + of_sw = (O @ r_sw.astype(np.uint8)) % 2 + if not np.array_equal(of_full, of_sw): + num_mismatches += 1 + assert num_mismatches == 0 + + def test_pymatching_parallel_edges_use_observable_faults(): # Same detector syndrome with different observable flips are distinct # logical fault mechanisms. After #610 these stay as separate DEM columns, diff --git a/libs/qec/unittests/backend-specific/stim/test_qec_stim.cpp b/libs/qec/unittests/backend-specific/stim/test_qec_stim.cpp index 32ba97843..6c030a00b 100644 --- a/libs/qec/unittests/backend-specific/stim/test_qec_stim.cpp +++ b/libs/qec/unittests/backend-specific/stim/test_qec_stim.cpp @@ -11,8 +11,10 @@ #include "cudaq.h" +#include "cudaq/qec/code.h" #include "cudaq/qec/decoder.h" #include "cudaq/qec/experiments.h" +#include "cudaq/qec/patch.h" TEST(QECCodeTester, checkRepetitionNoiseStim) { @@ -602,27 +604,27 @@ TEST(QECCodeTester, checkDemFromMemoryCircuit) { "1..............................", ".1.............................", "..1............................", - "1..1...........................", - ".1..1..........................", - "..1..1.........................", - "......1..1.....................", - ".......1..1....................", - "........1..1...................", - "...1........1..................", - "....1........1.................", - ".....1........1................", + "1.....1........................", + ".1.....1.......................", + "..1.....1......................", + "...1.....1.....................", + "....1.....1....................", + ".....1.....1...................", + "......1.....1..................", + ".......1.....1.................", + "........1.....1................", ".........1.....1...............", "..........1.....1..............", "...........1.....1.............", - "............1.....1............", - ".............1.....1...........", - "..............1.....1..........", - "...............1.....1.........", - "................1.....1........", - ".................1.....1.......", - "..................1.....1111...", - "...................1.....1.111.", - "....................1.....11.11"}; + "............1........1.........", + ".............1........1........", + "..............1........1.......", + "...............1..1............", + "................1..1...........", + ".................1..1..........", + ".....................1..1111...", + "......................1..1.111.", + ".......................1..11.11"}; // clang-format on check_matrix_bits("detector_error_matrix", dem.detector_error_matrix, expected_dem_str); @@ -653,3 +655,323 @@ TEST(QECCodeTester, checkDemFromMemoryCircuit) { check_matrix_bits("observables_flips_matrix", dem.observables_flips_matrix, expected_observables_flips_matrix_str); } + +// --------------------------------------------------------------------------- +// Test-only custom code: the Shor [[9,1,3]] code, a CSS code with 2 X-type +// and 6 Z-type stabilizers. The unequal stabilizer counts produce a +// non-uniform detector layout, exercising the boundary-aware round machinery. +// --------------------------------------------------------------------------- +namespace shor9_test { + +__qpu__ void prep0(cudaq::qec::patch logicalQubit) { + for (std::size_t i = 0; i < logicalQubit.data.size(); i++) + reset(logicalQubit.data[i]); +} + +// X-basis prep (|+...+>): a +1 eigenstate of the X-stabilizers, so X is the +// fixed (boundary) basis, exercising the X-type boundary. +__qpu__ void prepp(cudaq::qec::patch logicalQubit) { + prep0(logicalQubit); + h(logicalQubit.data); +} + +__qpu__ std::vector +stabilizer(cudaq::qec::patch logicalQubit, + const std::vector &x_stabilizers, + const std::vector &z_stabilizers) { + h(logicalQubit.ancx); + for (std::size_t xi = 0; xi < logicalQubit.ancx.size(); ++xi) + for (std::size_t di = 0; di < logicalQubit.data.size(); ++di) + if (x_stabilizers[xi * logicalQubit.data.size() + di] == 1) + cudaq::x(logicalQubit.ancx[xi], logicalQubit.data[di]); + h(logicalQubit.ancx); + + for (std::size_t zi = 0; zi < logicalQubit.ancz.size(); ++zi) + for (std::size_t di = 0; di < logicalQubit.data.size(); ++di) + if (z_stabilizers[zi * logicalQubit.data.size() + di] == 1) + cudaq::x(logicalQubit.data[di], logicalQubit.ancz[zi]); + + auto results = mz(logicalQubit.ancz, logicalQubit.ancx); + + for (std::size_t i = 0; i < logicalQubit.ancx.size(); i++) + reset(logicalQubit.ancx[i]); + for (std::size_t i = 0; i < logicalQubit.ancz.size(); i++) + reset(logicalQubit.ancz[i]); + return results; +} + +class shor9 : public cudaq::qec::code { +protected: + std::size_t get_num_data_qubits() const override { return 9; } + std::size_t get_num_ancilla_qubits() const override { return 8; } + std::size_t get_num_ancilla_x_qubits() const override { return 2; } + std::size_t get_num_ancilla_z_qubits() const override { return 6; } + std::size_t get_num_x_stabilizers() const override { return 2; } + std::size_t get_num_z_stabilizers() const override { return 6; } + +public: + shor9(const cudaqx::heterogeneous_map &) : code() { + operation_encodings.insert( + std::make_pair(cudaq::qec::operation::stabilizer_round, stabilizer)); + operation_encodings.insert( + std::make_pair(cudaq::qec::operation::prep0, prep0)); + operation_encodings.insert( + std::make_pair(cudaq::qec::operation::prepp, prepp)); + m_stabilizers = + fromPauliWords({"XXXXXXIII", "IIIXXXXXX", "ZZIIIIIII", "IZZIIIIII", + "IIIZZIIII", "IIIIZZIII", "IIIIIIZZI", "IIIIIIIZZ"}); + m_pauli_observables = fromPauliWords({"XXXIIIIII", "ZIIZIIZII"}); + } + + CUDAQ_EXTENSION_CUSTOM_CREATOR_FUNCTION( + shor9, static std::unique_ptr create( + const cudaqx::heterogeneous_map &options) { + return std::make_unique(options); + }) +}; + +CUDAQ_EXT_PT_REGISTER_TYPE(shor9) + +} // namespace shor9_test + +namespace { + +// Build the Shor [[9,1,3]] full both-basis DEM for a `num_rounds` memory +// experiment. `prep` selects the basis (prep0 => Z boundary, prepp => X). +cudaq::qec::detector_error_model +shor9_dem(std::size_t num_rounds, + cudaq::qec::operation prep = cudaq::qec::operation::prep0) { + auto shor = cudaq::qec::get_code("shor9"); + cudaq::noise_model noise; + noise.add_all_qubit_channel("mz", cudaq::bit_flip_channel(0.01)); + return cudaq::qec::dem_from_memory_circuit(*shor, prep, num_rounds, noise); +} + +// Build sliding_window parameters for a boundary-layout DEM. +cudaqx::heterogeneous_map +shor9_sliding_params(std::size_t window_size, std::size_t interior, + std::size_t numBoundary, + const std::vector &error_rates) { + cudaqx::heterogeneous_map inner_params; + inner_params.insert("dummy_param", 1); + cudaqx::heterogeneous_map params; + params.insert("window_size", window_size); + params.insert("step_size", static_cast(1)); + params.insert("num_syndromes_per_round", interior); + params.insert("num_boundary_syndromes", numBoundary); + params.insert("straddle_start_round", false); + params.insert("straddle_end_round", true); + params.insert("error_rate_vec", error_rates); + params.insert("inner_decoder_name", std::string("single_error_lut")); + params.insert("inner_decoder_params", inner_params); + return params; +} + +// For every DEM column, use it as a syndrome, decode it with `decode_fn`, and +// require the resulting observable-flip prediction to match the full decoder's. +template +void expectObservablesMatchFullDecoder( + const cudaq::qec::detector_error_model &dem, cudaq::qec::decoder &full, + DecodeFn decode_fn) { + const auto &H = dem.detector_error_matrix; + const auto &O = dem.observables_flips_matrix; + const std::size_t rows = H.shape()[0], cols = H.shape()[1]; + const std::size_t numObs = O.shape()[0]; + std::size_t mismatches = 0; + for (std::size_t col = 0; col < cols; ++col) { + std::vector syndrome(rows, 0.0); + for (std::size_t r = 0; r < rows; ++r) + syndrome[r] = H.at({r, col}); + + auto r_full = full.decode(syndrome).result; + auto r_sw = decode_fn(syndrome); + ASSERT_EQ(r_full.size(), r_sw.size()); + + // Compare observable flips O @ hard(result) (mod 2). + for (std::size_t o = 0; o < numObs; ++o) { + std::uint8_t of_full = 0, of_sw = 0; + for (std::size_t c = 0; c < r_full.size(); ++c) + if (O.at({o, c})) { + of_full ^= (r_full[c] > 0.5) ? 1 : 0; + of_sw ^= (r_sw[c] > 0.5) ? 1 : 0; + } + if (of_full != of_sw) + mismatches++; + } + } + EXPECT_EQ(mismatches, 0u) + << "sliding-window observable predictions disagree with the full decoder"; +} + +} // namespace + +TEST(QECCodeTester, checkDemFromMemoryCircuitShor9) { + auto shor = cudaq::qec::get_code("shor9"); + ASSERT_TRUE(shor != nullptr); + + const std::size_t numXStabs = shor->get_num_x_stabilizers(); + const std::size_t numZStabs = shor->get_num_z_stabilizers(); + EXPECT_EQ(numXStabs, 2u); + EXPECT_EQ(numZStabs, 6u); + + const std::size_t num_rounds = 4; + const std::uint32_t interiorWidth = numXStabs + numZStabs; + + // prep0 => Z-type boundary (numZStabs wide); prepp => X-type (numXStabs). + for (auto cfg : {std::make_pair(cudaq::qec::operation::prep0, numZStabs), + std::make_pair(cudaq::qec::operation::prepp, numXStabs)}) { + const auto prep = cfg.first; + const std::size_t numFixed = cfg.second; + + auto dem = shor9_dem(num_rounds, prep); + + // Non-uniform layout: a numFixed boundary, (num_rounds - 1) interior rounds + // of interiorWidth detectors, and a final numFixed boundary. + const std::size_t expected_rows = + 2 * numFixed + (num_rounds - 1) * interiorWidth; + EXPECT_EQ(dem.detector_error_matrix.shape()[0], expected_rows); + + // Column counts must be consistent across the three data structures. + const std::size_t num_cols = dem.detector_error_matrix.shape()[1]; + EXPECT_GT(num_cols, 0u); + EXPECT_EQ(dem.error_rates.size(), num_cols); + EXPECT_EQ(dem.observables_flips_matrix.shape()[1], num_cols); + EXPECT_EQ(dem.num_observables(), 1u); + + // Every retained column carries a nonzero rate and detector signature + // (remove_zero_syndrome_errors was requested). + for (std::size_t c = 0; c < num_cols; c++) { + EXPECT_GT(dem.error_rates[c], 0.0); + std::size_t weight = 0; + for (std::size_t r = 0; r < expected_rows; r++) + weight += dem.detector_error_matrix.at({r, c}); + EXPECT_GT(weight, 0u) << "column " << c << " has an empty syndrome"; + } + + // Columns must be ordered by their true rounds (first numFixed rows = + // round 0, then interiorWidth-wide interior rounds). + auto true_round = [&](std::uint32_t r) -> std::uint32_t { + if (r < numFixed) + return 0; + return 1 + (r - static_cast(numFixed)) / interiorWidth; + }; + + std::pair prev = {0, 0}; + for (std::size_t c = 0; c < num_cols; c++) { + std::uint32_t first_row = expected_rows, last_row = 0; + for (std::size_t r = 0; r < expected_rows; r++) + if (dem.detector_error_matrix.at({r, c})) { + first_row = std::min(first_row, r); + last_row = std::max(last_row, r); + } + std::pair key = {true_round(first_row), + true_round(last_row)}; + if (c > 0) + EXPECT_LE(prev, key) << "column " << c << " is out of true-round order"; + prev = key; + } + + // Calling the boundary-aware canonicalize again is stable. + dem.canonicalize_for_rounds(interiorWidth, static_cast(numFixed), + /*remove_zero_syndrome_errors=*/true); + EXPECT_EQ(dem.detector_error_matrix.shape()[0], expected_rows); + EXPECT_EQ(dem.detector_error_matrix.shape()[1], num_cols); + EXPECT_EQ(dem.error_rates.size(), num_cols); + } +} + +// Sliding-window decoding of the Shor [[9,1,3]] full both-basis (boundary) +// DEM. The sliding window internally zero-pads the boundary layers up to the +// interior width via num_boundary_syndromes, so its whole-block decode must +// agree with a full decoder. +TEST(QECCodeTester, checkSlidingWindowShor9Boundary) { + auto shor = cudaq::qec::get_code("shor9"); + const std::size_t numXStabs = shor->get_num_x_stabilizers(); + const std::size_t numZStabs = shor->get_num_z_stabilizers(); + const std::size_t interior = numXStabs + numZStabs; + const std::size_t num_rounds = 4; + + // prep0 => Z-type boundary; prepp => X-type boundary. + for (auto cfg : {std::make_pair(cudaq::qec::operation::prep0, numZStabs), + std::make_pair(cudaq::qec::operation::prepp, numXStabs)}) { + const auto prep = cfg.first; + const std::size_t numBoundary = cfg.second; + + auto dem = shor9_dem(num_rounds, prep); + const std::size_t rows = dem.detector_error_matrix.shape()[0]; + // Padded (uniform) layout adds (interior - numBoundary) rows at each + // boundary => num_rounds + 1 uniform rounds. + const std::size_t padded_rounds = + (rows + 2 * (interior - numBoundary)) / interior; + ASSERT_EQ(padded_rounds, num_rounds + 1); + + auto full = + cudaq::qec::decoder::get("single_error_lut", dem.detector_error_matrix); + // A single window spanning all rounds -- should match the full decoder. + auto sw = cudaq::qec::decoder::get( + "sliding_window", dem.detector_error_matrix, + shor9_sliding_params(padded_rounds, interior, numBoundary, + dem.error_rates)); + + expectObservablesMatchFullDecoder( + dem, *full, [&](const std::vector &syndrome) { + return sw->decode(syndrome).result; + }); + } +} + +// Real-time (streaming) sliding-window decoding of the Shor boundary layout: +// detector layers are fed one at a time. The two boundary layers arrive with +// numZStabs (6) values and the interior layers with numXStabs+numZStabs (8) +// values; the decoder pads the boundary layers on the fly. The streamed result +// must match a full decoder. +TEST(QECCodeTester, checkSlidingWindowShor9Streaming) { + auto shor = cudaq::qec::get_code("shor9"); + const std::size_t numXStabs = shor->get_num_x_stabilizers(); + const std::size_t numZStabs = shor->get_num_z_stabilizers(); + const std::size_t interior = numXStabs + numZStabs; + const std::size_t num_rounds = 4; + + // prep0 => Z-type boundary; prepp => X-type boundary. + for (auto cfg : {std::make_pair(cudaq::qec::operation::prep0, numZStabs), + std::make_pair(cudaq::qec::operation::prepp, numXStabs)}) { + const auto prep = cfg.first; + const std::size_t numBoundary = cfg.second; + + auto dem = shor9_dem(num_rounds, prep); + const std::size_t rows = dem.detector_error_matrix.shape()[0]; + const std::size_t padded_rounds = + (rows + 2 * (interior - numBoundary)) / interior; + + // Detector-layer sizes: [numBoundary | interior*(padded_rounds-2) | + // numBoundary]. + std::vector layer_sizes(padded_rounds, interior); + layer_sizes.front() = numBoundary; + layer_sizes.back() = numBoundary; + + auto full = + cudaq::qec::decoder::get("single_error_lut", dem.detector_error_matrix); + // A genuinely sliding configuration: window of 2 rounds, stepping by 1. + auto sw = cudaq::qec::decoder::get( + "sliding_window", dem.detector_error_matrix, + shor9_sliding_params(/*window_size=*/2, interior, numBoundary, + dem.error_rates)); + + expectObservablesMatchFullDecoder( + dem, *full, [&](const std::vector &syndrome) { + // Feed the syndrome one detector layer at a time (variable widths). + cudaq::qec::decoder_result streamed; + std::size_t off = 0; + for (auto ls : layer_sizes) { + std::vector layer(syndrome.begin() + off, + syndrome.begin() + off + ls); + off += ls; + auto r = sw->decode(layer); + if (!r.result.empty()) + streamed = std::move(r); // final layer yields the result + } + EXPECT_FALSE(streamed.result.empty()); + return streamed.result; + }); + } +} diff --git a/libs/qec/unittests/test_decoders.cpp b/libs/qec/unittests/test_decoders.cpp index 0f2f4c8d5..890257a05 100644 --- a/libs/qec/unittests/test_decoders.cpp +++ b/libs/qec/unittests/test_decoders.cpp @@ -1190,3 +1190,29 @@ TEST(SlidingWindowDecoder, BaseStreamingCopiesFirstRoundDetectors) { << "First-round detector copy runs, but the sliding window is not full " "yet so no final correction is committed."; } + +TEST(SlidingWindowDecoder, PreparePcmRejectsBadBoundaryLayout) { + // prepare_pcm runs in the constructor's member initializer list, so an + // inconsistent boundary layout must throw during construction. + auto params = [](std::size_t S, std::size_t B) { + cudaqx::heterogeneous_map p; + p.insert("window_size", std::size_t{1}); + p.insert("step_size", std::size_t{1}); + p.insert("num_syndromes_per_round", S); + p.insert("num_boundary_syndromes", B); + p.insert("error_rate_vec", std::vector{}); + p.insert("inner_decoder_name", std::string("single_error_lut")); + p.insert("inner_decoder_params", cudaqx::heterogeneous_map{}); + return p; + }; + + // Boundary wider than the interior (B > S). + cudaqx::tensor H4({4, 2}); + EXPECT_THROW(cudaq::qec::decoder::get("sliding_window", H4, params(2, 4)), + std::invalid_argument); + + // Row count inconsistent with a [B | K*S | B] layout (B < S). + cudaqx::tensor H10({10, 2}); + EXPECT_THROW(cudaq::qec::decoder::get("sliding_window", H10, params(8, 3)), + std::invalid_argument); +} diff --git a/libs/qec/unittests/test_qec.cpp b/libs/qec/unittests/test_qec.cpp index 47326f084..3284d234f 100644 --- a/libs/qec/unittests/test_qec.cpp +++ b/libs/qec/unittests/test_qec.cpp @@ -1701,6 +1701,40 @@ TEST(PCMUtilsTester, checkGetPCMForRounds) { } } +TEST(PCMUtilsTester, checkGetSortedColumnIndicesWithBoundary) { + // Boundary layout: interior width S = 8, boundary width B = 6. The + // boundary-aware overload must map detector rows to rounds as + // round(r) = (r < B) ? 0 : 1 + (r - B) / S, + // whereas the plain overload uses r / S. The two disagree on the transition + // rows [B, S), which is what this test pins down. + const std::uint32_t S = 8, B = 6; + + // Three columns whose correct (boundary) round order differs from the order + // the plain r/S mapping produces: + // c0 = {5, 12} : boundary span (0, 1), r/S span (0, 1) + // c1 = {6} : boundary span (1, 1), r/S span (0, 0) + // c2 = {0} : boundary span (0, 0), r/S span (0, 0) + std::vector> cols = {{5, 12}, {6}, {0}}; + + // Boundary order sorts by (first_round, last_round): c2(0,0) < c0(0,1) < + // c1(1,1), i.e. indices {2, 0, 1}. + auto order = cudaq::qec::get_sorted_pcm_column_indices(cols, S, B); + EXPECT_EQ(order, (std::vector{2, 0, 1})); + + // The plain r/S mapping would put c1 in round 0, yielding a different order + // ({2, 1, 0}); confirm the boundary argument actually changes the result. + auto uniform_order = cudaq::qec::get_sorted_pcm_column_indices(cols, S); + EXPECT_EQ(uniform_order, (std::vector{2, 1, 0})); + EXPECT_NE(order, uniform_order); + + // Invalid boundary arguments must throw rather than divide by zero or + // silently misassign rounds. + EXPECT_THROW(cudaq::qec::get_sorted_pcm_column_indices(cols, 0, B), + std::invalid_argument); + EXPECT_THROW(cudaq::qec::get_sorted_pcm_column_indices(cols, S, S + 1), + std::invalid_argument); +} + TEST(PCMUtilsTester, checkShufflePCMColumns) { std::size_t n_rounds = 4; std::size_t n_errs_per_round = 30; @@ -1987,6 +2021,21 @@ TEST(DetectorErrorModelTest, FailureOnEmptyErrorRatesCanonicalize) { EXPECT_THROW(dem.canonicalize_for_rounds(2), std::runtime_error); } +TEST(DetectorErrorModelTest, CanonicalizeBoundaryRejectsWideBoundary) { + // The boundary-aware overload must reject a boundary wider than the interior + // (num_boundary_syndromes > num_syndromes_per_round) + cudaq::qec::detector_error_model dem; + dem.detector_error_matrix = cudaqx::tensor({4, 2}); + dem.observables_flips_matrix = cudaqx::tensor({1, 2}); + dem.error_rates = {0.1, 0.2}; + + EXPECT_THROW( + dem.canonicalize_for_rounds(/*num_syndromes_per_round=*/2, + /*num_boundary_syndromes=*/3, + /*remove_zero_syndrome_errors=*/true), + std::invalid_argument); +} + TEST(DetectorErrorModelTest, CanonicalizeWithoutErrorIds) { // This test covers the std::numeric_limits::max() branch // when has_error_ids is false and duplicate columns need to be merged From b46b0708c138c66f61fd9dbd4d90186f61ad0ff2 Mon Sep 17 00:00:00 2001 From: Eliot Heinrich Date: Wed, 8 Jul 2026 15:56:20 -0700 Subject: [PATCH 2/2] Buffered sliding window syndromes in non-uniform buffer Signed-off-by: Eliot Heinrich --- docs/sphinx/api/qec/cpp_api.rst | 4 +- libs/qec/include/cudaq/qec/pcm_utils.h | 86 +++- libs/qec/lib/decoders/sliding_window.cpp | 488 ++++++----------------- libs/qec/lib/decoders/sliding_window.h | 85 +--- libs/qec/lib/pcm_utils.cpp | 78 ++-- libs/qec/unittests/test_decoders.cpp | 4 +- libs/qec/unittests/test_qec.cpp | 56 +++ 7 files changed, 347 insertions(+), 454 deletions(-) diff --git a/docs/sphinx/api/qec/cpp_api.rst b/docs/sphinx/api/qec/cpp_api.rst index fd752c534..f9e317487 100644 --- a/docs/sphinx/api/qec/cpp_api.rst +++ b/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 first_round_matrix) -.. doxygenfunction:: cudaq::qec::get_pcm_for_rounds(const cudaqx::tensor &, std::uint32_t, std::uint32_t, std::uint32_t, bool, bool); +.. doxygenfunction:: cudaq::qec::get_pcm_for_rounds(const cudaqx::tensor &, 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::uint32_t); +.. doxygenfunction:: cudaq::qec::get_sorted_pcm_column_indices(const std::vector> &, std::uint32_t, std::uint32_t); .. doxygenfunction:: cudaq::qec::get_sorted_pcm_column_indices(const cudaqx::tensor &, std::uint32_t); .. doxygenfunction:: cudaq::qec::pcm_extend_to_n_rounds(const cudaqx::tensor &, std::size_t, std::uint32_t); .. doxygenfunction:: cudaq::qec::pcm_from_sparse_vec(const std::vector& sparse_vec, std::size_t num_rows, std::size_t num_cols) .. doxygenfunction:: cudaq::qec::pcm_is_sorted(const cudaqx::tensor &, std::uint32_t); +.. doxygenfunction:: cudaq::qec::pcm_is_sorted(const std::vector> &, std::uint32_t, std::uint32_t); .. doxygenfunction:: cudaq::qec::pcm_to_sparse_vec(const cudaqx::tensor& pcm) .. doxygenfunction:: cudaq::qec::reorder_pcm_columns(const cudaqx::tensor &, const std::vector &, uint32_t, uint32_t); .. doxygenfunction:: cudaq::qec::shuffle_pcm_columns(const cudaqx::tensor &, std::mt19937_64 &&); diff --git a/libs/qec/include/cudaq/qec/pcm_utils.h b/libs/qec/include/cudaq/qec/pcm_utils.h index f4e925449..e106176e2 100644 --- a/libs/qec/include/cudaq/qec/pcm_utils.h +++ b/libs/qec/include/cudaq/qec/pcm_utils.h @@ -140,6 +140,18 @@ bool pcm_is_sorted(const cudaqx::tensor &pcm, bool pcm_is_sorted(const std::vector> &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> &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. @@ -182,6 +194,62 @@ simplify_pcm(const cudaqx::tensor &pcm, const std::vector &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. @@ -192,6 +260,8 @@ simplify_pcm(const cudaqx::tensor &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, std::uint32_t, std::uint32_t> @@ -199,13 +269,22 @@ get_pcm_for_rounds(const cudaqx::tensor &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` @@ -225,13 +304,16 @@ get_pcm_for_rounds(const cudaqx::tensor &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, 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. /// diff --git a/libs/qec/lib/decoders/sliding_window.cpp b/libs/qec/lib/decoders/sliding_window.cpp index c6c7796ab..0195dbfd3 100644 --- a/libs/qec/lib/decoders/sliding_window.cpp +++ b/libs/qec/lib/decoders/sliding_window.cpp @@ -17,11 +17,21 @@ namespace cudaq::qec { void sliding_window::validate_inputs() { uint32_t num_rows = H.num_rows(); - if (window_size < 1 || window_size > num_rounds) { + if (num_boundary_syndromes > num_syndromes_per_round) + throw std::invalid_argument( + "sliding_window constructor: num_boundary_syndromes must be <= " + "num_syndromes_per_round"); + // The detector rows must form a [B | K*S | B] layout. + if (num_rows < 2 * num_boundary_syndromes || + (num_rows - 2 * num_boundary_syndromes) % num_syndromes_per_round != 0) + throw std::invalid_argument( + "sliding_window constructor: number of PCM rows is inconsistent with " + "the given num_syndromes_per_round and num_boundary_syndromes"); + if (window_size < 1 || window_size > num_detector_layers) { throw std::invalid_argument( fmt::format("sliding_window constructor: window_size ({}) must " - "be between 1 and num_rounds ({})", - window_size, num_rounds)); + "be between 1 and the number of detector layers ({})", + window_size, num_detector_layers)); } if (step_size < 1 || step_size > window_size) { throw std::invalid_argument( @@ -29,21 +39,16 @@ void sliding_window::validate_inputs() { "be between 1 and window_size ({})", step_size, window_size)); } - if ((num_rounds - window_size) % step_size != 0) { + if ((num_detector_layers - window_size) % step_size != 0) { throw std::invalid_argument( - fmt::format("sliding_window constructor: num_rounds - " + fmt::format("sliding_window constructor: detector layers - " "window_size ({}) must be divisible by step_size ({})", - num_rounds - window_size, step_size)); + num_detector_layers - window_size, step_size)); } if (num_syndromes_per_round == 0) { throw std::invalid_argument("sliding_window constructor: " "num_syndromes_per_round must be non-zero"); } - if (num_rows % num_syndromes_per_round != 0) { - throw std::invalid_argument( - "sliding_window constructor: Number of rows in H must be divisible " - "by num_syndromes_per_round"); - } if (inner_decoder_name.empty()) { throw std::invalid_argument( "sliding_window constructor: inner_decoder_name must be non-empty"); @@ -59,211 +64,49 @@ void sliding_window::validate_inputs() { // Enforce topological column order. Ctor-time materialization only. if (!cudaq::qec::pcm_is_sorted(this->H.to_nested_csc(), - this->num_syndromes_per_round)) { + this->num_syndromes_per_round, + this->num_boundary_syndromes)) { throw std::invalid_argument("sliding_window constructor: PCM must be " "sorted. See cudaq::qec::simplify_pcm."); } } /// Helper function to initialize the window. -/// @param num_syndromes The number of syndromes to initialize the window for. -/// This will be 1 for non-batched mode. -void sliding_window::initialize_window(std::size_t num_syndromes) { +/// @param batch_size The number of independent syndromes (the batch size) to +/// initialize the window for. This will be 1 for non-batched mode. +void sliding_window::initialize_window(std::size_t batch_size) { // Initialize the syndrome mods and rw_results. auto t0 = std::chrono::high_resolution_clock::now(); window_proc_times_arr.fill(0.0); - syndrome_mods.resize(num_syndromes); - for (std::size_t s = 0; s < num_syndromes; ++s) { + syndrome_mods.resize(batch_size); + for (std::size_t s = 0; s < batch_size; ++s) { syndrome_mods[s].clear(); syndrome_mods[s].resize(this->syndrome_size); } rw_results.clear(); - rw_results.resize(num_syndromes); - for (std::size_t s = 0; s < num_syndromes; ++s) { + rw_results.resize(batch_size); + for (std::size_t s = 0; s < batch_size; ++s) { rw_results[s].converged = true; // Gets set to false if we fail to decode rw_results[s].result.resize(this->block_size); } - rolling_window.resize(num_syndromes); - for (std::size_t s = 0; s < num_syndromes; ++s) { - rolling_window[s].clear(); - rolling_window[s].resize(num_syndromes_per_window); - } window_proc_times.resize(num_windows); std::fill(window_proc_times.begin(), window_proc_times.end(), 0.0); - rw_next_write_index = 0; - rw_next_read_index = 0; - rw_filled = 0; - num_rounds_since_last_decode = 0; + this->batch_size = batch_size; + window_rounds.clear(); + window_start_round = 0; + rounds_since_last_reset = 0; + num_windows_decoded = 0; CUDA_QEC_DBG("Initializing window"); auto t1 = std::chrono::high_resolution_clock::now(); window_proc_times_arr[WindowProcTimes::INITIALIZE_WINDOW] = std::chrono::duration(t1 - t0).count() * 1000; } -/// Helper function to add a single syndrome to the rolling window (circular -/// buffer). -void sliding_window::add_syndrome_to_rolling_window( - const std::vector &syndrome, std::size_t syndrome_index, - bool update_next_write_index) { - // This assumes that the syndrome size evenly divides into the rolling - // window (of length num_syndromes_per_window), so verify that here. - if (num_syndromes_per_window % syndrome.size() != 0) { - throw std::invalid_argument( - fmt::format("add_syndrome_to_rolling_window: syndrome " - "size ({}) must evenly divide into the rolling " - "window size ({})", - syndrome.size(), num_syndromes_per_window)); - } - std::copy(syndrome.begin(), syndrome.end(), - rolling_window[syndrome_index].begin() + rw_next_write_index); - if (update_next_write_index) { - rw_next_write_index += syndrome.size(); - if (rw_next_write_index >= num_syndromes_per_window) - rw_next_write_index = 0; - } -} - -/// Helper function to add a batch of syndromes to the rolling window -/// (circular buffer). -void sliding_window::add_syndromes_to_rolling_window( - const std::vector> &syndromes) { - // Set update_next_write_index to false in the loop because we will update - // it once at the end. - for (std::size_t s = 0; s < syndromes.size(); ++s) { - add_syndrome_to_rolling_window(syndromes[s], s, - /*update_next_write_index=*/false); - if (syndromes[s].size() != syndromes[0].size()) { - throw std::invalid_argument( - fmt::format("add_syndromes_to_rolling_window: syndrome " - "size ({}) must be the same as the first syndrome " - "size ({})", - syndromes[s].size(), syndromes[0].size())); - } - } - rw_next_write_index += syndromes[0].size(); - if (rw_next_write_index >= num_syndromes_per_window) - rw_next_write_index = 0; -} - -/// Helper function to get a single syndrome from the rolling window -/// (unwrapping a circular buffer). -std::vector -sliding_window::get_syndrome_from_rolling_window(std::size_t syndrome_index) { - std::vector syndrome(num_syndromes_per_window); - // Copy from rw_next_read_index to the end of the buffer. - std::copy(rolling_window[syndrome_index].begin() + rw_next_read_index, - rolling_window[syndrome_index].end(), syndrome.begin()); - // Copy from the beginning of the rolling window to rw_next_read_index. - std::copy(rolling_window[syndrome_index].begin(), - rolling_window[syndrome_index].begin() + rw_next_read_index, - syndrome.end() - rw_next_read_index); - return syndrome; -} - -/// Helper function to get a batch of syndromes from the rolling window -/// (unwrapping a circular buffer). -std::vector> -sliding_window::get_syndromes_from_rolling_window() { - std::vector> syndromes(rolling_window.size()); - for (std::size_t s = 0; s < rolling_window.size(); ++s) { - syndromes[s] = get_syndrome_from_rolling_window(s); - } - return syndromes; -} - -/// Helper function to update the read index for the rolling window. -void sliding_window::update_rw_next_read_index() { - rw_next_read_index += step_size * num_syndromes_per_round; - if (rw_next_read_index >= num_syndromes_per_window) - rw_next_read_index -= num_syndromes_per_window; -} - -sparse_binary_matrix -sliding_window::prepare_pcm(const cudaq::qec::sparse_binary_matrix &H, - std::size_t num_syndromes_per_round, - std::size_t num_boundary_syndromes) { - const std::size_t S = num_syndromes_per_round; - const std::size_t B = num_boundary_syndromes; - - // Uniform layout (no boundary padding requested) - if (B == 0 || B == S) - return H.canonicalize().to_csc(); - - if (B > S) - throw std::invalid_argument( - "sliding_window: num_boundary_syndromes must be <= " - "num_syndromes_per_round"); - - const auto num_rows = static_cast(H.num_rows()); - // The boundary layout is [B | K*S | B]; verify it is consistent. - if (num_rows < 2 * B || (num_rows - 2 * B) % S != 0) - throw std::invalid_argument( - "sliding_window: number of PCM rows is inconsistent with the given " - "num_syndromes_per_round and num_boundary_syndromes"); - - const std::size_t pad = S - B; - // Shift rows past the initial boundary by (S - B): this zero-pads the initial - // boundary and leaves matching zero rows at the tail for the final boundary. - auto nested = H.to_nested_csc(); - for (auto &col : nested) - for (auto &r : col) - if (r >= static_cast(B)) - r += static_cast(pad); - - const std::size_t padded_num_rows = num_rows + 2 * pad; - return sparse_binary_matrix::from_nested_csc( - static_cast(padded_num_rows), - H.num_cols(), nested) - .canonicalize() - .to_csc(); -} - -std::vector -sliding_window::pad_syndrome(const std::vector &syndrome) const { - const std::size_t S = num_syndromes_per_round; - const std::size_t B = num_boundary_syndromes; - // Unpadded interior portion: syndrome[B, size - B) has length K * S. - const std::size_t num_interior = syndrome.size() - 2 * B; - - std::vector out(this->syndrome_size, 0.0); - - // Initial boundary round - for (std::size_t i = 0; i < B; ++i) - out[i] = syndrome[i]; - - // Interior rounds - for (std::size_t i = 0; i < num_interior; ++i) - out[S + i] = syndrome[B + i]; - - // Final boundary round - const std::size_t final_round_start = this->syndrome_size - S; - const std::size_t final_boundary_src = syndrome.size() - B; - for (std::size_t i = 0; i < B; ++i) - out[final_round_start + i] = syndrome[final_boundary_src + i]; - - return out; -} - -std::vector -sliding_window::pad_round(const std::vector &round) const { - // A boundary layer is [real | zeros]: the real values come first, matching - // where prepare_pcm placed the boundary rows within each padded round. - std::vector out(num_syndromes_per_round, 0.0); - std::copy(round.begin(), round.end(), out.begin()); - return out; -} - sliding_window::sliding_window(const cudaq::qec::sparse_binary_matrix &H, const cudaqx::heterogeneous_map ¶ms) // Canonical CSC is the steady-state contract for decode_window's column - // slices and for validate_inputs's per-column .front()/.back() reads. When - // a boundary layout is supplied (num_boundary_syndromes set), - // prepare_pcm zero-pads the narrower boundary layers up to the interior - // width so that this->H is uniform and all the round arithmetic below stays - // uniform. - : decoder( - prepare_pcm(H, params.get("num_syndromes_per_round", 0), - params.get("num_boundary_syndromes", 0))) { + // slices and for validate_inputs's per-column .front()/.back() reads. + : decoder(H.canonicalize().to_csc()) { // Fetch parameters from the params map. window_size = params.get("window_size", window_size); step_size = params.get("step_size", step_size); @@ -282,27 +125,25 @@ sliding_window::sliding_window(const cudaq::qec::sparse_binary_matrix &H, inner_decoder_params = params.get( "inner_decoder_params", inner_decoder_params); - // Guard the this->H.num_rows() / num_syndromes_per_round below. if (num_syndromes_per_round == 0) throw std::invalid_argument("sliding_window constructor: " "num_syndromes_per_round must be non-zero"); - // this->H is the (possibly padded) uniform matrix, so every round has - // num_syndromes_per_round rows. - num_rounds = this->H.num_rows() / num_syndromes_per_round; - num_windows = (num_rounds - window_size) / step_size + 1; - num_syndromes_per_window = num_syndromes_per_round * window_size; + // Treat a 0 boundary width as the uniform layout (B == S). + if (num_boundary_syndromes == 0) + num_boundary_syndromes = num_syndromes_per_round; - // Syndromes handed to decode() are in the original (unpadded) layout; the - // padded layout adds (num_syndromes_per_round - num_boundary_syndromes) - // rows at each of the two boundaries. - unpadded_syndrome_size = this->syndrome_size; - if (boundary_padding_active()) - unpadded_syndrome_size -= - 2 * (num_syndromes_per_round - num_boundary_syndromes); + const std::size_t num_rows = this->H.num_rows(); + // The [B | S...S | B] detector-layer layout drives all round<->row mapping. + layout = details::round_layout(num_syndromes_per_round, + num_boundary_syndromes, num_rows); + num_detector_layers = layout.num_rounds; // 2 boundary + K interior + num_windows = (num_detector_layers - window_size) / step_size + 1; validate_inputs(); + // Build the per-window inner decoders from the real (unpadded) sub-PCMs. The + // boundary-aware round layout is handled by get_pcm_for_rounds. // this->H is canonical CSC (ctor init list), so skip the per-call // canonicalize in get_pcm_for_rounds. for (std::size_t w = 0; w < num_windows; ++w) { @@ -310,7 +151,8 @@ sliding_window::sliding_window(const cudaq::qec::sparse_binary_matrix &H, std::size_t end_round = start_round + window_size - 1; auto [H_round, first_column, last_column] = cudaq::qec::get_pcm_for_rounds( this->H, num_syndromes_per_round, start_round, end_round, - straddle_start_round, straddle_end_round, /*pcm_is_canonical=*/true); + straddle_start_round, straddle_end_round, /*pcm_is_canonical=*/true, + num_boundary_syndromes); first_columns.push_back(first_column); // Slice the error vector to only include the current window. @@ -339,75 +181,10 @@ sliding_window::sliding_window(const cudaq::qec::sparse_binary_matrix &H, } decoder_result sliding_window::decode(const std::vector &syndrome) { - // Boundary-layout inputs are padded to the uniform layout and re-dispatched: - // a full unpadded block via pad_syndrome, a single boundary layer via - // pad_round. Interior layers and already-padded inputs fall through. - if (boundary_padding_active()) { - if (syndrome.size() == unpadded_syndrome_size) - return decode(pad_syndrome(syndrome)); - if (syndrome.size() == num_boundary_syndromes) - return decode(pad_round(syndrome)); - } - if (syndrome.size() == this->syndrome_size) { - auto t0 = std::chrono::high_resolution_clock::now(); - CUDA_QEC_DBG("Decoding whole block"); - // Decode the whole thing, iterating over windows manually. - decoder_result result; - std::vector syndrome_round(num_syndromes_per_round); - for (std::size_t r = 0; r < num_rounds; ++r) { - std::copy(syndrome.begin() + r * num_syndromes_per_round, - syndrome.begin() + (r + 1) * num_syndromes_per_round, - syndrome_round.begin()); - result = decode(syndrome_round); - // Note: result will be empty until the final loop iteration. - } - auto t1 = std::chrono::high_resolution_clock::now(); - std::chrono::duration diff = t1 - t0; - CUDA_QEC_INFO("Whole block time: {:.3f} ms", diff.count() * 1000); - return result; - } - // Else we're receiving a single round. - if (rw_filled == 0) { - initialize_window(/*num_syndromes=*/1); - } - if (this->rw_filled == num_syndromes_per_window) { - auto t0 = std::chrono::high_resolution_clock::now(); - CUDA_QEC_DBG("Window is full, sliding the window by one round"); - add_syndrome_to_rolling_window(syndrome, 0); - - auto t1 = std::chrono::high_resolution_clock::now(); - window_proc_times_arr[WindowProcTimes::SLIDE_WINDOW] += - std::chrono::duration(t1 - t0).count() * 1000; - } else { - // Just copy the data to the end of the rolling window. - auto t0 = std::chrono::high_resolution_clock::now(); - CUDA_QEC_DBG("Copying data to the end of the rolling window"); - add_syndrome_to_rolling_window(syndrome, 0); - this->rw_filled += num_syndromes_per_round; - auto t1 = std::chrono::high_resolution_clock::now(); - window_proc_times_arr[WindowProcTimes::COPY_DATA] += - std::chrono::duration(t1 - t0).count() * 1000; - } - num_rounds_since_last_decode++; - if (rw_filled == num_syndromes_per_window && - num_rounds_since_last_decode >= step_size) { - CUDA_QEC_DBG("Decoding window {}/{}", num_windows_decoded + 1, num_windows); - decode_window(); - num_rounds_since_last_decode = 0; - - num_windows_decoded++; - if (num_windows_decoded == num_windows) { - num_windows_decoded = 0; - rw_filled = 0; - // for (std::size_t w = 0; w < num_windows; ++w) { - // CUDA_QEC_DBG("Window {} time: {} ms", w, window_proc_times[w]); - // } - CUDA_QEC_DBG("Returning decoder_result"); - return std::move(this->rw_results[0]); - } - } - CUDA_QEC_DBG("Returning empty decoder_result"); - return decoder_result(); // empty return value + auto results = decode_batch({syndrome}); + if (results.empty()) + return decoder_result(); // empty until the final window + return std::move(results[0]); } std::vector sliding_window::decode_batch( @@ -416,69 +193,64 @@ std::vector sliding_window::decode_batch( CUDA_QEC_DBG("Returning empty decoder_result (no syndrome)"); return {}; } - // Boundary-layout inputs are padded up to the uniform layout, then - // re-dispatched: a full unpadded block via pad_syndrome, a single boundary - // detector layer via pad_round. Interior layers pass straight through. - if (boundary_padding_active()) { - const std::size_t sz = syndromes[0].size(); - if (sz == unpadded_syndrome_size || sz == num_boundary_syndromes) { - std::vector> padded(syndromes.size()); - for (std::size_t s = 0; s < syndromes.size(); ++s) - padded[s] = (sz == unpadded_syndrome_size) ? pad_syndrome(syndromes[s]) - : pad_round(syndromes[s]); - return decode_batch(padded); - } - } if (syndromes[0].size() == this->syndrome_size) { CUDA_QEC_DBG("Decoding whole block"); - // Decode the whole thing, iterating over windows manually. + // Decode the whole thing, feeding one detector layer at a time. std::vector results; std::vector> syndromes_round(syndromes.size()); - for (std::size_t r = 0; r < num_rounds; ++r) { + for (std::size_t r = 0; r < num_detector_layers; ++r) { + std::size_t round_start = layout.round_start(r); + std::size_t round_end = round_start + layout.round_width(r); for (std::size_t s = 0; s < syndromes.size(); ++s) { - syndromes_round[s].resize(num_syndromes_per_round); - std::copy(syndromes[s].begin() + r * num_syndromes_per_round, - syndromes[s].begin() + (r + 1) * num_syndromes_per_round, - syndromes_round[s].begin()); + syndromes_round[s].resize(round_end - round_start); + std::copy(syndromes[s].begin() + round_start, + syndromes[s].begin() + round_end, syndromes_round[s].begin()); } results = decode_batch(syndromes_round); } return results; } // Else we're receiving a single round. - if (rw_filled == 0) { + if (rounds_since_last_reset == 0) initialize_window(syndromes.size()); + + if (syndromes.size() != batch_size) + throw std::invalid_argument( + fmt::format("sliding_window: batch size changed mid-stream ({} vs {})", + syndromes.size(), batch_size)); + + const std::size_t expected = layout.round_width(rounds_since_last_reset); + for (const auto &r : syndromes) + if (r.size() != expected) + throw std::invalid_argument(fmt::format( + "sliding_window: round {} has width {} but expected {} for this " + "round in the boundary layout", + rounds_since_last_reset, r.size(), expected)); + + window_rounds.push_back(syndromes); + ++rounds_since_last_reset; + + if (window_rounds.size() < window_size) + return {}; + + // A full window is buffered; decode it. + CUDA_QEC_DBG("Decoding window {}/{}", num_windows_decoded + 1, num_windows); + decode_window(); + ++num_windows_decoded; + + if (num_windows_decoded == num_windows) { + // Final window decoded: hand back the accumulated results and reset. + auto results = std::move(rw_results); + window_rounds.clear(); + rounds_since_last_reset = 0; + num_windows_decoded = 0; + window_start_round = 0; + return results; } - if (this->rw_filled == num_syndromes_per_window) { - CUDA_QEC_DBG("Window is full, sliding the window by one round"); - // The window is full. Slide existing data to the left and write the new - // data at the end. - add_syndromes_to_rolling_window(syndromes); - num_rounds_since_last_decode++; - } else { - // Just copy the data to the end of the rolling window. - CUDA_QEC_DBG("Copying data to the end of the rolling window"); - add_syndromes_to_rolling_window(syndromes); - this->rw_filled += num_syndromes_per_round; - num_rounds_since_last_decode++; - } - if (rw_filled == num_syndromes_per_window && - num_rounds_since_last_decode >= step_size) { - CUDA_QEC_DBG("Decoding window {}/{}", num_windows_decoded + 1, num_windows); - decode_window(); - num_rounds_since_last_decode = 0; - num_windows_decoded++; - if (num_windows_decoded == num_windows) { - num_windows_decoded = 0; - rw_filled = 0; - // Dump the per window processing times. - // for (std::size_t w = 0; w < num_windows; ++w) { - // CUDA_QEC_DBG("Window {} time: {} ms", w, window_proc_times[w]); - // } - CUDA_QEC_DBG("Returning decoder_result"); - return std::move(this->rw_results); - } - } + + // Slide the window: drop the oldest step_size rounds. + window_rounds.erase(window_rounds.begin(), window_rounds.begin() + step_size); + window_start_round += step_size; CUDA_QEC_DBG("Returning empty decoder_result"); return std::vector(); // empty return value } @@ -490,50 +262,40 @@ std::vector sliding_window::decode_batch( void sliding_window::decode_window() { auto t0 = std::chrono::high_resolution_clock::now(); const auto &w = this->num_windows_decoded; - std::size_t syndrome_start = w * step_size * num_syndromes_per_round; - std::size_t syndrome_end = syndrome_start + num_syndromes_per_window - 1; - std::size_t syndrome_start_next_window = - (w + 1) * step_size * num_syndromes_per_round; - std::size_t syndrome_end_next_window = - syndrome_start_next_window + num_syndromes_per_round - 1; + // Detector range of window w's rounds. + std::size_t syndrome_start = layout.round_start(w * step_size); + std::size_t num_window_syndromes = + layout.round_start(w * step_size + window_size) - syndrome_start; auto t3 = std::chrono::high_resolution_clock::now(); - if (w > 0) { - // Modify the syndrome slice to account for the previous windows. - for (std::size_t s = 0; s < this->rolling_window.size(); ++s) { - std::size_t r2 = rw_next_read_index; - for (std::size_t r = 0; r < num_syndromes_per_window; ++r) { - auto &slice_val = this->rolling_window[s].at(r2); - slice_val = - static_cast(static_cast(slice_val) ^ - syndrome_mods[s].at(r + syndrome_start)); - r2++; - if (r2 >= num_syndromes_per_window) - r2 = 0; - } + std::vector> window_syndromes(batch_size); + for (std::size_t s = 0; s < batch_size; ++s) { + // Assemble each batch element's window syndrome by concatenating its + // buffered rounds + auto &syn = window_syndromes[s]; + syn.reserve(num_window_syndromes); + for (std::size_t slot = 0; slot < window_size; ++slot) + syn.insert(syn.end(), window_rounds[slot][s].begin(), + window_rounds[slot][s].end()); + if (w > 0) { + // Apply the accumulated syndrome mods from the previously committed + // windows. + for (std::size_t r = 0; r < num_window_syndromes; ++r) + syn[r] = static_cast(static_cast(syn[r]) ^ + syndrome_mods[s][syndrome_start + r]); } } auto t4 = std::chrono::high_resolution_clock::now(); - CUDA_QEC_DBG("Window {}: syndrome_start = {}, syndrome_end = {}, length1 = " - "{}, length2 = {}", - w, syndrome_start, syndrome_end, this->rolling_window[0].size(), - syndrome_end - syndrome_start + 1); - std::vector inner_results; - if (this->rolling_window.size() == 1) { - inner_results.push_back( - inner_decoders[w]->decode(get_syndrome_from_rolling_window(0))); - } else { - inner_results = - inner_decoders[w]->decode_batch(get_syndromes_from_rolling_window()); - } - // We've grabbed data from the rolling window, so we need to update the - // read index for the next call to decode_window. - update_rw_next_read_index(); + CUDA_QEC_DBG("Window {}: syndrome_start = {}, num_window_syndromes = {}", w, + syndrome_start, num_window_syndromes); + + std::vector inner_results = + inner_decoders[w]->decode_batch(window_syndromes); if (!inner_results[0].converged) { CUDA_QEC_DBG("Window {}: inner decoder failed to converge", w); } auto t5 = std::chrono::high_resolution_clock::now(); - std::vector> window_results(this->rolling_window.size()); - for (std::size_t s = 0; s < this->rolling_window.size(); ++s) { + std::vector> window_results(batch_size); + for (std::size_t s = 0; s < batch_size; ++s) { this->rw_results[s].converged &= inner_results[s].converged; cudaq::qec::convert_vec_soft_to_hard(inner_results[s].result, window_results[s]); @@ -546,7 +308,7 @@ void sliding_window::decode_window() { auto this_window_first_column = first_columns[w]; auto num_to_commit = next_window_first_column - this_window_first_column; CUDA_QEC_DBG(" Committing {} bits from window {}", num_to_commit, w); - for (std::size_t s = 0; s < this->rolling_window.size(); ++s) { + for (std::size_t s = 0; s < batch_size; ++s) { for (std::size_t c = 0; c < num_to_commit; ++c) { rw_results[s].result[c + this_window_first_column] = window_results[s][c]; @@ -555,9 +317,13 @@ void sliding_window::decode_window() { // Back out committed errors from the next window's syndrome by flipping // the rows where the corresponding H columns have a 1. Read directly off // the canonical CSC arrays. + std::size_t syndrome_start_next_window = + layout.round_start((w + 1) * step_size); + std::size_t syndrome_end_next_window = + layout.round_start((w + 1) * step_size + 1) - 1; const auto &h_ptr = this->H.ptr(); const auto &h_indices = this->H.indices(); - for (std::size_t s = 0; s < this->rolling_window.size(); ++s) { + for (std::size_t s = 0; s < batch_size; ++s) { for (std::size_t c = 0; c < num_to_commit; ++c) { if (rw_results[s].result[c + this_window_first_column]) { // Flip next-round syndrome bits where PCM has a 1 in this column. @@ -578,7 +344,7 @@ void sliding_window::decode_window() { auto this_window_first_column = first_columns[w]; auto num_to_commit = window_results[0].size(); CUDA_QEC_DBG(" Committing {} bits from window {}", num_to_commit, w); - for (std::size_t s = 0; s < this->rolling_window.size(); ++s) { + for (std::size_t s = 0; s < batch_size; ++s) { for (std::size_t c = 0; c < num_to_commit; ++c) { rw_results[s].result[c + this_window_first_column] = window_results[s][c]; diff --git a/libs/qec/lib/decoders/sliding_window.h b/libs/qec/lib/decoders/sliding_window.h index 7e7a23709..8b3ee0b10 100644 --- a/libs/qec/lib/decoders/sliding_window.h +++ b/libs/qec/lib/decoders/sliding_window.h @@ -9,6 +9,7 @@ #pragma once #include "cudaq/qec/decoder.h" +#include "cudaq/qec/pcm_utils.h" #include namespace cudaq::qec { @@ -29,7 +30,9 @@ class sliding_window : public decoder { std::size_t step_size = 1; /// The number of syndromes per round (the interior-layer width). std::size_t num_syndromes_per_round = 0; - /// The number of boundary measurements. + /// The width of the first/last (boundary) rounds. A caller-supplied 0 means + /// "uniform layout"; the constructor normalizes it to num_syndromes_per_round + /// so the rest of the class always sees a concrete boundary width. std::size_t num_boundary_syndromes = 0; /// When forming a window, should error mechanisms that span the start round /// and any preceding rounds be included? @@ -46,14 +49,15 @@ class sliding_window : public decoder { // Derived parameters. std::size_t num_windows = 0; - std::size_t num_rounds = 0; - std::size_t num_syndromes_per_window = 0; - /// Number of detector rows in an unpadded (boundary-layout) syndrome. - std::size_t unpadded_syndrome_size = 0; - std::size_t num_rounds_since_last_decode = 0; + /// Detector-row layers in the [B | S...S | B] layout (minimum 2 for a memory + /// circuit with no interior rounds) + std::size_t num_detector_layers = 0; std::vector> inner_decoders; std::vector first_columns; + // Boundary-aware detector-layer layout ([B | S...S | B]); maps between layers + // and detector rows. Shared with the get_pcm_for_rounds machinery. + details::round_layout layout; // Enum type for timing data. enum WindowProcTimes { INITIALIZE_WINDOW, // 0 @@ -68,12 +72,10 @@ class sliding_window : public decoder { }; // State data - std::vector> - rolling_window; // [batch_size, num_syndromes_per_window] - // rolling window read and write indices (circular buffer) - std::size_t rw_next_write_index = 0; // [0, num_syndromes_per_window) - std::size_t rw_next_read_index = 0; // [0, num_syndromes_per_window) - std::size_t rw_filled = 0; + std::vector>> window_rounds; + std::size_t window_start_round = 0; + std::size_t rounds_since_last_reset = 0; + std::size_t batch_size = 1; std::size_t num_windows_decoded = 0; std::vector> syndrome_mods; // [batch_size, syndrome_size] std::vector rw_results; // [batch_size] @@ -84,62 +86,15 @@ class sliding_window : public decoder { /// @brief Validate constructor inputs void validate_inputs(); - /// @brief Initialize the window - /// @param num_syndromes The number of syndromes to initialize the window for - void initialize_window(std::size_t num_syndromes); + /// @brief Reset per-syndrome streaming state and size the result/mod buffers. + /// @param batch_size The number of independent syndromes to decode together + /// (1 for non-batched mode). + void initialize_window(std::size_t batch_size); - /// @brief Add a single syndrome to the rolling window (circular buffer) - void add_syndrome_to_rolling_window(const std::vector &syndrome, - std::size_t syndrome_index, - bool update_next_write_index = true); - - /// @brief Add a batch of syndromes to the rolling window (circular buffer) - void add_syndromes_to_rolling_window( - const std::vector> &syndromes); - - /// @brief Get a single syndrome from the rolling window (unwrapping circular - /// buffer) - std::vector - get_syndrome_from_rolling_window(std::size_t syndrome_index); - - /// @brief Get a batch of syndromes from the rolling window (unwrapping - /// circular buffer) - std::vector> get_syndromes_from_rolling_window(); - - /// @brief Update the read index for the rolling window - void update_rw_next_read_index(); - - /// @brief Decode a single window (internal helper) + /// @brief Decode the active window from the rolling buffer, commit, and back + /// out committed errors into the next window's syndrome mods. void decode_window(); - /// @brief Whether the decoder is zero-padding narrower boundary layers up to - /// the interior width (i.e. a boundary-layout was supplied). - bool boundary_padding_active() const { - return num_boundary_syndromes != 0 && - num_boundary_syndromes != num_syndromes_per_round; - } - - /// @brief Pad an all-rounds syndrome from the unpadded boundary - /// layout up to the uniform padded layout by inserting zeros for the - /// boundary-layer padding rows. - std::vector pad_syndrome(const std::vector &syndrome) const; - - /// @brief Pad a single boundary detector layer up to the interior width by - /// appending zeros. - std::vector pad_round(const std::vector &round) const; - - /// @brief Zero-pad the two boundary layers of @p H up to the interior width - /// so all rounds are uniform. The first/last @p num_boundary_syndromes - /// rows are the narrower boundary layers and interior layers have - /// @p num_syndromes_per_round rows. Returns the canonical CSC form expected - /// by the base decoder. If no boundary padding is needed (@p - /// num_boundary_syndromes is 0 or equals @p num_syndromes_per_round), this - /// is just H.canonicalize().to_csc(). - static sparse_binary_matrix - prepare_pcm(const cudaq::qec::sparse_binary_matrix &H, - std::size_t num_syndromes_per_round, - std::size_t num_boundary_syndromes); - public: /// @brief Constructor /// @param H The full parity check matrix for all rounds diff --git a/libs/qec/lib/pcm_utils.cpp b/libs/qec/lib/pcm_utils.cpp index db244ef77..c07b0b9f9 100644 --- a/libs/qec/lib/pcm_utils.cpp +++ b/libs/qec/lib/pcm_utils.cpp @@ -118,6 +118,7 @@ void validate_reorder_pcm_columns_column_order( /// Validate the round-related parameters shared by both get_pcm_for_rounds /// overloads. Throws std::invalid_argument on violation. void validate_pcm_for_rounds_params(std::uint32_t num_syndromes_per_round, + std::uint32_t num_boundary_syndromes, std::uint32_t start_round, std::uint32_t end_round, std::size_t num_rows) { @@ -130,10 +131,23 @@ void validate_pcm_for_rounds_params(std::uint32_t num_syndromes_per_round, "get_pcm_for_rounds: num_syndromes_per_round must be less than the " "number of rows in PCM"); } - const std::size_t first_row_to_keep = - static_cast(start_round) * num_syndromes_per_round; - const std::size_t last_row_to_keep = - static_cast(end_round + 1) * num_syndromes_per_round - 1; + if (num_boundary_syndromes > num_syndromes_per_round) { + throw std::invalid_argument( + "get_pcm_for_rounds: num_boundary_syndromes must be <= " + "num_syndromes_per_round"); + } + const cudaq::qec::details::round_layout layout( + num_syndromes_per_round, num_boundary_syndromes, num_rows); + if (layout.boundary && + (num_rows < 2 * num_boundary_syndromes || + (num_rows - 2 * num_boundary_syndromes) % num_syndromes_per_round != + 0)) { + throw std::invalid_argument( + "get_pcm_for_rounds: number of PCM rows is inconsistent with the given " + "num_syndromes_per_round and num_boundary_syndromes"); + } + const std::size_t first_row_to_keep = layout.round_start(start_round); + const std::size_t last_row_to_keep = layout.round_start(end_round + 1) - 1; if (first_row_to_keep >= num_rows) { throw std::invalid_argument( "get_pcm_for_rounds: first_row_to_keep is greater than the number of " @@ -151,7 +165,7 @@ void validate_pcm_for_rounds_params(std::uint32_t num_syndromes_per_round, /// `get_pcm_for_rounds`). void select_pcm_columns_for_round_range( const std::vector> &row_indices, - std::uint32_t num_syndromes_per_round, std::uint32_t start_round, + const cudaq::qec::details::round_layout &layout, std::uint32_t start_round, std::uint32_t end_round, bool straddle_start_round, bool straddle_end_round, std::vector &columns_in_range_out, std::uint32_t &first_column, std::uint32_t &last_column) { @@ -160,8 +174,8 @@ void select_pcm_columns_for_round_range( const auto &col = row_indices[c]; if (col.empty()) continue; - auto first_round = col.front() / num_syndromes_per_round; - auto last_round = col.back() / num_syndromes_per_round; + auto first_round = layout.row_to_round(col.front()); + auto last_round = layout.row_to_round(col.back()); if (last_round < start_round || first_round > end_round) continue; @@ -338,6 +352,18 @@ bool pcm_is_sorted(const std::vector> &sparse_pcm, return true; } +bool pcm_is_sorted(const std::vector> &sparse_pcm, + std::uint32_t num_syndromes_per_round, + std::uint32_t num_boundary_syndromes) { + auto column_indices = get_sorted_pcm_column_indices( + sparse_pcm, num_syndromes_per_round, num_boundary_syndromes); + auto num_cols = sparse_pcm.size(); + for (std::size_t c = 0; c < num_cols; c++) + if (column_indices[c] != c) + return false; + return true; +} + /// @brief Check if a PCM is sorted. /// @param pcm The PCM to check. /// @param num_syndromes_per_round The number of syndromes per round. @@ -550,20 +576,23 @@ std::tuple, std::uint32_t, std::uint32_t> get_pcm_for_rounds(const cudaqx::tensor &pcm, std::uint32_t num_syndromes_per_round, std::uint32_t start_round, std::uint32_t end_round, - bool straddle_start_round, bool straddle_end_round) { - validate_pcm_for_rounds_params(num_syndromes_per_round, start_round, - end_round, pcm.shape()[0]); - auto first_row_to_keep = start_round * num_syndromes_per_round; - auto last_row_to_keep = (end_round + 1) * num_syndromes_per_round - 1; + bool straddle_start_round, bool straddle_end_round, + std::uint32_t num_boundary_syndromes) { + validate_pcm_for_rounds_params(num_syndromes_per_round, + num_boundary_syndromes, start_round, end_round, + pcm.shape()[0]); + const details::round_layout layout(num_syndromes_per_round, + num_boundary_syndromes, pcm.shape()[0]); + auto first_row_to_keep = layout.round_start(start_round); + auto last_row_to_keep = layout.round_start(end_round + 1) - 1; // Get a sparse representation of the PCM. auto row_indices = dense_to_sparse(pcm); std::vector columns_in_range; std::uint32_t first_column = 0, last_column = 0; select_pcm_columns_for_round_range( - row_indices, num_syndromes_per_round, start_round, end_round, - straddle_start_round, straddle_end_round, columns_in_range, first_column, - last_column); + row_indices, layout, start_round, end_round, straddle_start_round, + straddle_end_round, columns_in_range, first_column, last_column); return std::make_tuple(reorder_pcm_columns(pcm, columns_in_range, first_row_to_keep, @@ -576,13 +605,17 @@ 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, bool straddle_end_round, - bool pcm_is_canonical) { - validate_pcm_for_rounds_params(num_syndromes_per_round, start_round, - end_round, pcm.num_rows()); + bool pcm_is_canonical, + std::uint32_t num_boundary_syndromes) { + validate_pcm_for_rounds_params(num_syndromes_per_round, + num_boundary_syndromes, start_round, end_round, + pcm.num_rows()); + const details::round_layout layout(num_syndromes_per_round, + num_boundary_syndromes, pcm.num_rows()); auto first_row_to_keep = - static_cast(start_round * num_syndromes_per_round); + static_cast(layout.round_start(start_round)); auto last_row_to_keep = - static_cast((end_round + 1) * num_syndromes_per_round - 1); + static_cast(layout.round_start(end_round + 1) - 1); // select_pcm_columns_for_round_range reads .front()/.back() as min/max row, // which requires sorted-per-column input. Canonicalize unless the caller @@ -592,9 +625,8 @@ get_pcm_for_rounds(const sparse_binary_matrix &pcm, std::vector columns_in_range; std::uint32_t first_column = 0, last_column = 0; select_pcm_columns_for_round_range( - row_indices, num_syndromes_per_round, start_round, end_round, - straddle_start_round, straddle_end_round, columns_in_range, first_column, - last_column); + row_indices, layout, start_round, end_round, straddle_start_round, + straddle_end_round, columns_in_range, first_column, last_column); const std::uint32_t num_rows_to_copy = last_row_to_keep - first_row_to_keep + 1; diff --git a/libs/qec/unittests/test_decoders.cpp b/libs/qec/unittests/test_decoders.cpp index 890257a05..f98a5fb27 100644 --- a/libs/qec/unittests/test_decoders.cpp +++ b/libs/qec/unittests/test_decoders.cpp @@ -1192,8 +1192,8 @@ TEST(SlidingWindowDecoder, BaseStreamingCopiesFirstRoundDetectors) { } TEST(SlidingWindowDecoder, PreparePcmRejectsBadBoundaryLayout) { - // prepare_pcm runs in the constructor's member initializer list, so an - // inconsistent boundary layout must throw during construction. + // The boundary layout is validated during construction, so an inconsistent + // boundary layout must throw. auto params = [](std::size_t S, std::size_t B) { cudaqx::heterogeneous_map p; p.insert("window_size", std::size_t{1}); diff --git a/libs/qec/unittests/test_qec.cpp b/libs/qec/unittests/test_qec.cpp index 3284d234f..ccb72f048 100644 --- a/libs/qec/unittests/test_qec.cpp +++ b/libs/qec/unittests/test_qec.cpp @@ -1735,6 +1735,62 @@ TEST(PCMUtilsTester, checkGetSortedColumnIndicesWithBoundary) { std::invalid_argument); } +TEST(PCMUtilsTester, checkPCMIsSortedWithBoundary) { + // Same boundary layout as above (S = 8, B = 6). Rows < B are round 0; the + // boundary and uniform mappings disagree on rows [B, S). + const std::uint32_t S = 8, B = 6; + + // Columns in boundary-sorted order: spans (0,0) < (0,1) < (1,1). + std::vector> sorted = {{0}, {5, 12}, {6}}; + EXPECT_TRUE(cudaq::qec::pcm_is_sorted(sorted, S, B)); + + // Swapping the last two breaks the boundary order: {6} is round 1 and must + // follow {5, 12} (which ends in round 1), so this is not boundary-sorted. + std::vector> unsorted = {{0}, {6}, {5, 12}}; + EXPECT_FALSE(cudaq::qec::pcm_is_sorted(unsorted, S, B)); + + // The uniform overload disagrees: both {6} and {5, 12} + // start in round 0, so `unsorted` *is* uniform-sorted. + EXPECT_TRUE(cudaq::qec::pcm_is_sorted(unsorted, S)); +} + +TEST(PCMUtilsTester, checkGetPCMForRoundsWithBoundary) { + // Boundary layout [B | S | S | B] with B = 2, S = 4 -> 12 rows, 4 rounds: + // round 0: rows [0,2), round 1: [2,6), round 2: [6,10), round 3: [10,12). + const std::uint32_t S = 4, B = 2; + const std::size_t n_rows = 12; + + // Identity PCM: column j triggers exactly detector row j. + cudaqx::tensor pcm({n_rows, n_rows}); + for (std::size_t j = 0; j < n_rows; ++j) + pcm.at({j, j}) = 1; + + // The leading boundary round is B rows wide, not S: the boundary overload + // keeps rows [0,B) while the uniform interpretation keeps rows [0,S). + auto [sub_b, fc_b, lc_b] = cudaq::qec::get_pcm_for_rounds( + pcm, S, /*start_round=*/0, /*end_round=*/0, /*straddle_start=*/false, + /*straddle_end=*/false, /*num_boundary_syndromes=*/B); + EXPECT_EQ(sub_b.shape()[0], B); + auto [sub_u, fc_u, lc_u] = cudaq::qec::get_pcm_for_rounds( + pcm, S, /*start_round=*/0, /*end_round=*/0); + EXPECT_EQ(sub_u.shape()[0], S); + + // Spanning every round returns all rows (confirms the round count is right). + auto [sub_all, fc_all, lc_all] = cudaq::qec::get_pcm_for_rounds( + pcm, S, /*start_round=*/0, /*end_round=*/3, false, false, B); + EXPECT_EQ(sub_all.shape()[0], n_rows); + + // Invalid boundary arguments / inconsistent row counts must throw. + EXPECT_THROW( + cudaq::qec::get_pcm_for_rounds(pcm, S, 0, 0, false, false, S + 1), + std::invalid_argument); + // 13 rows cannot be split as 2*B + K*S for B = 2, S = 4. + cudaqx::tensor bad_pcm({13, 13}); + EXPECT_THROW( + cudaq::qec::get_pcm_for_rounds(bad_pcm, S, 0, 0, false, false, B), + std::invalid_argument); +} + TEST(PCMUtilsTester, checkShufflePCMColumns) { std::size_t n_rounds = 4; std::size_t n_errs_per_round = 30;