diff --git a/libs/qec/include/cudaq/qec/code.h b/libs/qec/include/cudaq/qec/code.h index f7b9ea1dc..5b18182d1 100644 --- a/libs/qec/include/cudaq/qec/code.h +++ b/libs/qec/include/cudaq/qec/code.h @@ -201,6 +201,50 @@ class code : public cudaqx::extension_point { /// @return Tensor representing Lz cudaqx::tensor get_observables_z() const; + /// @brief Get the inlined feedback matrix for detector construction + /// + /// Shape is [numCols x numCols], where numCols = + /// get_num_ancilla_qubits() (one measurement record per ancilla per + /// round, in [Z][X] order). Entries must be exactly 0 or 1. Entry + /// (j, k) = 1 means the cross-round + /// detector comparing record j between consecutive rounds additionally + /// XORs record k of the earlier round, and the final boundary detector + /// for record j additionally XORs record k of the last round. + /// @return Tensor representing the inlined feedback matrix. An empty + /// tensor (the default) means no feedback is applied and detectors are + /// formed from same-record comparisons only. + virtual cudaqx::tensor get_inlined_feedback() const; + + /// @brief Get the inlined feedback matrix for logical observables measured + /// in the Z basis + /// + /// Shape is [num_observables x numCols], where numCols = + /// get_num_ancilla_qubits() (one measurement record per ancilla per + /// round, in [Z][X] order). Entries must be exactly 0 or 1. Entry + /// (m, k) = 1 means logical observable m + /// additionally XORs record k of every round. This getter is consumed when + /// the memory experiment measures its observables in the Z basis (prep0 / + /// prep1 state preparations). + /// @return Tensor representing the Z-basis observable inlined feedback + /// matrix. An empty tensor (the default) means no observable feedback is + /// applied in the Z basis. + virtual cudaqx::tensor get_observable_inlined_feedback_z() const; + + /// @brief Get the inlined feedback matrix for logical observables measured + /// in the X basis + /// + /// Shape is [num_observables x numCols], where numCols = + /// get_num_ancilla_qubits() (one measurement record per ancilla per + /// round, in [Z][X] order). Entries must be exactly 0 or 1. Entry + /// (m, k) = 1 means logical observable m + /// additionally XORs record k of every round. This getter is consumed when + /// the memory experiment measures its observables in the X basis (prepp / + /// prepm state preparations). + /// @return Tensor representing the X-basis observable inlined feedback + /// matrix. An empty tensor (the default) means no observable feedback is + /// applied in the X basis. + virtual cudaqx::tensor get_observable_inlined_feedback_x() const; + /// @brief Get the stabilizer generators /// @return Reference to stabilizers const std::vector &get_stabilizers() const { diff --git a/libs/qec/include/cudaq/qec/detector_error_model.h b/libs/qec/include/cudaq/qec/detector_error_model.h index 5f3e870ef..6ab74e207 100644 --- a/libs/qec/include/cudaq/qec/detector_error_model.h +++ b/libs/qec/include/cudaq/qec/detector_error_model.h @@ -1,5 +1,5 @@ /****************************************************************-*- C++ -*-**** - * Copyright (c) 2025 NVIDIA Corporation & Affiliates. * + * Copyright (c) 2025 - 2026 NVIDIA Corporation & Affiliates. * * All rights reserved. * * * * This source code and the accompanying materials are made available under * diff --git a/libs/qec/include/cudaq/qec/experiments.h b/libs/qec/include/cudaq/qec/experiments.h index e4ebe184b..0eeb702b9 100644 --- a/libs/qec/include/cudaq/qec/experiments.h +++ b/libs/qec/include/cudaq/qec/experiments.h @@ -93,7 +93,10 @@ dem_sampling(const cudaqx::tensor &check_matrix, std::size_t numShots, /// type matching `statePrep`'s basis, since only that type is deterministic /// at the circuit's endpoints), then one detector block per each of the /// `numRounds - 1` inter-round transitions, then `numFixed` more boundary -/// detectors. `dataResults` has shape `(numShots, numDataQubits)`. +/// detectors. `dataResults` has shape `(numShots, numDataQubits)` and contains +/// the raw final data-qubit measurements. In particular, observable inlined +/// feedback is represented in a generated detector error model, but is not +/// folded into `dataResults` as a corrected logical-observable value. std::tuple, cudaqx::tensor> sample_memory_circuit(const code &code, operation statePrep, std::size_t numShots, std::size_t numRounds, diff --git a/libs/qec/lib/code.cpp b/libs/qec/lib/code.cpp index 19693b8c5..39be5efed 100644 --- a/libs/qec/lib/code.cpp +++ b/libs/qec/lib/code.cpp @@ -61,6 +61,18 @@ cudaqx::tensor code::get_observables_z() const { return to_parity_matrix(m_pauli_observables, stabilizer_type::Z); } +cudaqx::tensor code::get_inlined_feedback() const { + return cudaqx::tensor(); +} + +cudaqx::tensor code::get_observable_inlined_feedback_z() const { + return cudaqx::tensor(); +} + +cudaqx::tensor code::get_observable_inlined_feedback_x() const { + return cudaqx::tensor(); +} + std::unique_ptr get_code(const std::string &name, const std::vector &stab, const heterogeneous_map options) { diff --git a/libs/qec/lib/device/inlined_feedback.h b/libs/qec/lib/device/inlined_feedback.h new file mode 100644 index 000000000..702202f25 --- /dev/null +++ b/libs/qec/lib/device/inlined_feedback.h @@ -0,0 +1,161 @@ +/****************************************************************-*- C++ -*-**** + * Copyright (c) 2024 - 2025 NVIDIA Corporation & Affiliates. * + * All rights reserved. * + * * + * This source code and the accompanying materials are made available under * + * the terms of the Apache License 2.0 which accompanies this distribution. * + ******************************************************************************/ + +#pragma once + +#include "cudaq.h" + +namespace cudaq::qec { + +/// @brief Number of nonzero entries in row @p row of a row-major 0/1 matrix +/// that has @p num_cols columns (i.e. the weight of one feedback/observable/ +/// stabilizer row). +inline __qpu__ std::size_t +row_support_weight(const std::vector &matrix_flat, std::size_t row, + std::size_t num_cols) { + std::size_t weight = 0; + for (std::size_t k = 0; k < num_cols; ++k) + if (matrix_flat[row * num_cols + k] != 0) + weight++; + return weight; +} + +/// @brief Records for the cross-round detector of syndrome record @p j: the +/// earlier-vs-current comparison `{prev[j], curr[j]}` followed by the +/// earlier-round herald records selected by the CSR row for j. +inline __qpu__ std::vector +cross_round_detector_records(const std::vector &prev, + const std::vector &curr, + std::size_t j, + const std::vector &feedback_indices, + const std::vector &feedback_offsets) { + std::size_t weight = feedback_offsets[j + 1] - feedback_offsets[j]; + std::vector det(2 + weight); + det[0] = prev[j]; + det[1] = curr[j]; + std::size_t idx = 2; + for (std::size_t i = feedback_offsets[j]; i < feedback_offsets[j + 1]; ++i) + det[idx++] = prev[feedback_indices[i]]; + return det; +} + +/// @brief Records for the boundary detector of syndrome record @p record_row: +/// `{last_syndrome[record_row]}`, then the data-qubit readouts +/// `data_results[q]` for each data qubit q (ascending) in the stabilizer +/// support `stabilizers[row_base + q]`, then the last-round herald records +/// selected by the CSR row for @p record_row. Empty feedback offsets mean the +/// herald part contributes nothing, so one helper serves both feedback and +/// legacy branches. +inline __qpu__ std::vector boundary_detector_records( + const std::vector &last_syndrome, + std::size_t record_row, + const std::vector &data_results, + const std::vector &stabilizers, std::size_t row_base, + std::size_t num_data, const std::vector &feedback_indices, + const std::vector &feedback_offsets) { + std::size_t support_weight = 0; + for (std::size_t q = 0; q < num_data; ++q) + if (stabilizers[row_base + q] != 0) + support_weight++; + std::size_t fb_weight = + feedback_offsets.size() > 0 + ? feedback_offsets[record_row + 1] - feedback_offsets[record_row] + : 0; + std::vector support(1 + support_weight + fb_weight); + support[0] = last_syndrome[record_row]; + std::size_t idx = 1; + for (std::size_t q = 0; q < num_data; ++q) + if (stabilizers[row_base + q] != 0) + support[idx++] = data_results[q]; + if (feedback_offsets.size() > 0) + for (std::size_t i = feedback_offsets[record_row]; + i < feedback_offsets[record_row + 1]; ++i) + support[idx++] = last_syndrome[feedback_indices[i]]; + return support; +} + +/// @brief Round-major slice offsets into the observable-feedback buffer. Entry +/// m is the start index of observable m's slice; each observable occupies +/// `num_rounds * (obs_feedback_offsets[m+1] - obs_feedback_offsets[m])` +/// records. The returned vector has `num_observables + 1` entries; the last +/// entry is the total buffer size. Returns an empty vector when the input +/// offsets are empty (no observable feedback declared). +inline __qpu__ std::vector observable_feedback_record_offsets( + const std::vector &obs_feedback_offsets, + std::size_t num_observables, std::size_t num_rounds) { + // Sized (not default) construction: the __qpu__ dialect forbids the default + // std::vector constructor. Empty CSR offsets yield a size-0 vector. + std::size_t num_offsets = + obs_feedback_offsets.size() > 0 ? num_observables + 1 : 0; + std::vector offsets(num_offsets); + if (obs_feedback_offsets.size() == 0) + return offsets; + std::size_t total = 0; + for (std::size_t m = 0; m < num_observables; ++m) { + offsets[m] = total; + total += + num_rounds * (obs_feedback_offsets[m + 1] - obs_feedback_offsets[m]); + } + offsets[num_observables] = total; + return offsets; +} + +/// @brief Write the records that observable feedback collects during @p round +/// into the round-major @p obs_fb_records buffer. For each observable m, the +/// round-r block starts at `record_offsets[m] + round * row_weight(m)` and +/// holds the syndrome records selected by CSR row m. Call with round 0 for the +/// first round's syndrome and with the running round index for each subsequent +/// round. +inline __qpu__ void collect_observable_feedback_round( + std::vector &obs_fb_records, + const std::vector &record_offsets, + const std::vector &syndrome, std::size_t round, + const std::vector &obs_feedback_indices, + const std::vector &obs_feedback_offsets, + std::size_t num_observables) { + for (std::size_t m = 0; m < num_observables; ++m) { + std::size_t weight = obs_feedback_offsets[m + 1] - obs_feedback_offsets[m]; + std::size_t idx = record_offsets[m] + round * weight; + for (std::size_t i = obs_feedback_offsets[m]; + i < obs_feedback_offsets[m + 1]; ++i) + obs_fb_records[idx++] = syndrome[obs_feedback_indices[i]]; + } +} + +/// @brief Support records for logical observable @p obs: the feedback records +/// collected in @p obs_fb_records (round-major, occupying +/// `[offsets[obs], offsets[obs + 1])`), followed by the data-qubit readouts +/// `data_results[q]` for each data qubit q (ascending) with +/// `obs_matrix_flat(obs, q) != 0`. When @p offsets is empty (no observable +/// feedback declared) the feedback prefix is empty and only the data support +/// is returned. +inline __qpu__ std::vector observable_support_records( + std::size_t obs, const std::vector &obs_matrix_flat, + std::size_t num_data, + const std::vector &data_results, + const std::vector &obs_fb_records, + const std::vector &offsets) { + std::size_t support_weight = + row_support_weight(obs_matrix_flat, obs, num_data); + std::size_t fb_count = 0; + std::size_t base = 0; + if (offsets.size() > 0) { + base = offsets[obs]; + fb_count = offsets[obs + 1] - offsets[obs]; + } + std::vector obs_support(fb_count + support_weight); + for (std::size_t i = 0; i < fb_count; ++i) + obs_support[i] = obs_fb_records[base + i]; + std::size_t idx = fb_count; + for (std::size_t q = 0; q < num_data; ++q) + if (obs_matrix_flat[obs * num_data + q] != 0) + obs_support[idx++] = data_results[q]; + return obs_support; +} + +} // namespace cudaq::qec diff --git a/libs/qec/lib/device/memory_circuit.cpp b/libs/qec/lib/device/memory_circuit.cpp index b00754a28..b79a63a87 100644 --- a/libs/qec/lib/device/memory_circuit.cpp +++ b/libs/qec/lib/device/memory_circuit.cpp @@ -6,19 +6,22 @@ * the terms of the Apache License 2.0 which accompanies this distribution. * ******************************************************************************/ #include "memory_circuit.h" -#include +#include "inlined_feedback.h" namespace cudaq::qec { -__qpu__ void memory_circuit(const code::stabilizer_round &stabilizer_round, - const code::one_qubit_encoding &statePrep, - std::size_t num_data, std::size_t numAncx, - std::size_t numAncz, std::size_t num_rounds, - const std::vector &x_stabilizers, - const std::vector &z_stabilizers, - const std::vector &obs_matrix_flat, - std::size_t num_observables, - bool measure_in_x_basis) { +__qpu__ void +memory_circuit(const code::stabilizer_round &stabilizer_round, + const code::one_qubit_encoding &statePrep, std::size_t num_data, + std::size_t numAncx, std::size_t numAncz, std::size_t num_rounds, + const std::vector &x_stabilizers, + const std::vector &z_stabilizers, + const std::vector &obs_matrix_flat, + std::size_t num_observables, bool measure_in_x_basis, + const std::vector &feedback_indices, + const std::vector &feedback_offsets, + const std::vector &obs_feedback_indices, + const std::vector &obs_feedback_offsets) { // Allocate the data and ancilla qubits cudaq::qvector data(num_data), xstab_anc(numAncx), zstab_anc(numAncz); @@ -39,10 +42,44 @@ __qpu__ void memory_circuit(const code::stabilizer_round &stabilizer_round, cudaq::detector(final_syndrome[fixed_offset + i]); } + std::size_t numCols = numAncx + numAncz; + + // Observable inlined feedback accumulation. Nested std::vector is not + // supported in __qpu__ code, so instead of one record vector per observable + // we use a single flat buffer with per-observable round-major slices; see + // observable_feedback_record_offsets / collect_observable_feedback_round in + // inlined_feedback.h. + std::vector obs_fb_record_offsets = + observable_feedback_record_offsets(obs_feedback_offsets, num_observables, + num_rounds); + std::size_t obs_fb_total = obs_feedback_offsets.size() > 0 + ? obs_fb_record_offsets[num_observables] + : 0; + std::vector obs_fb_records(obs_fb_total); + + // Collect the first-round feedback records for each observable. + if (obs_feedback_offsets.size() > 0) + collect_observable_feedback_round(obs_fb_records, obs_fb_record_offsets, + final_syndrome, 0, obs_feedback_indices, + obs_feedback_offsets, num_observables); + // Generate syndrome data for (std::size_t round = 1; round < num_rounds; ++round) { auto syndrome = stabilizer_round(logical, x_stabilizers, z_stabilizers); - cudaq::detectors(final_syndrome, syndrome); + if (obs_feedback_offsets.size() > 0) + collect_observable_feedback_round(obs_fb_records, obs_fb_record_offsets, + syndrome, round, obs_feedback_indices, + obs_feedback_offsets, num_observables); + if (feedback_offsets.size() == 0) { + cudaq::detectors(final_syndrome, syndrome); + } else { + // Cross-round detector for record j: earlier vs current round record, + // augmented with the earlier-round records declared in row j of the + // feedback matrix. + for (std::size_t j = 0; j < numCols; ++j) + cudaq::detector(cross_round_detector_records( + final_syndrome, syndrome, j, feedback_indices, feedback_offsets)); + } final_syndrome = syndrome; } @@ -51,48 +88,23 @@ __qpu__ void memory_circuit(const code::stabilizer_round &stabilizer_round, } auto data_results = mz(data); - // Emit one logical_observable per row of the observable matrix. - for (std::size_t obs = 0; obs < num_observables; ++obs) { - std::size_t support_weight = 0; - for (std::size_t q = 0; q < num_data; ++q) { - if (obs_matrix_flat[obs * num_data + q] != 0) - support_weight++; - } - std::vector obs_support(support_weight); - std::size_t idx = 0; - for (std::size_t q = 0; q < num_data; ++q) { - if (obs_matrix_flat[obs * num_data + q] != 0) - obs_support[idx++] = data_results[q]; - } - cudaq::logical_observable(obs_support); - } + // Emit one logical_observable per row of the observable matrix: the per-round + // feedback records followed by the data-qubit support. + for (std::size_t obs = 0; obs < num_observables; ++obs) + cudaq::logical_observable( + observable_support_records(obs, obs_matrix_flat, num_data, data_results, + obs_fb_records, obs_fb_record_offsets)); // For each stabilizer, form detectors from data qubit readout connected with - // final stabilizer round. + // final stabilizer round. With inlined feedback, the boundary detector for + // record (fixed_offset + x) is extended with the declared last-round records. const std::vector &stabilizers = measure_in_x_basis ? x_stabilizers : z_stabilizers; - for (std::size_t x = 0; x < num_fixed_measurements; ++x) { - std::size_t row_base = x * num_data; - - std::size_t support_weight = 0; - for (std::size_t q = 0; q < num_data; ++q) { - if (stabilizers[row_base + q] != 0) { - support_weight++; - } - } - - std::vector support(support_weight + 1); - support[0] = final_syndrome[fixed_offset + x]; - std::size_t support_idx = 1; - for (std::size_t q = 0; q < num_data; ++q) { - if (stabilizers[row_base + q] != 0) { - support[support_idx++] = data_results[q]; - } - } - - cudaq::detector(support); - } + for (std::size_t x = 0; x < num_fixed_measurements; ++x) + cudaq::detector(boundary_detector_records( + final_syndrome, fixed_offset + x, data_results, stabilizers, + x * num_data, num_data, feedback_indices, feedback_offsets)); } } // namespace cudaq::qec diff --git a/libs/qec/lib/device/memory_circuit.h b/libs/qec/lib/device/memory_circuit.h index 1f0280d9a..7c777cca5 100644 --- a/libs/qec/lib/device/memory_circuit.h +++ b/libs/qec/lib/device/memory_circuit.h @@ -29,13 +29,23 @@ namespace cudaq::qec { /// (num_observables × numData entries, values 0/1). /// @param num_observables Number of rows in the observable matrix (k). /// @param measure_in_x_basis Performing X- or Z-memory circuit -__qpu__ void memory_circuit(const code::stabilizer_round &stabilizer_round, - const code::one_qubit_encoding &statePrep, - std::size_t numData, std::size_t numAncx, - std::size_t numAncz, std::size_t numRounds, - const std::vector &x_stabilizers, - const std::vector &z_stabilizers, - const std::vector &obs_matrix_flat, - std::size_t num_observables, - bool measure_in_x_basis); +/// @param feedback_indices CSR column indices for detector feedback. +/// @param feedback_offsets CSR row offsets for detector feedback. Size is 0 +/// (no feedback; legacy detector structure) or numCols + 1 with +/// numCols = numAncx + numAncz. +/// @param obs_feedback_indices CSR column indices for observable feedback. +/// @param obs_feedback_offsets CSR row offsets for observable feedback. Size +/// is 0 (no feedback) or num_observables + 1. +__qpu__ void +memory_circuit(const code::stabilizer_round &stabilizer_round, + const code::one_qubit_encoding &statePrep, std::size_t numData, + std::size_t numAncx, std::size_t numAncz, std::size_t numRounds, + const std::vector &x_stabilizers, + const std::vector &z_stabilizers, + const std::vector &obs_matrix_flat, + std::size_t num_observables, bool measure_in_x_basis, + const std::vector &feedback_indices, + const std::vector &feedback_offsets, + const std::vector &obs_feedback_indices, + const std::vector &obs_feedback_offsets); } // namespace cudaq::qec diff --git a/libs/qec/lib/experiments.cpp b/libs/qec/lib/experiments.cpp index 4f1c28e9c..4cf1a1219 100644 --- a/libs/qec/lib/experiments.cpp +++ b/libs/qec/lib/experiments.cpp @@ -8,6 +8,7 @@ #include "cudaq/qec/experiments.h" #include "device/memory_circuit.h" +#include "inlined_feedback_layout.h" #include "cudaq/algorithms/dem.h" #include "cudaq/qec/dem_sampling.h" #include "cudaq/qec/pcm_utils.h" @@ -309,21 +310,34 @@ sample_memory_circuit(const code &code, operation statePrep, const std::size_t numCols = numAncx + numAncz; + auto feedback_layout = build_inlined_feedback_layout( + code.get_inlined_feedback(), + is_z_prep ? code.get_observable_inlined_feedback_z() + : code.get_observable_inlined_feedback_x(), + numCols, num_obs, "get_inlined_feedback()", + is_z_prep ? "get_observable_inlined_feedback_z()" + : "get_observable_inlined_feedback_x()"); + // Obtain the Measurement-to-Detector (M2D) sparse matrix. // m2d.rows[d] = set of chronological measurement indices whose XOR = detector // d. cudaq::M2DSparseMatrix m2d; cudaq::M2OSparseMatrix m2o; - cudaq::dem_from_kernel(memory_circuit, &noise, /*options=*/{}, m2d, m2o, - stabRound, prep, numData, numAncx, numAncz, numRounds, - xVec, zVec, obs_flat, num_obs, !is_z_prep); + cudaq::dem_from_kernel( + memory_circuit, &noise, /*options=*/{}, m2d, m2o, stabRound, prep, + numData, numAncx, numAncz, numRounds, xVec, zVec, obs_flat, num_obs, + !is_z_prep, feedback_layout.detector_indices, + feedback_layout.detector_offsets, feedback_layout.observable_indices, + feedback_layout.observable_offsets); // Sample the memory circuit and collect all raw measurements. cudaq::sample_options opts{ .shots = numShots, .noise = noise, .explicit_measurements = true}; - auto result = cudaq::sample(opts, memory_circuit, stabRound, prep, numData, - numAncx, numAncz, numRounds, xVec, zVec, obs_flat, - num_obs, !is_z_prep); + auto result = cudaq::sample( + opts, memory_circuit, stabRound, prep, numData, numAncx, numAncz, + numRounds, xVec, zVec, obs_flat, num_obs, !is_z_prep, + feedback_layout.detector_indices, feedback_layout.detector_offsets, + feedback_layout.observable_indices, feedback_layout.observable_offsets); // mzTable[shot, meas_idx] = raw 0/1 outcome; shape (numShots, // numMeasPerShot). Measurement layout per shot: numRounds*numCols ancilla, @@ -479,11 +493,22 @@ dem_from_memory_circuit(const code &code, operation statePrep, std::vector obs_flat(logical_obs.data(), logical_obs.data() + logical_obs.size()); + const std::size_t numAnc = numAncx + numAncz; + auto feedback_layout = build_inlined_feedback_layout( + code.get_inlined_feedback(), + is_z_prep ? code.get_observable_inlined_feedback_z() + : code.get_observable_inlined_feedback_x(), + numAnc, num_obs, "get_inlined_feedback()", + is_z_prep ? "get_observable_inlined_feedback_z()" + : "get_observable_inlined_feedback_x()"); + cudaq::dem_options dem_opts; dem_opts.decompose_errors = decompose_errors; auto dem_text = cudaq::dem_from_kernel( memory_circuit, &noise, dem_opts, stabRound, prep, numData, numAncx, - numAncz, numRounds, xVec, zVec, obs_flat, num_obs, !is_z_prep); + numAncz, numRounds, xVec, zVec, obs_flat, num_obs, !is_z_prep, + feedback_layout.detector_indices, feedback_layout.detector_offsets, + feedback_layout.observable_indices, feedback_layout.observable_offsets); auto dem = cudaq::qec::dem_from_stim_text(dem_text, decompose_errors); const auto numXStabs = code.get_num_x_stabilizers(); diff --git a/libs/qec/lib/inlined_feedback_layout.h b/libs/qec/lib/inlined_feedback_layout.h new file mode 100644 index 000000000..0afc2c28d --- /dev/null +++ b/libs/qec/lib/inlined_feedback_layout.h @@ -0,0 +1,90 @@ +/******************************************************************************* + * Copyright (c) 2026 NVIDIA Corporation & Affiliates. * + * All rights reserved. * + * * + * This source code and the accompanying materials are made available under * + * the terms of the Apache License 2.0 which accompanies this distribution. * + ******************************************************************************/ + +#pragma once + +#include "cuda-qx/core/tensor.h" +#include +#include +#include +#include +#include + +namespace cudaq::qec::details { + +struct inlined_feedback_layout { + std::vector detector_indices; + std::vector detector_offsets; + std::vector observable_indices; + std::vector observable_offsets; +}; + +inline void validate_inlined_feedback_tensor( + const cudaqx::tensor &tensor, std::size_t expected_rows, + std::size_t expected_cols, const std::string &name) { + if (tensor.size() == 0) + return; + + if (tensor.rank() != 2 || tensor.shape()[0] != expected_rows || + tensor.shape()[1] != expected_cols) { + std::string actual; + for (std::size_t d = 0; d < tensor.rank(); ++d) + actual += (d ? ", " : "") + std::to_string(tensor.shape()[d]); + throw std::runtime_error(name + " has invalid shape [" + actual + + "] - expected [" + std::to_string(expected_rows) + + ", " + std::to_string(expected_cols) + + "] or an empty tensor."); + } + + for (std::size_t i = 0; i < tensor.size(); ++i) { + const auto value = tensor.data()[i]; + if (value > 1) + throw std::runtime_error( + name + " has non-binary value " + + std::to_string(static_cast(value)) + " at flat index " + + std::to_string(i) + "; entries must be exactly 0 or 1."); + } +} + +inline void feedback_rows_to_csr(const cudaqx::tensor &matrix, + std::size_t num_rows, std::size_t num_cols, + const std::string &name, + std::vector &indices, + std::vector &offsets) { + validate_inlined_feedback_tensor(matrix, num_rows, num_cols, name); + if (matrix.size() == 0) + return; + + offsets.resize(num_rows + 1); + offsets[0] = 0; + for (std::size_t row = 0; row < num_rows; ++row) { + for (std::size_t col = 0; col < num_cols; ++col) + if (matrix.at({row, col}) != 0) + indices.push_back(col); + offsets[row + 1] = indices.size(); + } +} + +inline inlined_feedback_layout build_inlined_feedback_layout( + const cudaqx::tensor &feedback, + const cudaqx::tensor &observable_feedback, + std::size_t num_syndromes_per_round, std::size_t num_observables, + const std::string &feedback_name = "inlined feedback", + const std::string &observable_feedback_name = + "observable inlined feedback") { + inlined_feedback_layout layout; + feedback_rows_to_csr(feedback, num_syndromes_per_round, + num_syndromes_per_round, feedback_name, + layout.detector_indices, layout.detector_offsets); + feedback_rows_to_csr(observable_feedback, num_observables, + num_syndromes_per_round, observable_feedback_name, + layout.observable_indices, layout.observable_offsets); + return layout; +} + +} // namespace cudaq::qec::details diff --git a/libs/qec/python/bindings/py_code.cpp b/libs/qec/python/bindings/py_code.cpp index 2d0ba7e24..c10b1c17e 100644 --- a/libs/qec/python/bindings/py_code.cpp +++ b/libs/qec/python/bindings/py_code.cpp @@ -29,11 +29,70 @@ using namespace cudaqx; namespace cudaq::qec { +/// @brief Convert the array returned by a Python inlined-feedback getter to +/// an owning cudaqx::tensor. +/// @details Exact binary values from other dtypes or nested lists are accepted, +/// then copied into a C-contiguous uint8 numpy array. Validation happens before +/// conversion so fractional, negative, and overflowing values cannot silently +/// truncate or wrap. toTensor only borrows the numpy buffer, so the data is +/// copied into the returned tensor to outlive the temporary Python array. +static cudaqx::tensor feedbackArrayToTensor(nb::object array, + const char *methodName) { + auto np = nb::module_::import_("numpy"); + nb::object input = np.attr("asarray")(array); + if (nb::cast(input.attr("ndim")) != 2) + throw std::runtime_error(std::string(methodName) + + " must return a 2-D array."); + + nb::object isBinary = np.attr("logical_or")(np.attr("equal")(input, 0), + np.attr("equal")(input, 1)); + if (!nb::cast(np.attr("all")(isBinary).attr("item")())) + throw std::runtime_error(std::string(methodName) + + " must contain only exact 0 or 1 entries."); + + nb::object coerced = np.attr("ascontiguousarray")(input, np.attr("uint8")); + auto ndarray = nb::cast>(coerced); + auto borrowed = cudaqx::toTensor(ndarray); + cudaqx::tensor owned; + owned.copy(borrowed.data(), borrowed.shape()); + return owned; +} + class PyCode : public qec::code { public: - NB_TRAMPOLINE(qec::code, 6); + NB_TRAMPOLINE(qec::code, 9); protected: + // Trampoline methods for the optional inlined-feedback virtuals: forward to + // a Python override when one exists (converting its numpy result), else + // fall back to the base-class empty default. This is NB_OVERRIDE with a + // custom return conversion, since cudaqx::tensor has no type caster. + cudaqx::tensor get_inlined_feedback() const override { + nb::detail::ticket nb_ticket(nb_trampoline, "get_inlined_feedback", false); + if (nb_ticket.key.is_valid()) + return feedbackArrayToTensor(nb_trampoline.base().attr(nb_ticket.key)(), + "get_inlined_feedback"); + return NBBase::get_inlined_feedback(); + } + + cudaqx::tensor get_observable_inlined_feedback_z() const override { + nb::detail::ticket nb_ticket(nb_trampoline, + "get_observable_inlined_feedback_z", false); + if (nb_ticket.key.is_valid()) + return feedbackArrayToTensor(nb_trampoline.base().attr(nb_ticket.key)(), + "get_observable_inlined_feedback_z"); + return NBBase::get_observable_inlined_feedback_z(); + } + + cudaqx::tensor get_observable_inlined_feedback_x() const override { + nb::detail::ticket nb_ticket(nb_trampoline, + "get_observable_inlined_feedback_x", false); + if (nb_ticket.key.is_valid()) + return feedbackArrayToTensor(nb_trampoline.base().attr(nb_ticket.key)(), + "get_observable_inlined_feedback_x"); + return NBBase::get_observable_inlined_feedback_x(); + } + // Trampoline methods for pure virtual functions std::size_t get_num_data_qubits() const override { NB_OVERRIDE_PURE(get_num_data_qubits); @@ -186,6 +245,38 @@ class PyCodeHandle : public qec::code { return m_py_operation_encodings; } + bool declares_inlined_feedback() const { + return nb::hasattr(pyCode, "get_inlined_feedback") || + nb::hasattr(pyCode, "get_observable_inlined_feedback_z") || + nb::hasattr(pyCode, "get_observable_inlined_feedback_x"); + } + + /// @brief Forward the optional inlined-feedback declarations to the + /// registered Python code when it defines them; otherwise return the + /// base-class empty default (identity detector layout). + cudaqx::tensor get_inlined_feedback() const override { + if (!nb::hasattr(pyCode, "get_inlined_feedback")) + return code::get_inlined_feedback(); + return feedbackArrayToTensor(pyCode.attr("get_inlined_feedback")(), + "get_inlined_feedback"); + } + + cudaqx::tensor get_observable_inlined_feedback_z() const override { + if (!nb::hasattr(pyCode, "get_observable_inlined_feedback_z")) + return code::get_observable_inlined_feedback_z(); + return feedbackArrayToTensor( + pyCode.attr("get_observable_inlined_feedback_z")(), + "get_observable_inlined_feedback_z"); + } + + cudaqx::tensor get_observable_inlined_feedback_x() const override { + if (!nb::hasattr(pyCode, "get_observable_inlined_feedback_x")) + return code::get_observable_inlined_feedback_x(); + return feedbackArrayToTensor( + pyCode.attr("get_observable_inlined_feedback_x")(), + "get_observable_inlined_feedback_x"); + } + protected: // Trampoline methods for pure virtual functions std::size_t get_num_data_qubits() const override { @@ -242,6 +333,22 @@ static auto copyCodeMatrixToPyArray(const cudaqx::tensor &tensor, return cudaq::python::copyCUDAQXTensorToPyArray(tensor); } + +static void ensure_feedback_sampling_target_supported(const code &self) { + auto *pyHandle = dynamic_cast(&self); + if (!pyHandle || !pyHandle->declares_inlined_feedback()) + return; + + const auto targetName = nb::cast( + nb::module_::import_("cudaq").attr("get_target")().attr("name")); + if (targetName == "qpp-cpu") + throw std::runtime_error( + "Memory-circuit sampling for Python codes that declare inlined " + "feedback is not supported on target 'qpp-cpu': this target does not " + "preserve the required cross-round mid-circuit measurement " + "correlations. Use cudaq.set_target('stim'). Detector-error-model " + "generation remains supported."); +} } // namespace // Registry to store code factory functions @@ -548,6 +655,7 @@ void bindCode(nb::module_ &mod) { "sample_memory_circuit", [](code &code, std::size_t numShots, std::size_t numRounds, std::optional noise = std::nullopt) { + ensure_feedback_sampling_target_supported(code); auto [synd, dataRes] = noise ? sample_memory_circuit(code, numShots, numRounds, *noise) : sample_memory_circuit(code, numShots, numRounds); @@ -555,12 +663,17 @@ void bindCode(nb::module_ &mod) { cudaq::python::copyCUDAQXTensorToPyArray(synd), cudaq::python::copyCUDAQXTensorToPyArray(dataRes)); }, - "Sample the memory circuit of the code", nb::arg("code"), - nb::arg("numShots"), nb::arg("numRounds"), nb::arg("noise") = nb::none()); + "Sample the memory circuit of the code. Returns (syndromes, data), " + "where data contains raw final data-qubit measurements; observable " + "inlined feedback is represented in generated detector error models and " + "is not folded into data as a corrected logical-observable value.", + nb::arg("code"), nb::arg("numShots"), nb::arg("numRounds"), + nb::arg("noise") = nb::none()); qecmod.def( "sample_memory_circuit", [](code &code, operation op, std::size_t numShots, std::size_t numRounds, std::optional noise = std::nullopt) { + ensure_feedback_sampling_target_supported(code); auto [synd, dataRes] = noise ? sample_memory_circuit(code, op, numShots, numRounds, *noise) : sample_memory_circuit(code, op, numShots, numRounds); @@ -571,7 +684,10 @@ void bindCode(nb::module_ &mod) { "Sample the memory circuit of the code with a specific initial " "operation. Returns (syndromes, data): syndromes has numFixed " "boundary detectors (basis matching op), then one block per " - "inter-round transition, then numFixed more boundary detectors.", + "inter-round transition, then numFixed more boundary detectors. data " + "contains raw final data-qubit measurements; observable inlined " + "feedback is represented in generated detector error models and is not " + "folded into data as a corrected logical-observable value.", nb::arg("code"), nb::arg("op"), nb::arg("numShots"), nb::arg("numRounds"), nb::arg("noise") = nb::none()); @@ -580,6 +696,7 @@ void bindCode(nb::module_ &mod) { std::size_t numRounds, std::optional noise = std::nullopt) { + ensure_feedback_sampling_target_supported(code); cudaq::noise_model emptyNoise; auto [synd, dataRes] = samplingFn(code, op, numShots, numRounds, noise ? *noise : emptyNoise); diff --git a/libs/qec/python/tests/test_code.py b/libs/qec/python/tests/test_code.py index aa1603774..7c641a7db 100644 --- a/libs/qec/python/tests/test_code.py +++ b/libs/qec/python/tests/test_code.py @@ -337,5 +337,219 @@ def test_version(): assert "CUDA-Q QEC" in qec.__version__ +# Kernels for the inlined-feedback toy code: two data qubits with stabilizers +# XX and ZZ, both extracted through a single superdense (Bell-pair) ancilla +# pair (ancx[0] carries the XX record, ancz[0] the ZZ record). The round ends +# with an *uncorrected* record-conditioned byproduct: after the Bell decode, +# ancx[0] sits in the computational basis holding its future measurement +# outcome r_X, so the trailing CX(ancx[0], data[1]) is the unitary equivalent +# of the classically controlled byproduct X^{r_X} on data[1] that a hardware +# frame update would otherwise track. Same gadget and feedback matrices as the +# C++ toy in unittests/feedback_toy_device.cpp / test_qec.cpp (see the full +# derivation there): feedback = [[0, 1], [0, 0]], observable feedback = +# [[0, 1]], with records per round in [Z][X] order. +@cudaq.kernel +def _feedback_toy_prep0(q: patch): + reset(q.data[0]) + reset(q.data[1]) + + +@cudaq.kernel +def _feedback_toy_round(q: patch, x_stabilizers: list[int], + z_stabilizers: list[int]) -> list[cudaq.measure_handle]: + # Bell-prepare the superdense ancilla pair. + h(q.ancx[0]) + x.ctrl(q.ancx[0], q.ancz[0]) + # Couple XX through the X ancilla and ZZ through the Z ancilla. + x.ctrl(q.ancx[0], q.data[0]) + x.ctrl(q.ancx[0], q.data[1]) + x.ctrl(q.data[0], q.ancz[0]) + x.ctrl(q.data[1], q.ancz[0]) + # Decode the Bell pair. + x.ctrl(q.ancx[0], q.ancz[0]) + h(q.ancx[0]) + # Uncorrected record-conditioned byproduct: X^{r_X} on data[1]. + x.ctrl(q.ancx[0], q.data[1]) + # Records in [Z][X] order. + results = mz([*q.ancz, *q.ancx]) + reset(q.ancz[0]) + reset(q.ancx[0]) + return results + + +def _register_feedback_toy(name, + detector_feedback=None, + observable_feedback_z=None, + observable_feedback_x=None): + """Register one variant of the common Python feedback toy.""" + + class FeedbackToy: + + def __init__(self): + qec.Code.__init__(self) + self.stabilizers = [ + cudaq.SpinOperator.from_word(w) for w in ["XX", "ZZ"] + ] + self.pauli_observables = [cudaq.SpinOperator.from_word("ZZ")] + self.operation_encodings = { + qec.operation.prep0: _feedback_toy_prep0, + qec.operation.stabilizer_round: _feedback_toy_round + } + + def get_num_data_qubits(self): + return 2 + + def get_num_ancilla_qubits(self): + return 2 + + def get_num_ancilla_x_qubits(self): + return 1 + + def get_num_ancilla_z_qubits(self): + return 1 + + def get_num_x_stabilizers(self): + return 1 + + def get_num_z_stabilizers(self): + return 1 + + if detector_feedback is not None: + FeedbackToy.get_inlined_feedback = lambda self: detector_feedback() + if observable_feedback_z is not None: + FeedbackToy.get_observable_inlined_feedback_z = \ + lambda self: observable_feedback_z() + if observable_feedback_x is not None: + FeedbackToy.get_observable_inlined_feedback_x = \ + lambda self: observable_feedback_x() + qec.code(name)(FeedbackToy) + + +def _valid_detector_feedback(): + return np.array([[0, 1], [0, 0]], dtype=np.uint8) + + +def _valid_observable_feedback(): + # Exercise the bridge's dtype coercion: int64 instead of uint8. + return np.array([[0, 1]]) + + +def test_python_inlined_feedback_toy(): + _register_feedback_toy('py-feedback-toy', _valid_detector_feedback, + _valid_observable_feedback) + + # Sample on stim, mirroring the C++ unit tests (test_qec links the stim + # target). qpp-cpu cannot preserve the cross-round correlations required by + # this Python kernel and is rejected explicitly by the binding below. + cudaq.set_target('stim') + try: + code = qec.get_code('py-feedback-toy') + numShots = 20 + for numRounds in (1, 2, 3, 4): + syndromes, data = qec.sample_memory_circuit(code, + numShots=numShots, + numRounds=numRounds) + + # 1 first-round boundary + 2 * (numRounds - 1) cross-round + 1 + # final boundary detectors, all deterministic (zero) in the + # noiseless circuit. + assert syndromes.shape == (numShots, 2 * numRounds) + assert not np.any(syndromes) + + # data contains raw final measurements; observable feedback is + # folded into the DEM observable, not into this returned tensor. + assert data.shape == (numShots, 2) + finally: + cudaq.reset_target() + + +def test_python_inlined_feedback_rejects_unsupported_qpp_sampling(): + _register_feedback_toy('py-feedback-toy-qpp', _valid_detector_feedback, + _valid_observable_feedback) + cudaq.set_target('qpp-cpu') + try: + with pytest.raises(RuntimeError, + match="Use cudaq.set_target\\('stim'\\)"): + qec.sample_memory_circuit(qec.get_code('py-feedback-toy-qpp'), + numShots=1, + numRounds=1) + finally: + cudaq.reset_target() + + +def test_python_inlined_feedback_toy_negative_control(): + # The identical toy without the feedback declarations: the uncorrected + # byproduct makes the record-0 detectors non-deterministic and stim must + # reject the circuit. This proves the Python-declared feedback is doing + # the work in the positive test above. + _register_feedback_toy('py-feedback-toy-nofb') + + # stim's determinism analysis (a std::invalid_argument surfaced as + # ValueError) rejects the circuit inside dem_from_kernel, before any + # sampling happens. + with pytest.raises(ValueError, match="non-deterministic detectors"): + qec.sample_memory_circuit(qec.get_code('py-feedback-toy-nofb'), + numShots=20, + numRounds=4) + + +@pytest.mark.parametrize("invalid_target", ["detector", "observable"]) +def test_python_inlined_feedback_rejects_non_binary_values(invalid_target): + + def invalid_detector_feedback(): + # All of these values previously truncated or wrapped during the eager + # uint8 conversion in the Python bridge. + return np.array([[2, -1], [0.5, 256]]) + + def invalid_observable_feedback(): + return np.array([[0, 2]]) + + name = f'py-feedback-toy-non-binary-{invalid_target}' + _register_feedback_toy( + name, invalid_detector_feedback if invalid_target == "detector" else + _valid_detector_feedback, invalid_observable_feedback + if invalid_target == "observable" else _valid_observable_feedback) + + cudaq.set_target('stim') + try: + with pytest.raises(RuntimeError, + match="must contain only exact 0 or 1 entries"): + qec.sample_memory_circuit(qec.get_code(name), + numShots=1, + numRounds=1) + finally: + cudaq.reset_target() + + +def test_python_inlined_feedback_observable_basis_selection(): + # The bridge routes get_observable_inlined_feedback_x to the X-basis getter, + # so it is never queried on the Z (prep0) path. Plant a deliberately + # wrong-shape _x matrix alongside a valid _z matrix: were the bridge to feed + # _x into the Z path, flatten_feedback_tensor's shape check would reject the + # circuit. The toy registers only prep0, so this proves selection from the + # Z path (the wrong-shape _x is inert there). + def wrong_shape_observable_feedback(): + # Wrong shape for [num_obs=1 x numCols=2]; must never reach the Z-basis + # flatten on the prep0 path. + return np.zeros((3, 5), dtype=np.uint8) + + _register_feedback_toy('py-feedback-toy-basis-select', + _valid_detector_feedback, _valid_observable_feedback, + wrong_shape_observable_feedback) + + cudaq.set_target('stim') + try: + numShots, numRounds = 20, 4 + syndromes, data = qec.sample_memory_circuit( + qec.get_code('py-feedback-toy-basis-select'), + numShots=numShots, + numRounds=numRounds) + assert syndromes.shape == (numShots, 2 * numRounds) + assert not np.any(syndromes) + assert data.shape == (numShots, 2) + finally: + cudaq.reset_target() + + if __name__ == "__main__": pytest.main() diff --git a/libs/qec/unittests/CMakeLists.txt b/libs/qec/unittests/CMakeLists.txt index 359258a51..dcae72f26 100644 --- a/libs/qec/unittests/CMakeLists.txt +++ b/libs/qec/unittests/CMakeLists.txt @@ -62,7 +62,16 @@ add_dependencies(CUDAQXQECUnitTests test_decoding_server_core) gtest_discover_tests(test_decoding_server_core) add_executable(test_qec test_qec.cpp) -target_link_libraries(test_qec PRIVATE GTest::gtest_main cudaq-qec cudaq::cudaq-stim-target) +# The inlined-feedback toy code's kernels are nvq++-compiled (like the +# registered codes' *_device.cpp files) so the memory-circuit kernel can +# resolve them through the kernel registry; cudaq::cudaq and nvqir provide +# the runtime symbols that the device object references. +cudaqx_add_device_code(test_qec + SOURCES + feedback_toy_device.cpp +) +target_link_directories(test_qec PRIVATE ${CUDAQ_INSTALL_DIR}/lib) +target_link_libraries(test_qec PRIVATE GTest::gtest_main cudaq-qec cudaq::cudaq nvqir cudaq::cudaq-stim-target) add_dependencies(CUDAQXQECUnitTests test_qec) gtest_discover_tests(test_qec) diff --git a/libs/qec/unittests/feedback_toy_device.cpp b/libs/qec/unittests/feedback_toy_device.cpp new file mode 100644 index 000000000..0eb4a657e --- /dev/null +++ b/libs/qec/unittests/feedback_toy_device.cpp @@ -0,0 +1,64 @@ +/****************************************************************-*- C++ -*-**** + * Copyright (c) 2026 NVIDIA Corporation & Affiliates. * + * All rights reserved. * + * * + * This source code and the accompanying materials are made available under * + * the terms of the Apache License 2.0 which accompanies this distribution. * + ******************************************************************************/ + +#include "cudaq.h" +#include "cudaq/qec/patch.h" + +// Device kernels for the inlined-feedback toy code used by test_qec.cpp +// (see the QECCodeTester.checkInlinedFeedbackToy* tests). This file is +// nvq++-compiled, like the registered codes' *_device.cpp files, so the +// kernels are visible to the kernel registry that the memory-circuit +// machinery uses to resolve qkernel arguments. + +namespace cudaq::qec::feedback_toy { + +__qpu__ void prep0(patch p) { + for (std::size_t i = 0; i < p.data.size(); i++) + reset(p.data[i]); +} + +__qpu__ void prepp(patch p) { + for (std::size_t i = 0; i < p.data.size(); i++) { + reset(p.data[i]); + h(p.data[i]); + } +} + +// A superdense (Bell-pair) round for a 2-data-qubit code with stabilizers XX +// and ZZ that deliberately ends with an *uncorrected* record-conditioned +// byproducts: after the Bell decode, the ancillas sit in the computational +// basis holding their future measurement outcomes r_X and r_Z. The trailing +// controlled gates are the unitary equivalent of classically controlled +// X^{r_X} Z^{r_Z} byproducts on data[1] that a hardware frame update would +// otherwise track. The stabilizer supports are hardcoded (a single XX and ZZ +// plaquette), so the x/z_stabilizers arguments are intentionally unused. +__qpu__ std::vector +stabilizer_round(patch p, const std::vector &x_stabilizers, + const std::vector &z_stabilizers) { + // Bell-prepare the superdense ancilla pair. + h(p.ancx[0]); + cudaq::x(p.ancx[0], p.ancz[0]); + // Couple XX through the X ancilla and ZZ through the Z ancilla. + cudaq::x(p.ancx[0], p.data[0]); + cudaq::x(p.ancx[0], p.data[1]); + cudaq::x(p.data[0], p.ancz[0]); + cudaq::x(p.data[1], p.ancz[0]); + // Decode the Bell pair. + cudaq::x(p.ancx[0], p.ancz[0]); + h(p.ancx[0]); + // Uncorrected record-conditioned byproducts on data[1]. + cudaq::x(p.ancx[0], p.data[1]); + cudaq::z(p.ancz[0], p.data[1]); + // Records in [Z][X] order. + auto results = mz(p.ancz, p.ancx); + reset(p.ancz[0]); + reset(p.ancx[0]); + return results; +} + +} // namespace cudaq::qec::feedback_toy diff --git a/libs/qec/unittests/test_qec.cpp b/libs/qec/unittests/test_qec.cpp index 01d3aff31..ea99cedfc 100644 --- a/libs/qec/unittests/test_qec.cpp +++ b/libs/qec/unittests/test_qec.cpp @@ -11,9 +11,12 @@ #include #include #include +#include #include #include "cuda-qx/core/library_utils.h" +#include "../lib/inlined_feedback_layout.h" +#include "cudaq/algorithms/dem.h" #include "cudaq/qec/codes/surface_code.h" #include "cudaq/qec/experiments.h" #include "cudaq/qec/pcm_utils.h" @@ -2297,3 +2300,560 @@ TEST(PluginLoaderTester, checkCleanupPluginsEdgeCases) { // not do anything. cudaq::qec::cleanup_plugins(cudaq::qec::PluginType::CODE); } + +TEST(QECCodeTester, checkInlinedFeedbackDefaults) { + auto steane = cudaq::qec::get_code("steane"); + auto feedback = steane->get_inlined_feedback(); + EXPECT_EQ(feedback.rank(), 0); + EXPECT_EQ(feedback.size(), 0); + auto obs_feedback_z = steane->get_observable_inlined_feedback_z(); + EXPECT_EQ(obs_feedback_z.rank(), 0); + EXPECT_EQ(obs_feedback_z.size(), 0); + auto obs_feedback_x = steane->get_observable_inlined_feedback_x(); + EXPECT_EQ(obs_feedback_x.rank(), 0); + EXPECT_EQ(obs_feedback_x.size(), 0); +} + +// Device kernels for the inlined-feedback toy code, nvq++-compiled in +// feedback_toy_device.cpp (host-compiled kernels cannot be resolved through +// the kernel registry that the memory-circuit kernel uses for its qkernel +// arguments). +namespace cudaq::qec::feedback_toy { +__qpu__ void prep0(patch p); +__qpu__ void prepp(patch p); +__qpu__ std::vector +stabilizer_round(patch p, const std::vector &x_stabilizers, + const std::vector &z_stabilizers); +} // namespace cudaq::qec::feedback_toy + +namespace { + +// --------------------------------------------------------------------------- +// Toy code for the inlined-feedback tests below. +// +// Two data qubits with stabilizers XX and ZZ and logical observables XX/ZZ. +// Both stabilizers are extracted through a single superdense (Bell-pair) +// ancilla pair: ancx[0] carries the XX record and ancz[0] the ZZ record. The +// round (defined in feedback_toy_device.cpp) deliberately ends with uncorrected +// record-conditioned X^{r_X} Z^{r_Z} byproducts on data[1]. +// +// Derivation of the feedback matrices (records per round in [Z][X] order: +// record 0 = ancz outcome r_Z, record 1 = ancx outcome r_X): +// - In a Z-basis memory, X^{r_X} changes the later ZZ parity. Record 0's +// detector therefore includes the earlier record 1, and the ZZ observable +// includes record 1 from every round. +// - In an X-basis memory, Z^{r_Z} changes the later XX parity. Record 1's +// detector therefore includes the earlier record 0, and the XX observable +// includes record 0 from every round. +// Thus detector feedback is [[0,1],[1,0]], Z-observable feedback is [[0,1]], +// and X-observable feedback is [[1,0]]. +// --------------------------------------------------------------------------- + +class feedback_toy_code : public cudaq::qec::code { +protected: + std::size_t get_num_data_qubits() const override { return 2; } + std::size_t get_num_ancilla_qubits() const override { return 2; } + std::size_t get_num_ancilla_x_qubits() const override { return 1; } + std::size_t get_num_ancilla_z_qubits() const override { return 1; } + std::size_t get_num_x_stabilizers() const override { return 1; } + std::size_t get_num_z_stabilizers() const override { return 1; } + + bool declare_feedback; + uint8_t detector_feedback_value; + uint8_t observable_feedback_value; + +public: + explicit feedback_toy_code(bool declare_feedback, + uint8_t detector_feedback_value = 1, + uint8_t observable_feedback_value = 1) + : declare_feedback(declare_feedback), + detector_feedback_value(detector_feedback_value), + observable_feedback_value(observable_feedback_value) { + operation_encodings.insert(std::make_pair(cudaq::qec::operation::prep0, + cudaq::qec::feedback_toy::prep0)); + operation_encodings.insert(std::make_pair(cudaq::qec::operation::prepp, + cudaq::qec::feedback_toy::prepp)); + operation_encodings.insert( + std::make_pair(cudaq::qec::operation::stabilizer_round, + cudaq::qec::feedback_toy::stabilizer_round)); + m_stabilizers = fromPauliWords({"XX", "ZZ"}); + m_pauli_observables = fromPauliWords({"XX", "ZZ"}); + } + + cudaqx::tensor get_inlined_feedback() const override { + if (!declare_feedback) + return code::get_inlined_feedback(); + cudaqx::tensor feedback({2, 2}); + feedback.at({0, 1}) = detector_feedback_value; + feedback.at({1, 0}) = detector_feedback_value; + return feedback; + } + + cudaqx::tensor get_observable_inlined_feedback_z() const override { + if (!declare_feedback) + return code::get_observable_inlined_feedback_z(); + cudaqx::tensor feedback({1, 2}); + feedback.at({0, 1}) = observable_feedback_value; + return feedback; + } + + cudaqx::tensor get_observable_inlined_feedback_x() const override { + if (!declare_feedback) + return code::get_observable_inlined_feedback_x(); + cudaqx::tensor feedback({1, 2}); + feedback.at({0, 0}) = observable_feedback_value; + return feedback; + } +}; + +// Basis-selection variant of the toy: returns caller-planted observable- +// feedback matrices for the two bases so a test can put a deliberately +// wrong-shape matrix on the basis it expects sample_memory_circuit NOT to +// query. The detector feedback (get_inlined_feedback) stays valid, so the +// underlying Z-memory circuit is otherwise identical to the positive toy. +class feedback_toy_basis_select_code : public feedback_toy_code { + cudaqx::tensor m_obs_z; + cudaqx::tensor m_obs_x; + +public: + feedback_toy_basis_select_code(cudaqx::tensor obs_z, + cudaqx::tensor obs_x) + : feedback_toy_code(/*declare_feedback=*/true), m_obs_z(std::move(obs_z)), + m_obs_x(std::move(obs_x)) {} + + cudaqx::tensor get_observable_inlined_feedback_z() const override { + return m_obs_z; + } + cudaqx::tensor get_observable_inlined_feedback_x() const override { + return m_obs_x; + } +}; + +} // namespace + +TEST(QECCodeTester, checkInlinedFeedbackToyMemoryCircuit) { + feedback_toy_code toy(/*declare_feedback=*/true); + const std::size_t numShots = 20; + for (const std::size_t numRounds : {1UL, 2UL, 3UL, 4UL}) { + auto [syndromes, data] = cudaq::qec::sample_memory_circuit( + toy, cudaq::qec::operation::prep0, numShots, numRounds); + + // 1 first-round boundary + 2 * (numRounds - 1) cross-round + 1 final + // boundary detectors. + ASSERT_EQ(syndromes.rank(), 2); + ASSERT_EQ(syndromes.shape()[0], numShots); + ASSERT_EQ(syndromes.shape()[1], 2 * numRounds); + for (std::size_t shot = 0; shot < numShots; ++shot) + for (std::size_t d = 0; d < syndromes.shape()[1]; ++d) + EXPECT_EQ(syndromes.at({shot, d}), 0) + << "rounds " << numRounds << ", shot " << shot << ", detector " + << d; + + // data contains raw final measurements. In particular, it does not fold + // get_observable_inlined_feedback_z() into a corrected logical result. + ASSERT_EQ(data.rank(), 2); + ASSERT_EQ(data.shape()[0], numShots); + ASSERT_EQ(data.shape()[1], 2); + } +} + +TEST(QECCodeTester, checkInlinedFeedbackToyXMemoryCircuit) { + feedback_toy_code toy(/*declare_feedback=*/true); + const std::size_t numShots = 20; + for (const std::size_t numRounds : {1UL, 2UL, 3UL, 4UL}) { + auto [syndromes, data] = cudaq::qec::sample_memory_circuit( + toy, cudaq::qec::operation::prepp, numShots, numRounds); + + ASSERT_EQ(syndromes.rank(), 2); + ASSERT_EQ(syndromes.shape()[0], numShots); + ASSERT_EQ(syndromes.shape()[1], 2 * numRounds); + for (std::size_t shot = 0; shot < numShots; ++shot) + for (std::size_t d = 0; d < syndromes.shape()[1]; ++d) + EXPECT_EQ(syndromes.at({shot, d}), 0) + << "rounds " << numRounds << ", shot " << shot << ", detector " + << d; + + ASSERT_EQ(data.rank(), 2); + ASSERT_EQ(data.shape()[0], numShots); + ASSERT_EQ(data.shape()[1], 2); + } +} + +TEST(QECCodeTester, checkInlinedFeedbackObservableBasisSelection) { + auto valid_z_obs = [] { + cudaqx::tensor t({1, 2}); + t.at({0, 1}) = 1; + return t; + }; + auto valid_x_obs = [] { + cudaqx::tensor t({1, 2}); + t.at({0, 0}) = 1; + return t; + }; + // Non-empty and the wrong shape for [num_obs=1 x numCols=2], so any code that + // selects it trips the host layout's shape check. + auto wrong_shape_obs = [] { return cudaqx::tensor({3, 5}); }; + + // A wrong-shape matrix for the unselected basis must not affect the run. + { + feedback_toy_basis_select_code toy(valid_z_obs(), wrong_shape_obs()); + const std::size_t numShots = 20; + const std::size_t numRounds = 4; + auto [syndromes, data] = cudaq::qec::sample_memory_circuit( + toy, cudaq::qec::operation::prep0, numShots, numRounds); + ASSERT_EQ(syndromes.rank(), 2); + ASSERT_EQ(syndromes.shape()[0], numShots); + ASSERT_EQ(syndromes.shape()[1], 2 * numRounds); + for (std::size_t shot = 0; shot < numShots; ++shot) + for (std::size_t d = 0; d < syndromes.shape()[1]; ++d) + EXPECT_EQ(syndromes.at({shot, d}), 0) + << "shot " << shot << ", detector " << d; + ASSERT_EQ(data.shape()[0], numShots); + ASSERT_EQ(data.shape()[1], 2); + } + { + feedback_toy_basis_select_code toy(wrong_shape_obs(), valid_x_obs()); + const std::size_t numShots = 20; + const std::size_t numRounds = 4; + auto [syndromes, data] = cudaq::qec::sample_memory_circuit( + toy, cudaq::qec::operation::prepp, numShots, numRounds); + ASSERT_EQ(syndromes.rank(), 2); + ASSERT_EQ(syndromes.shape()[0], numShots); + ASSERT_EQ(syndromes.shape()[1], 2 * numRounds); + for (std::size_t shot = 0; shot < numShots; ++shot) + for (std::size_t d = 0; d < syndromes.shape()[1]; ++d) + EXPECT_EQ(syndromes.at({shot, d}), 0) + << "shot " << shot << ", detector " << d; + ASSERT_EQ(data.shape()[0], numShots); + ASSERT_EQ(data.shape()[1], 2); + } + + // A wrong-shape matrix for the selected basis must be rejected with the + // corresponding getter name. + { + feedback_toy_basis_select_code toy(wrong_shape_obs(), valid_x_obs()); + try { + cudaq::qec::sample_memory_circuit(toy, cudaq::qec::operation::prep0, + /*numShots=*/20, /*numRounds=*/4); + FAIL() << "expected sample_memory_circuit to reject the wrong-shape " + "Z-basis observable feedback"; + } catch (const std::exception &e) { + EXPECT_NE(std::string(e.what()).find( + "get_observable_inlined_feedback_z() has invalid shape"), + std::string::npos) + << "unexpected failure reason: " << e.what(); + } + } + { + feedback_toy_basis_select_code toy(valid_z_obs(), wrong_shape_obs()); + try { + cudaq::qec::sample_memory_circuit(toy, cudaq::qec::operation::prepp, + /*numShots=*/20, /*numRounds=*/4); + FAIL() << "expected sample_memory_circuit to reject the wrong-shape " + "X-basis observable feedback"; + } catch (const std::exception &e) { + EXPECT_NE(std::string(e.what()).find( + "get_observable_inlined_feedback_x() has invalid shape"), + std::string::npos) + << "unexpected failure reason: " << e.what(); + } + } +} + +TEST(QECCodeTester, checkInlinedFeedbackRejectsNonBinaryDeclarations) { + auto expect_non_binary_rejection = [](const feedback_toy_code &toy, + const std::string &getter_name) { + try { + cudaq::qec::sample_memory_circuit(toy, cudaq::qec::operation::prep0, + /*numShots=*/1, /*numRounds=*/1); + FAIL() << "expected sample_memory_circuit to reject " << getter_name; + } catch (const std::exception &e) { + EXPECT_NE( + std::string(e.what()).find(getter_name + " has non-binary value 2"), + std::string::npos) + << "unexpected failure reason: " << e.what(); + } + }; + + expect_non_binary_rejection(feedback_toy_code(/*declare_feedback=*/true, + /*detector_feedback_value=*/2), + "get_inlined_feedback()"); + expect_non_binary_rejection( + feedback_toy_code(/*declare_feedback=*/true, + /*detector_feedback_value=*/1, + /*observable_feedback_value=*/2), + "get_observable_inlined_feedback_z()"); +} + +TEST(QECCodeTester, checkInlinedFeedbackToyNegativeControl) { + // The identical toy without the feedback declaration: the uncorrected + // byproducts make the matching-basis detectors non-deterministic and stim + // must reject both circuits. + feedback_toy_code toy(/*declare_feedback=*/false); + for (auto state_prep : + {cudaq::qec::operation::prep0, cudaq::qec::operation::prepp}) { + try { + cudaq::qec::sample_memory_circuit(toy, state_prep, + /*numShots=*/20, /*numRounds=*/4); + FAIL() + << "expected sample_memory_circuit to reject the feedback-less toy"; + } catch (const std::exception &e) { + EXPECT_NE(std::string(e.what()).find("non-deterministic detectors"), + std::string::npos) + << "unexpected failure reason: " << e.what(); + } + } +} + +TEST(QECCodeTester, checkInlinedFeedbackToyDem) { + // The DEM path threads the same feedback declaration through + // dem_from_kernel; with circuit-level noise it must produce a valid model + // (a noiseless model has no error mechanisms and is rejected downstream). + feedback_toy_code toy(/*declare_feedback=*/true); + const std::size_t numRounds = 3; + cudaq::noise_model noise; + noise.add_all_qubit_channel("x", cudaq::qec::two_qubit_depolarization(0.01), + /*num_controls=*/1); + auto dem = cudaq::qec::dem_from_memory_circuit( + toy, cudaq::qec::operation::prep0, numRounds, noise); + // 1 first-round boundary + 2 * (numRounds - 1) cross-round + 1 final + // boundary detectors. + EXPECT_EQ(dem.num_detectors(), 2 * numRounds); + EXPECT_GT(dem.num_error_mechanisms(), 0); + EXPECT_EQ(dem.num_observables(), 1); +} + +TEST(InlinedFeedbackLayout, EmptyTensorsGiveEmptyLayout) { + cudaqx::tensor empty; + auto layout = + cudaq::qec::details::build_inlined_feedback_layout(empty, empty, 2, 1); + EXPECT_TRUE(layout.detector_indices.empty()); + EXPECT_TRUE(layout.detector_offsets.empty()); + EXPECT_TRUE(layout.observable_indices.empty()); + EXPECT_TRUE(layout.observable_offsets.empty()); +} + +TEST(InlinedFeedbackLayout, ToyMatricesCsr) { + // fb = [[0,1],[0,0]]: record 0 heralded by record 1; row 1 empty. + cudaqx::tensor fb({2, 2}); + fb.at({0, 1}) = 1; + // obs_fb = [[0,1]]: observable 0 XORs record 1 every round. + cudaqx::tensor obs_fb({1, 2}); + obs_fb.at({0, 1}) = 1; + auto layout = + cudaq::qec::details::build_inlined_feedback_layout(fb, obs_fb, 2, 1); + EXPECT_EQ(layout.detector_indices, (std::vector{1})); + EXPECT_EQ(layout.detector_offsets, (std::vector{0, 1, 1})); + EXPECT_EQ(layout.observable_indices, (std::vector{1})); + EXPECT_EQ(layout.observable_offsets, (std::vector{0, 1})); +} + +TEST(InlinedFeedbackLayout, MultiRowWeights) { + // 3x3 with row weights 2, 0, 1. + cudaqx::tensor fb({3, 3}); + fb.at({0, 0}) = 1; + fb.at({0, 2}) = 1; + fb.at({2, 1}) = 1; + cudaqx::tensor empty; + auto layout = + cudaq::qec::details::build_inlined_feedback_layout(fb, empty, 3, 0); + EXPECT_EQ(layout.detector_indices, (std::vector{0, 2, 1})); + EXPECT_EQ(layout.detector_offsets, (std::vector{0, 2, 2, 3})); + EXPECT_TRUE(layout.observable_offsets.empty()); +} + +TEST(InlinedFeedbackLayout, AllZeroMatrixGivesZeroWeightRows) { + cudaqx::tensor fb({2, 2}); + cudaqx::tensor empty; + auto layout = + cudaq::qec::details::build_inlined_feedback_layout(fb, empty, 2, 0); + EXPECT_TRUE(layout.detector_indices.empty()); + EXPECT_EQ(layout.detector_offsets, (std::vector{0, 0, 0})); +} + +TEST(InlinedFeedbackLayout, ThrowsOnWrongShape) { + cudaqx::tensor bad({2, 3}); // expected [2 x 2] + cudaqx::tensor empty; + EXPECT_THROW( + cudaq::qec::details::build_inlined_feedback_layout(bad, empty, 2, 0), + std::runtime_error); + cudaqx::tensor bad_obs({2, 2}); // expected [1 x 2] + EXPECT_THROW( + cudaq::qec::details::build_inlined_feedback_layout(empty, bad_obs, 2, 1), + std::runtime_error); +} + +TEST(InlinedFeedbackLayout, ThrowsOnNonBinaryEntries) { + cudaqx::tensor empty; + + cudaqx::tensor bad_detector({2, 2}); + bad_detector.at({1, 0}) = 2; + try { + cudaq::qec::details::build_inlined_feedback_layout(bad_detector, empty, 2, + 0); + FAIL() << "expected non-binary detector feedback to be rejected"; + } catch (const std::exception &e) { + EXPECT_NE(std::string(e.what()).find("non-binary value 2 at flat index 2"), + std::string::npos) + << "unexpected failure reason: " << e.what(); + } + + cudaqx::tensor bad_observable({1, 2}); + bad_observable.at({0, 1}) = 255; + try { + cudaq::qec::details::build_inlined_feedback_layout(empty, bad_observable, 2, + 1); + FAIL() << "expected non-binary observable feedback to be rejected"; + } catch (const std::exception &e) { + EXPECT_NE( + std::string(e.what()).find("non-binary value 255 at flat index 1"), + std::string::npos) + << "unexpected failure reason: " << e.what(); + } +} + +// The entry-point memory-circuit kernel (defined in device/memory_circuit.cpp); +// forward-declared here - like the feedback_toy kernels above - so the +// conformance test below can drive dem_from_kernel directly. +namespace cudaq::qec { +__qpu__ void +memory_circuit(const code::stabilizer_round &stabilizer_round, + const code::one_qubit_encoding &statePrep, std::size_t numData, + std::size_t numAncx, std::size_t numAncz, std::size_t numRounds, + const std::vector &x_stabilizers, + const std::vector &z_stabilizers, + const std::vector &obs_matrix_flat, + std::size_t num_observables, bool measure_in_x_basis, + const std::vector &feedback_indices, + const std::vector &feedback_offsets, + const std::vector &obs_feedback_indices, + const std::vector &obs_feedback_offsets); +} // namespace cudaq::qec + +// Conformance: the device fold inside memory_circuit and the host-side rule in +// build_inlined_feedback_layout must agree on which records compose each +// detector and observable. Build M2D and M2O from the kernel and compare them, +// row by row, against the matrices reconstructed from the host layout plus the +// known emission order. If either fold drifts, this test fails. +static void +checkHostLayoutMatchesKernelM2DAndM2O(cudaq::qec::operation statePrep, + bool measureInXBasis) { + SCOPED_TRACE(measureInXBasis ? "X memory" : "Z memory"); + feedback_toy_code toy(/*declare_feedback=*/true); + const std::size_t numRounds = 3; + const std::size_t numData = 2, numAncx = 1, numAncz = 1; + const std::size_t numCols = numAncx + numAncz; // records per round + const std::size_t num_obs = 1; + + // --- Kernel-derived M2D, mirroring experiments.cpp's sample path. --- + auto &prep = + toy.get_operation(statePrep); + auto &stabRound = toy.get_operation( + cudaq::qec::operation::stabilizer_round); + auto parity_x = toy.get_parity_x(); + auto parity_z = toy.get_parity_z(); + std::vector xVec(parity_x.data(), + parity_x.data() + parity_x.size()); + std::vector zVec(parity_z.data(), + parity_z.data() + parity_z.size()); + auto logical_obs = + measureInXBasis ? toy.get_observables_x() : toy.get_observables_z(); + std::vector obs_flat(logical_obs.data(), + logical_obs.data() + logical_obs.size()); + + auto feedback = toy.get_inlined_feedback(); + auto obs_feedback = measureInXBasis ? toy.get_observable_inlined_feedback_x() + : toy.get_observable_inlined_feedback_z(); + auto layout = cudaq::qec::details::build_inlined_feedback_layout( + feedback, obs_feedback, numCols, num_obs); + + cudaq::noise_model noise; + cudaq::M2DSparseMatrix m2d; + cudaq::M2OSparseMatrix m2o; + cudaq::dem_from_kernel( + cudaq::qec::memory_circuit, &noise, /*options=*/{}, m2d, m2o, stabRound, + prep, numData, numAncx, numAncz, numRounds, xVec, zVec, obs_flat, num_obs, + measureInXBasis, layout.detector_indices, layout.detector_offsets, + layout.observable_indices, layout.observable_offsets); + + // --- Host-layout-derived expected M2D. Global record index for round r, + // column c is r*numCols + c; the readout for data qubit q is at + // numRounds*numCols + q. --- + auto rec = [numCols](std::size_t round, std::size_t col) { + return round * numCols + col; + }; + auto herald_cols = [](const std::vector &indices, + const std::vector &offsets, + std::size_t row) { + std::vector cols; + if (!offsets.empty()) + cols.assign(indices.begin() + offsets[row], + indices.begin() + offsets[row + 1]); + return cols; + }; + + std::vector> expected; + const std::size_t fixedOffset = measureInXBasis ? numAncz : 0; + const auto &fixedStabilizers = measureInXBasis ? xVec : zVec; + // Round-0 fixed detector for the stabilizer type matching the memory basis. + expected.push_back({rec(0, fixedOffset)}); + // Interior cross-round detectors, emitted each round for j = 0..numCols-1: + // {prev_j, curr_j} plus the layout-row-j heralds taken from the earlier + // round. + for (std::size_t r = 1; r < numRounds; ++r) { + for (std::size_t j = 0; j < numCols; ++j) { + std::set row{rec(r - 1, j), rec(r, j)}; + for (auto k : + herald_cols(layout.detector_indices, layout.detector_offsets, j)) + row.insert(rec(r - 1, k)); + expected.push_back(row); + } + } + // Final boundary detector for the fixed record: the last-round fixed record, + // data readouts in the stabilizer support, and the layout-row heralds from + // the last round. + { + std::set row{rec(numRounds - 1, fixedOffset)}; + for (std::size_t q = 0; q < numData; ++q) + if (fixedStabilizers[q] != 0) + row.insert(numRounds * numCols + q); + for (auto k : herald_cols(layout.detector_indices, layout.detector_offsets, + fixedOffset)) + row.insert(rec(numRounds - 1, k)); + expected.push_back(row); + } + + ASSERT_EQ(m2d.rows.size(), 2 * numRounds); + ASSERT_EQ(m2d.rows.size(), expected.size()); + for (std::size_t d = 0; d < expected.size(); ++d) { + std::set actual(m2d.rows[d].begin(), m2d.rows[d].end()); + EXPECT_EQ(actual, expected[d]) << "detector " << d; + } + + // Observable m contains every per-round record selected by observable row m + // followed by the final data measurements in logical-observable row m. + ASSERT_EQ(m2o.rows.size(), num_obs); + for (std::size_t obs = 0; obs < num_obs; ++obs) { + std::set obs_expected; + for (std::size_t r = 0; r < numRounds; ++r) + for (auto k : herald_cols(layout.observable_indices, + layout.observable_offsets, obs)) + obs_expected.insert(rec(r, k)); + for (std::size_t q = 0; q < numData; ++q) + if (obs_flat[obs * numData + q] != 0) + obs_expected.insert(numRounds * numCols + q); + + std::set actual(m2o.rows[obs].begin(), m2o.rows[obs].end()); + EXPECT_EQ(actual, obs_expected) << "observable " << obs; + } +} + +TEST(InlinedFeedbackLayout, HostLayoutMatchesKernelM2DAndM2O) { + checkHostLayoutMatchesKernelM2DAndM2O(cudaq::qec::operation::prep0, + /*measureInXBasis=*/false); +} + +TEST(InlinedFeedbackLayout, XHostLayoutMatchesKernelM2DAndM2O) { + checkHostLayoutMatchesKernelM2DAndM2O(cudaq::qec::operation::prepp, + /*measureInXBasis=*/true); +}