From e80a08545b55927137a4997b1c66d1d127e6648f Mon Sep 17 00:00:00 2001 From: kvmto Date: Thu, 9 Jul 2026 17:43:28 +0200 Subject: [PATCH 1/4] Support code-declared inlined feedback in memory circuits Codes may declare record->detector feedback corrections as data via two new virtuals (get_inlined_feedback / get_observable_inlined_feedback, empty default = existing behavior); memory_circuit folds the declared records into its detector and observable definitions, so measurement schemes with readout byproducts get deterministic detectors through sample_memory_circuit without a physical correction gate. Co-Authored-By: Claude Fable 5 Signed-off-by: kvmto --- libs/qec/include/cudaq/qec/code.h | 24 +++ libs/qec/lib/code.cpp | 8 + libs/qec/lib/device/memory_circuit.cpp | 117 +++++++++++++-- libs/qec/lib/device/memory_circuit.h | 14 +- libs/qec/lib/experiments.cpp | 49 +++++- libs/qec/python/bindings/py_code.cpp | 61 +++++++- libs/qec/python/tests/test_code.py | 152 +++++++++++++++++++ libs/qec/unittests/CMakeLists.txt | 11 +- libs/qec/unittests/feedback_toy_device.cpp | 56 +++++++ libs/qec/unittests/test_qec.cpp | 164 +++++++++++++++++++++ 10 files changed, 639 insertions(+), 17 deletions(-) create mode 100644 libs/qec/unittests/feedback_toy_device.cpp diff --git a/libs/qec/include/cudaq/qec/code.h b/libs/qec/include/cudaq/qec/code.h index f7b9ea1dc..3a9eefce7 100644 --- a/libs/qec/include/cudaq/qec/code.h +++ b/libs/qec/include/cudaq/qec/code.h @@ -201,6 +201,30 @@ 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). 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 + /// + /// Shape is [num_observables x numCols], where numCols = + /// get_num_ancilla_qubits() (one measurement record per ancilla per + /// round, in [Z][X] order). Entry (m, k) = 1 means logical observable m + /// additionally XORs record k of every round. + /// @return Tensor representing the observable inlined feedback matrix. + /// An empty tensor (the default) means no feedback is applied to the + /// logical observables. + virtual cudaqx::tensor get_observable_inlined_feedback() const; + /// @brief Get the stabilizer generators /// @return Reference to stabilizers const std::vector &get_stabilizers() const { diff --git a/libs/qec/lib/code.cpp b/libs/qec/lib/code.cpp index 19693b8c5..bb6b720a6 100644 --- a/libs/qec/lib/code.cpp +++ b/libs/qec/lib/code.cpp @@ -61,6 +61,14 @@ 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() 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/memory_circuit.cpp b/libs/qec/lib/device/memory_circuit.cpp index b00754a28..b2d8e9ae3 100644 --- a/libs/qec/lib/device/memory_circuit.cpp +++ b/libs/qec/lib/device/memory_circuit.cpp @@ -18,7 +18,9 @@ __qpu__ void memory_circuit(const code::stabilizer_round &stabilizer_round, const std::vector &z_stabilizers, const std::vector &obs_matrix_flat, std::size_t num_observables, - bool measure_in_x_basis) { + bool measure_in_x_basis, + const std::vector &feedback_flat, + const std::vector &obs_feedback_flat) { // Allocate the data and ancilla qubits cudaq::qvector data(num_data), xstab_anc(numAncx), zstab_anc(numAncz); @@ -39,10 +41,77 @@ __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 prefix + // offsets: observable m's feedback records occupy + // [obs_fb_offset[m], obs_fb_offset[m] + num_rounds * obs_fb_weight[m]), + // laid out round-major within each observable's slice. + std::vector obs_fb_weight(num_observables); + std::vector obs_fb_offset(num_observables); + std::size_t obs_fb_total = 0; + if (obs_feedback_flat.size() > 0) { + for (std::size_t m = 0; m < num_observables; ++m) { + std::size_t weight = 0; + for (std::size_t k = 0; k < numCols; ++k) { + if (obs_feedback_flat[m * numCols + k] != 0) + weight++; + } + obs_fb_weight[m] = weight; + obs_fb_offset[m] = obs_fb_total; + obs_fb_total += num_rounds * weight; + } + } + std::vector obs_fb_records(obs_fb_total); + + // Collect the first-round feedback records for each observable. + if (obs_feedback_flat.size() > 0) { + for (std::size_t m = 0; m < num_observables; ++m) { + std::size_t idx = obs_fb_offset[m]; + for (std::size_t k = 0; k < numCols; ++k) { + if (obs_feedback_flat[m * numCols + k] != 0) + obs_fb_records[idx++] = final_syndrome[k]; + } + } + } + // 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_flat.size() > 0) { + for (std::size_t m = 0; m < num_observables; ++m) { + std::size_t idx = obs_fb_offset[m] + round * obs_fb_weight[m]; + for (std::size_t k = 0; k < numCols; ++k) { + if (obs_feedback_flat[m * numCols + k] != 0) + obs_fb_records[idx++] = syndrome[k]; + } + } + } + if (feedback_flat.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) { + std::size_t fb_weight = 0; + for (std::size_t k = 0; k < numCols; ++k) { + if (feedback_flat[j * numCols + k] != 0) + fb_weight++; + } + std::vector det(2 + fb_weight); + det[0] = final_syndrome[j]; + det[1] = syndrome[j]; + std::size_t det_idx = 2; + for (std::size_t k = 0; k < numCols; ++k) { + if (feedback_flat[j * numCols + k] != 0) + det[det_idx++] = final_syndrome[k]; + } + cudaq::detector(det); + } + } final_syndrome = syndrome; } @@ -58,13 +127,27 @@ __qpu__ void memory_circuit(const code::stabilizer_round &stabilizer_round, 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]; + if (obs_feedback_flat.size() == 0) { + 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); + } else { + // Feedback records from every round, then the data-qubit support. + std::size_t fb_count = num_rounds * obs_fb_weight[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[obs_fb_offset[obs] + 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]; + } + cudaq::logical_observable(obs_support); } - cudaq::logical_observable(obs_support); } // For each stabilizer, form detectors from data qubit readout connected with @@ -82,7 +165,17 @@ __qpu__ void memory_circuit(const code::stabilizer_round &stabilizer_round, } } - std::vector support(support_weight + 1); + // With inlined feedback, the boundary detector for record + // (fixed_offset + x) is extended with the declared last-round records. + std::size_t fb_weight = 0; + if (feedback_flat.size() > 0) { + for (std::size_t k = 0; k < numCols; ++k) { + if (feedback_flat[(fixed_offset + x) * numCols + k] != 0) + fb_weight++; + } + } + + std::vector support(support_weight + 1 + fb_weight); support[0] = final_syndrome[fixed_offset + x]; std::size_t support_idx = 1; for (std::size_t q = 0; q < num_data; ++q) { @@ -90,6 +183,12 @@ __qpu__ void memory_circuit(const code::stabilizer_round &stabilizer_round, support[support_idx++] = data_results[q]; } } + if (feedback_flat.size() > 0) { + for (std::size_t k = 0; k < numCols; ++k) { + if (feedback_flat[(fixed_offset + x) * numCols + k] != 0) + support[support_idx++] = final_syndrome[k]; + } + } cudaq::detector(support); } diff --git a/libs/qec/lib/device/memory_circuit.h b/libs/qec/lib/device/memory_circuit.h index 1f0280d9a..a15b25842 100644 --- a/libs/qec/lib/device/memory_circuit.h +++ b/libs/qec/lib/device/memory_circuit.h @@ -29,6 +29,16 @@ 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 +/// @param feedback_flat Row-major flattened inlined feedback matrix for +/// detectors. Size is 0 (no feedback; legacy detector structure) or +/// numCols*numCols with numCols = numAncx + numAncz. Entry (j, k) = 1 +/// means the cross-round detector for record j additionally XORs +/// record k of the earlier round, and the final boundary detector for +/// record j additionally XORs record k of the last round. +/// @param obs_feedback_flat Row-major flattened inlined feedback matrix for +/// logical observables. Size is 0 (no feedback) or +/// num_observables*numCols. Entry (m, k) = 1 means logical observable +/// m additionally XORs record k of every round. __qpu__ void memory_circuit(const code::stabilizer_round &stabilizer_round, const code::one_qubit_encoding &statePrep, std::size_t numData, std::size_t numAncx, @@ -37,5 +47,7 @@ __qpu__ void memory_circuit(const code::stabilizer_round &stabilizer_round, const std::vector &z_stabilizers, const std::vector &obs_matrix_flat, std::size_t num_observables, - bool measure_in_x_basis); + bool measure_in_x_basis, + const std::vector &feedback_flat, + const std::vector &obs_feedback_flat); } // namespace cudaq::qec diff --git a/libs/qec/lib/experiments.cpp b/libs/qec/lib/experiments.cpp index 9dbb9d0d1..efb9add83 100644 --- a/libs/qec/lib/experiments.cpp +++ b/libs/qec/lib/experiments.cpp @@ -254,6 +254,29 @@ stabilizer_detector_ranges(std::size_t numRounds, std::size_t numZStabs, return ranges; } +/// @brief Validate and flatten (row-major) an inlined feedback tensor +/// declared by a code. An empty tensor (the identity default) flattens to an +/// empty vector; a non-empty tensor must have shape +/// [expected_rows, expected_cols]. +static std::vector +flatten_feedback_tensor(const cudaqx::tensor &t, + std::size_t expected_rows, std::size_t expected_cols, + const std::string &name) { + if (t.size() == 0) + return {}; + if (t.rank() != 2 || t.shape()[0] != expected_rows || + t.shape()[1] != expected_cols) { + std::string actual; + for (std::size_t d = 0; d < t.rank(); ++d) + actual += (d ? ", " : "") + std::to_string(t.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."); + } + return std::vector(t.data(), t.data() + t.size()); +} + /// @brief Given a memory circuit setup, sample syndrome and data qubit /// measurements. This is the main driver function that all of the /// sample_memory_circuit overloads invoke. @@ -309,6 +332,12 @@ sample_memory_circuit(const code &code, operation statePrep, const std::size_t numCols = numAncx + numAncz; + auto feedback_flat = flatten_feedback_tensor( + code.get_inlined_feedback(), numCols, numCols, "get_inlined_feedback()"); + auto obs_feedback_flat = + flatten_feedback_tensor(code.get_observable_inlined_feedback(), num_obs, + numCols, "get_observable_inlined_feedback()"); + // Obtain the Measurement-to-Detector (M2D) sparse matrix. // m2d.rows[d] = set of chronological measurement indices whose XOR = detector // d. @@ -316,14 +345,16 @@ sample_memory_circuit(const code &code, operation statePrep, 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); + xVec, zVec, obs_flat, num_obs, !is_z_prep, + feedback_flat, obs_feedback_flat); // 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_flat, obs_feedback_flat); // mzTable[shot, meas_idx] = raw 0/1 outcome; shape (numShots, // numMeasPerShot). Measurement layout per shot: numRounds*numCols ancilla, @@ -479,11 +510,19 @@ 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_flat = flatten_feedback_tensor( + code.get_inlined_feedback(), numAnc, numAnc, "get_inlined_feedback()"); + auto obs_feedback_flat = + flatten_feedback_tensor(code.get_observable_inlined_feedback(), num_obs, + numAnc, "get_observable_inlined_feedback()"); + 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_flat, obs_feedback_flat); 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/python/bindings/py_code.cpp b/libs/qec/python/bindings/py_code.cpp index 8e3c493d7..3c63816b6 100644 --- a/libs/qec/python/bindings/py_code.cpp +++ b/libs/qec/python/bindings/py_code.cpp @@ -29,11 +29,52 @@ using namespace cudaqx; namespace cudaq::qec { +/// @brief Convert the array returned by a Python inlined-feedback getter to +/// an owning cudaqx::tensor. +/// @details The result is coerced to a C-contiguous uint8 numpy array first, +/// so any integer-convertible input (other dtypes, nested lists) is accepted. +/// 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 coerced = np.attr("ascontiguousarray")(array, np.attr("uint8")); + auto ndarray = nb::cast>(coerced); + if (ndarray.ndim() != 2) + throw std::runtime_error(std::string(methodName) + + " must return a 2-D array."); + 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, 8); 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() const override { + nb::detail::ticket nb_ticket(nb_trampoline, + "get_observable_inlined_feedback", false); + if (nb_ticket.key.is_valid()) + return feedbackArrayToTensor(nb_trampoline.base().attr(nb_ticket.key)(), + "get_observable_inlined_feedback"); + return NBBase::get_observable_inlined_feedback(); + } + // Trampoline methods for pure virtual functions std::size_t get_num_data_qubits() const override { NB_OVERRIDE_PURE(get_num_data_qubits); @@ -186,6 +227,24 @@ class PyCodeHandle : public qec::code { return m_py_operation_encodings; } + /// @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() const override { + if (!nb::hasattr(pyCode, "get_observable_inlined_feedback")) + return code::get_observable_inlined_feedback(); + return feedbackArrayToTensor( + pyCode.attr("get_observable_inlined_feedback")(), + "get_observable_inlined_feedback"); + } + protected: // Trampoline methods for pure virtual functions std::size_t get_num_data_qubits() const override { diff --git a/libs/qec/python/tests/test_code.py b/libs/qec/python/tests/test_code.py index 0280d037a..ca8a0f3d1 100644 --- a/libs/qec/python/tests/test_code.py +++ b/libs/qec/python/tests/test_code.py @@ -315,5 +315,157 @@ 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 test_python_inlined_feedback_toy(): + + @qec.code('py-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 + + def get_inlined_feedback(self): + return np.array([[0, 1], [0, 0]], dtype=np.uint8) + + def get_observable_inlined_feedback(self): + # Exercise the bridge's dtype coercion: int64 instead of uint8. + return np.array([[0, 1]]) + + # Sample on stim, mirroring the C++ unit tests (test_qec links the stim + # target). The toy's X record is genuinely random in round 1, and the + # default (qpp) target does not preserve cross-round mid-circuit + # measurement correlations when sampling Python kernels - independent of + # the inlined-feedback machinery, but this toy (unlike the steane + # example, whose records are all deterministically 0) exposes it. + cudaq.set_target('stim') + + numShots, numRounds = 20, 4 + syndromes, data = qec.sample_memory_circuit(qec.get_code('py-feedback-toy'), + 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) + + # The individual data qubits are randomized by the XX measurement, but + # the ZZ parity must close deterministically every shot. + assert data.shape == (numShots, 2) + assert not np.any(data[:, 0] ^ data[:, 1]) + 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. + @qec.code('py-feedback-toy-nofb') + class FeedbackToyNoFeedback: + + 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 + + # 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) + + if __name__ == "__main__": pytest.main() diff --git a/libs/qec/unittests/CMakeLists.txt b/libs/qec/unittests/CMakeLists.txt index 344d05449..588834825 100644 --- a/libs/qec/unittests/CMakeLists.txt +++ b/libs/qec/unittests/CMakeLists.txt @@ -50,7 +50,16 @@ add_dependencies(CUDAQXQECUnitTests test_decoders_yaml) gtest_discover_tests(test_decoders_yaml) 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..50b5e2a03 --- /dev/null +++ b/libs/qec/unittests/feedback_toy_device.cpp @@ -0,0 +1,56 @@ +/****************************************************************-*- 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]); +} + +// A superdense (Bell-pair) round for a 2-data-qubit code with stabilizers XX +// and ZZ that deliberately 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. 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 byproduct: X^{r_X} on data[1]. + cudaq::x(p.ancx[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 47326f084..0358784ed 100644 --- a/libs/qec/unittests/test_qec.cpp +++ b/libs/qec/unittests/test_qec.cpp @@ -2192,3 +2192,167 @@ 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 = steane->get_observable_inlined_feedback(); + EXPECT_EQ(obs_feedback.rank(), 0); + EXPECT_EQ(obs_feedback.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__ 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 observable 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 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. +// +// Derivation of the feedback matrices (records per round in [Z][X] order: +// record 0 = ancz outcome r_Z, record 1 = ancx outcome r_X): +// - r_Z(t) reads the ZZ parity of the data entering round t, i.e. the +// accumulated byproduct parity b_t, with b_1 = 0 and +// b_{t+1} = b_t XOR r_X(t). +// - r_X(t) reads XX: uniformly random in round 1 of a Z-basis memory and +// constant afterwards (the X byproduct commutes with XX). +// - Cross-round detector for record 0 without feedback: +// r_Z(t) XOR r_Z(t+1) = r_X(t) -> non-deterministic (the negative-control +// test relies on stim rejecting exactly this). XOR-ing in the earlier +// round's record 1 fixes it: feedback(0, 1) = 1. +// - Record 1 compares deterministically on its own: feedback row 1 = 0. +// - Final boundary detector for record 0: +// r_Z(T) XOR d0 XOR d1 = b_T XOR b_{T+1} = r_X(T); the same +// feedback(0, 1) entry extends the boundary detector with the last +// round's record 1. +// - The ZZ observable closes as d0 XOR d1 = b_{T+1} = XOR of every round's +// record 1: observable_feedback(0, 1) = 1. +// Note that the byproduct must be heralded by the *other* (X) record: a +// byproduct conditioned on the Z record - which is identically 0 in the +// noiseless circuit - would never fire and could not break determinism. +// Both matrices were validated against stim's determinism analysis for +// 1 through 6 rounds. +// --------------------------------------------------------------------------- + +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; + +public: + explicit feedback_toy_code(bool declare_feedback) + : declare_feedback(declare_feedback) { + operation_encodings.insert(std::make_pair( + cudaq::qec::operation::prep0, cudaq::qec::feedback_toy::prep0)); + 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({"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}) = 1; + return feedback; + } + + cudaqx::tensor get_observable_inlined_feedback() const override { + if (!declare_feedback) + return code::get_observable_inlined_feedback(); + cudaqx::tensor feedback({1, 2}); + feedback.at({0, 1}) = 1; + return feedback; + } +}; + +} // namespace + +TEST(QECCodeTester, checkInlinedFeedbackToyMemoryCircuit) { + feedback_toy_code toy(/*declare_feedback=*/true); + 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); + + // 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) + << "shot " << shot << ", detector " << d; + + // The individual data qubits are randomized by the XX measurement, but the + // ZZ parity must close deterministically every shot. + ASSERT_EQ(data.rank(), 2); + ASSERT_EQ(data.shape()[0], numShots); + ASSERT_EQ(data.shape()[1], 2); + for (std::size_t shot = 0; shot < numShots; ++shot) + EXPECT_EQ(data.at({shot, 0}) ^ data.at({shot, 1}), 0) << "shot " << shot; +} + +TEST(QECCodeTester, checkInlinedFeedbackToyNegativeControl) { + // The identical toy without the feedback declaration: the uncorrected + // byproduct makes the record-0 detectors non-deterministic and stim must + // reject the circuit. This proves the declared feedback is doing the work + // in the positive test above. + feedback_toy_code toy(/*declare_feedback=*/false); + try { + cudaq::qec::sample_memory_circuit( + toy, cudaq::qec::operation::prep0, /*numShots=*/20, /*numRounds=*/4); + FAIL() << "expected sample_memory_circuit to reject the feedback-less toy"; + } catch (const std::exception &e) { + // Fail-for-the-right-reason check: stim's determinism analysis must be + // what rejects the circuit (it flags the record-0 cross-round and final + // boundary detectors). + 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); + 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=*/4, noise); + EXPECT_GT(dem.num_detectors(), 0); + EXPECT_GT(dem.num_error_mechanisms(), 0); + EXPECT_EQ(dem.num_observables(), 1); +} From 8bffc7548d4819a3da45f9dde7f58b1050bf30bf Mon Sep 17 00:00:00 2001 From: kvmto Date: Sat, 11 Jul 2026 01:17:28 +0200 Subject: [PATCH 2/4] Extract inlined-feedback composition rule into host-side apply_inlined_feedback The record->detector composition rule now lives in detector_error_model as a CSR layout builder; the memory_circuit kernel is a thin emitter walking the layout. Behavior is bit-identical; all suites pass unchanged. Signed-off-by: kvmto --- .../include/cudaq/qec/detector_error_model.h | 43 +++++++ libs/qec/lib/detector_error_model.cpp | 45 +++++++ libs/qec/lib/device/memory_circuit.cpp | 111 +++++++----------- libs/qec/lib/device/memory_circuit.h | 30 +++-- libs/qec/lib/experiments.cpp | 60 +++------- libs/qec/python/tests/test_code.py | 38 +++--- libs/qec/unittests/test_qec.cpp | 69 ++++++++++- 7 files changed, 254 insertions(+), 142 deletions(-) diff --git a/libs/qec/include/cudaq/qec/detector_error_model.h b/libs/qec/include/cudaq/qec/detector_error_model.h index 64729f3d0..3899dd2f7 100644 --- a/libs/qec/include/cudaq/qec/detector_error_model.h +++ b/libs/qec/include/cudaq/qec/detector_error_model.h @@ -82,6 +82,49 @@ struct detector_error_model { bool remove_zero_syndrome_errors = false); }; +/// @brief Record-to-detector composition layout derived from a code's +/// declared inlined-feedback matrices (CSR-style). +/// +/// Entry range [detector_offsets[j], detector_offsets[j+1]) of +/// detector_indices lists the herald record columns XOR-ed into record j's +/// cross-round detector (records of the earlier round) and into its final +/// boundary detector (records of the last round). Entry range +/// [observable_offsets[m], observable_offsets[m+1]) of observable_indices +/// lists the record columns XOR-ed into logical observable m on every round. +/// Empty offsets vectors mean no feedback is declared for that target +/// (identity detector/observable structure). +struct inlined_feedback_layout { + std::vector detector_indices; + std::vector detector_offsets; + std::vector observable_indices; + std::vector observable_offsets; +}; + +/// @brief Build the record-to-detector composition layout from a code's +/// declared inlined-feedback matrices. +/// +/// @param feedback Detector feedback matrix, shape +/// [num_syndromes_per_round x num_syndromes_per_round] with 0/1 +/// entries, or an empty tensor for none. Entry (j, k) = 1 means the +/// cross-round detector for record j additionally XORs record k of +/// the earlier round, and the final boundary detector for record j +/// additionally XORs record k of the last round. +/// @param obs_feedback Observable feedback matrix, shape +/// [num_observables x num_syndromes_per_round] with 0/1 entries, or +/// an empty tensor for none. Entry (m, k) = 1 means logical +/// observable m additionally XORs record k of every round. +/// @param num_syndromes_per_round Records per round (number of ancilla +/// qubits). +/// @param num_observables Number of logical observables. +/// @return The CSR layout; see inlined_feedback_layout. +/// @throws std::runtime_error if a non-empty tensor's shape does not match +/// the expected dimensions. +inlined_feedback_layout +apply_inlined_feedback(const cudaqx::tensor &feedback, + const cudaqx::tensor &obs_feedback, + std::size_t num_syndromes_per_round, + std::size_t num_observables); + /// Parse the Stim DEM string @p dem_text into detector/observable flip /// matrices and error rates. DEM-native decoders should consume raw DEM text /// instead. By default (@p use_decomp_suggestions = false) the '^' separators diff --git a/libs/qec/lib/detector_error_model.cpp b/libs/qec/lib/detector_error_model.cpp index 79a89b798..aa1c683fc 100644 --- a/libs/qec/lib/detector_error_model.cpp +++ b/libs/qec/lib/detector_error_model.cpp @@ -296,4 +296,49 @@ void detector_error_model::canonicalize_for_rounds( this->observables_flips_matrix, final_column_order); } +/// @brief Convert one 0/1 feedback matrix into a CSR row layout. An empty +/// tensor (the identity default) leaves @p indices and @p offsets empty; a +/// non-empty tensor must have shape [num_rows, num_cols]. +static void feedback_rows_to_csr(const cudaqx::tensor &matrix, + std::size_t num_rows, std::size_t num_cols, + const char *name, + std::vector &indices, + std::vector &offsets) { + if (matrix.rank() == 0 || matrix.size() == 0) + return; + if (matrix.rank() != 2 || matrix.shape()[0] != num_rows || + matrix.shape()[1] != num_cols) { + std::string actual; + for (std::size_t d = 0; d < matrix.rank(); ++d) + actual += (d ? ", " : "") + std::to_string(matrix.shape()[d]); + throw std::runtime_error( + std::string(name) + " has invalid shape [" + actual + "] - expected [" + + std::to_string(num_rows) + ", " + std::to_string(num_cols) + + "] or an empty tensor."); + } + offsets.resize(num_rows + 1); + offsets[0] = 0; + for (std::size_t r = 0; r < num_rows; ++r) { + for (std::size_t c = 0; c < num_cols; ++c) + if (matrix.at({r, c}) != 0) + indices.push_back(c); + offsets[r + 1] = indices.size(); + } +} + +inlined_feedback_layout +apply_inlined_feedback(const cudaqx::tensor &feedback, + const cudaqx::tensor &obs_feedback, + std::size_t num_syndromes_per_round, + std::size_t num_observables) { + inlined_feedback_layout layout; + feedback_rows_to_csr(feedback, num_syndromes_per_round, + num_syndromes_per_round, "inlined feedback", + layout.detector_indices, layout.detector_offsets); + feedback_rows_to_csr(obs_feedback, num_observables, num_syndromes_per_round, + "observable inlined feedback", layout.observable_indices, + layout.observable_offsets); + return layout; +} + } // namespace cudaq::qec diff --git a/libs/qec/lib/device/memory_circuit.cpp b/libs/qec/lib/device/memory_circuit.cpp index b2d8e9ae3..0fa0db3c2 100644 --- a/libs/qec/lib/device/memory_circuit.cpp +++ b/libs/qec/lib/device/memory_circuit.cpp @@ -19,8 +19,10 @@ __qpu__ void memory_circuit(const code::stabilizer_round &stabilizer_round, const std::vector &obs_matrix_flat, std::size_t num_observables, bool measure_in_x_basis, - const std::vector &feedback_flat, - const std::vector &obs_feedback_flat) { + const std::vector &fb_indices, + const std::vector &fb_offsets, + const std::vector &obs_fb_indices, + const std::vector &obs_fb_offsets) { // Allocate the data and ancilla qubits cudaq::qvector data(num_data), xstab_anc(numAncx), zstab_anc(numAncz); @@ -43,72 +45,51 @@ __qpu__ void memory_circuit(const code::stabilizer_round &stabilizer_round, 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 prefix - // offsets: observable m's feedback records occupy - // [obs_fb_offset[m], obs_fb_offset[m] + num_rounds * obs_fb_weight[m]), - // laid out round-major within each observable's slice. - std::vector obs_fb_weight(num_observables); - std::vector obs_fb_offset(num_observables); - std::size_t obs_fb_total = 0; - if (obs_feedback_flat.size() > 0) { - for (std::size_t m = 0; m < num_observables; ++m) { - std::size_t weight = 0; - for (std::size_t k = 0; k < numCols; ++k) { - if (obs_feedback_flat[m * numCols + k] != 0) - weight++; - } - obs_fb_weight[m] = weight; - obs_fb_offset[m] = obs_fb_total; - obs_fb_total += num_rounds * weight; - } - } + // Observable inlined-feedback accumulation. Nested std::vector is not + // supported in __qpu__ code, so observable m's per-round feedback records + // occupy the round-major slice starting at num_rounds * obs_fb_offsets[m] + // of a single flat buffer; the round-r block of observable m starts at + // num_rounds * obs_fb_offsets[m] + r * row_weight(m) with + // row_weight(m) = obs_fb_offsets[m + 1] - obs_fb_offsets[m]. + bool has_obs_fb = obs_fb_offsets.size() > 0; + std::size_t obs_fb_total = + has_obs_fb ? num_rounds * obs_fb_offsets[num_observables] : 0; std::vector obs_fb_records(obs_fb_total); // Collect the first-round feedback records for each observable. - if (obs_feedback_flat.size() > 0) { + if (has_obs_fb) { for (std::size_t m = 0; m < num_observables; ++m) { - std::size_t idx = obs_fb_offset[m]; - for (std::size_t k = 0; k < numCols; ++k) { - if (obs_feedback_flat[m * numCols + k] != 0) - obs_fb_records[idx++] = final_syndrome[k]; - } + std::size_t idx = num_rounds * obs_fb_offsets[m]; + for (std::size_t i = obs_fb_offsets[m]; i < obs_fb_offsets[m + 1]; ++i) + obs_fb_records[idx++] = final_syndrome[obs_fb_indices[i]]; } } // Generate syndrome data for (std::size_t round = 1; round < num_rounds; ++round) { auto syndrome = stabilizer_round(logical, x_stabilizers, z_stabilizers); - if (obs_feedback_flat.size() > 0) { + if (has_obs_fb) { for (std::size_t m = 0; m < num_observables; ++m) { - std::size_t idx = obs_fb_offset[m] + round * obs_fb_weight[m]; - for (std::size_t k = 0; k < numCols; ++k) { - if (obs_feedback_flat[m * numCols + k] != 0) - obs_fb_records[idx++] = syndrome[k]; - } + std::size_t row_weight = obs_fb_offsets[m + 1] - obs_fb_offsets[m]; + std::size_t idx = num_rounds * obs_fb_offsets[m] + round * row_weight; + for (std::size_t i = obs_fb_offsets[m]; i < obs_fb_offsets[m + 1]; ++i) + obs_fb_records[idx++] = syndrome[obs_fb_indices[i]]; } } - if (feedback_flat.size() == 0) { + if (fb_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. + // augmented with the earlier-round herald records in row j of the + // feedback layout. for (std::size_t j = 0; j < numCols; ++j) { - std::size_t fb_weight = 0; - for (std::size_t k = 0; k < numCols; ++k) { - if (feedback_flat[j * numCols + k] != 0) - fb_weight++; - } - std::vector det(2 + fb_weight); + std::size_t begin = fb_offsets[j]; + std::size_t end = fb_offsets[j + 1]; + std::vector det(2 + end - begin); det[0] = final_syndrome[j]; det[1] = syndrome[j]; - std::size_t det_idx = 2; - for (std::size_t k = 0; k < numCols; ++k) { - if (feedback_flat[j * numCols + k] != 0) - det[det_idx++] = final_syndrome[k]; - } + for (std::size_t i = begin; i < end; ++i) + det[2 + (i - begin)] = final_syndrome[fb_indices[i]]; cudaq::detector(det); } } @@ -127,7 +108,7 @@ __qpu__ void memory_circuit(const code::stabilizer_round &stabilizer_round, if (obs_matrix_flat[obs * num_data + q] != 0) support_weight++; } - if (obs_feedback_flat.size() == 0) { + if (!has_obs_fb) { std::vector obs_support(support_weight); std::size_t idx = 0; for (std::size_t q = 0; q < num_data; ++q) { @@ -137,10 +118,12 @@ __qpu__ void memory_circuit(const code::stabilizer_round &stabilizer_round, cudaq::logical_observable(obs_support); } else { // Feedback records from every round, then the data-qubit support. - std::size_t fb_count = num_rounds * obs_fb_weight[obs]; + std::size_t row_weight = obs_fb_offsets[obs + 1] - obs_fb_offsets[obs]; + std::size_t fb_count = num_rounds * row_weight; + std::size_t base = num_rounds * obs_fb_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[obs_fb_offset[obs] + 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) @@ -166,16 +149,16 @@ __qpu__ void memory_circuit(const code::stabilizer_round &stabilizer_round, } // With inlined feedback, the boundary detector for record - // (fixed_offset + x) is extended with the declared last-round records. - std::size_t fb_weight = 0; - if (feedback_flat.size() > 0) { - for (std::size_t k = 0; k < numCols; ++k) { - if (feedback_flat[(fixed_offset + x) * numCols + k] != 0) - fb_weight++; - } + // (fixed_offset + x) is extended with its declared last-round records. + std::size_t fb_begin = 0; + std::size_t fb_end = 0; + if (fb_offsets.size() > 0) { + fb_begin = fb_offsets[fixed_offset + x]; + fb_end = fb_offsets[fixed_offset + x + 1]; } - std::vector support(support_weight + 1 + fb_weight); + std::vector support(support_weight + 1 + + (fb_end - fb_begin)); support[0] = final_syndrome[fixed_offset + x]; std::size_t support_idx = 1; for (std::size_t q = 0; q < num_data; ++q) { @@ -183,12 +166,8 @@ __qpu__ void memory_circuit(const code::stabilizer_round &stabilizer_round, support[support_idx++] = data_results[q]; } } - if (feedback_flat.size() > 0) { - for (std::size_t k = 0; k < numCols; ++k) { - if (feedback_flat[(fixed_offset + x) * numCols + k] != 0) - support[support_idx++] = final_syndrome[k]; - } - } + for (std::size_t i = fb_begin; i < fb_end; ++i) + support[support_idx++] = final_syndrome[fb_indices[i]]; cudaq::detector(support); } diff --git a/libs/qec/lib/device/memory_circuit.h b/libs/qec/lib/device/memory_circuit.h index a15b25842..e21ac943b 100644 --- a/libs/qec/lib/device/memory_circuit.h +++ b/libs/qec/lib/device/memory_circuit.h @@ -29,16 +29,20 @@ 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 -/// @param feedback_flat Row-major flattened inlined feedback matrix for -/// detectors. Size is 0 (no feedback; legacy detector structure) or -/// numCols*numCols with numCols = numAncx + numAncz. Entry (j, k) = 1 -/// means the cross-round detector for record j additionally XORs -/// record k of the earlier round, and the final boundary detector for -/// record j additionally XORs record k of the last round. -/// @param obs_feedback_flat Row-major flattened inlined feedback matrix for -/// logical observables. Size is 0 (no feedback) or -/// num_observables*numCols. Entry (m, k) = 1 means logical observable -/// m additionally XORs record k of every round. +/// @param fb_indices CSR herald-column indices for the detector inlined +/// feedback layout. `fb_indices[fb_offsets[j] .. fb_offsets[j+1])` +/// lists the herald record columns XORed into record j's cross-round +/// detector (records of the earlier round) and into its final boundary +/// detector (records of the last round). +/// @param fb_offsets CSR row offsets for the detector inlined feedback layout. +/// Size is 0 (no feedback; legacy `cudaq::detectors` structure) or +/// numCols + 1 with numCols = numAncx + numAncz. +/// @param obs_fb_indices CSR record-column indices for the observable inlined +/// feedback layout. `obs_fb_indices[obs_fb_offsets[m] .. +/// obs_fb_offsets[m+1])` lists the record columns XORed into logical +/// observable m on every round. +/// @param obs_fb_offsets CSR row offsets for the observable inlined feedback +/// layout. 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, @@ -48,6 +52,8 @@ __qpu__ void memory_circuit(const code::stabilizer_round &stabilizer_round, const std::vector &obs_matrix_flat, std::size_t num_observables, bool measure_in_x_basis, - const std::vector &feedback_flat, - const std::vector &obs_feedback_flat); + const std::vector &fb_indices, + const std::vector &fb_offsets, + const std::vector &obs_fb_indices, + const std::vector &obs_fb_offsets); } // namespace cudaq::qec diff --git a/libs/qec/lib/experiments.cpp b/libs/qec/lib/experiments.cpp index efb9add83..39a465475 100644 --- a/libs/qec/lib/experiments.cpp +++ b/libs/qec/lib/experiments.cpp @@ -254,29 +254,6 @@ stabilizer_detector_ranges(std::size_t numRounds, std::size_t numZStabs, return ranges; } -/// @brief Validate and flatten (row-major) an inlined feedback tensor -/// declared by a code. An empty tensor (the identity default) flattens to an -/// empty vector; a non-empty tensor must have shape -/// [expected_rows, expected_cols]. -static std::vector -flatten_feedback_tensor(const cudaqx::tensor &t, - std::size_t expected_rows, std::size_t expected_cols, - const std::string &name) { - if (t.size() == 0) - return {}; - if (t.rank() != 2 || t.shape()[0] != expected_rows || - t.shape()[1] != expected_cols) { - std::string actual; - for (std::size_t d = 0; d < t.rank(); ++d) - actual += (d ? ", " : "") + std::to_string(t.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."); - } - return std::vector(t.data(), t.data() + t.size()); -} - /// @brief Given a memory circuit setup, sample syndrome and data qubit /// measurements. This is the main driver function that all of the /// sample_memory_circuit overloads invoke. @@ -332,29 +309,29 @@ sample_memory_circuit(const code &code, operation statePrep, const std::size_t numCols = numAncx + numAncz; - auto feedback_flat = flatten_feedback_tensor( - code.get_inlined_feedback(), numCols, numCols, "get_inlined_feedback()"); - auto obs_feedback_flat = - flatten_feedback_tensor(code.get_observable_inlined_feedback(), num_obs, - numCols, "get_observable_inlined_feedback()"); + auto fb_layout = apply_inlined_feedback( + code.get_inlined_feedback(), code.get_observable_inlined_feedback(), + numCols, num_obs); // 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, - feedback_flat, obs_feedback_flat); + 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, fb_layout.detector_indices, fb_layout.detector_offsets, + fb_layout.observable_indices, fb_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, feedback_flat, obs_feedback_flat); + auto result = cudaq::sample( + opts, memory_circuit, stabRound, prep, numData, numAncx, numAncz, + numRounds, xVec, zVec, obs_flat, num_obs, !is_z_prep, + fb_layout.detector_indices, fb_layout.detector_offsets, + fb_layout.observable_indices, fb_layout.observable_offsets); // mzTable[shot, meas_idx] = raw 0/1 outcome; shape (numShots, // numMeasPerShot). Measurement layout per shot: numRounds*numCols ancilla, @@ -511,18 +488,17 @@ dem_from_memory_circuit(const code &code, operation statePrep, logical_obs.data() + logical_obs.size()); const std::size_t numAnc = numAncx + numAncz; - auto feedback_flat = flatten_feedback_tensor( - code.get_inlined_feedback(), numAnc, numAnc, "get_inlined_feedback()"); - auto obs_feedback_flat = - flatten_feedback_tensor(code.get_observable_inlined_feedback(), num_obs, - numAnc, "get_observable_inlined_feedback()"); + auto fb_layout = apply_inlined_feedback( + code.get_inlined_feedback(), code.get_observable_inlined_feedback(), + numAnc, num_obs); 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, - feedback_flat, obs_feedback_flat); + fb_layout.detector_indices, fb_layout.detector_offsets, + fb_layout.observable_indices, fb_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/python/tests/test_code.py b/libs/qec/python/tests/test_code.py index ca8a0f3d1..0ba39a425 100644 --- a/libs/qec/python/tests/test_code.py +++ b/libs/qec/python/tests/test_code.py @@ -403,22 +403,28 @@ def get_observable_inlined_feedback(self): # the inlined-feedback machinery, but this toy (unlike the steane # example, whose records are all deterministically 0) exposes it. cudaq.set_target('stim') - - numShots, numRounds = 20, 4 - syndromes, data = qec.sample_memory_circuit(qec.get_code('py-feedback-toy'), - 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) - - # The individual data qubits are randomized by the XX measurement, but - # the ZZ parity must close deterministically every shot. - assert data.shape == (numShots, 2) - assert not np.any(data[:, 0] ^ data[:, 1]) - cudaq.reset_target() + try: + # numRounds >= 2 exercises the final boundary detector (same-round + # herald) distinctly from the cross-round detectors (earlier-round + # herald). + numShots, numRounds = 20, 4 + syndromes, data = qec.sample_memory_circuit( + qec.get_code('py-feedback-toy'), + 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) + + # The individual data qubits are randomized by the XX measurement, + # but the ZZ parity must close deterministically every shot. + assert data.shape == (numShots, 2) + assert not np.any(data[:, 0] ^ data[:, 1]) + finally: + cudaq.reset_target() def test_python_inlined_feedback_toy_negative_control(): diff --git a/libs/qec/unittests/test_qec.cpp b/libs/qec/unittests/test_qec.cpp index 0358784ed..429f960b3 100644 --- a/libs/qec/unittests/test_qec.cpp +++ b/libs/qec/unittests/test_qec.cpp @@ -2268,8 +2268,8 @@ class feedback_toy_code : public cudaq::qec::code { public: explicit feedback_toy_code(bool declare_feedback) : declare_feedback(declare_feedback) { - operation_encodings.insert(std::make_pair( - cudaq::qec::operation::prep0, cudaq::qec::feedback_toy::prep0)); + operation_encodings.insert(std::make_pair(cudaq::qec::operation::prep0, + cudaq::qec::feedback_toy::prep0)); operation_encodings.insert( std::make_pair(cudaq::qec::operation::stabilizer_round, cudaq::qec::feedback_toy::stabilizer_round)); @@ -2329,8 +2329,8 @@ TEST(QECCodeTester, checkInlinedFeedbackToyNegativeControl) { // in the positive test above. feedback_toy_code toy(/*declare_feedback=*/false); try { - cudaq::qec::sample_memory_circuit( - toy, cudaq::qec::operation::prep0, /*numShots=*/20, /*numRounds=*/4); + cudaq::qec::sample_memory_circuit(toy, cudaq::qec::operation::prep0, + /*numShots=*/20, /*numRounds=*/4); FAIL() << "expected sample_memory_circuit to reject the feedback-less toy"; } catch (const std::exception &e) { // Fail-for-the-right-reason check: stim's determinism analysis must be @@ -2347,12 +2347,69 @@ TEST(QECCodeTester, checkInlinedFeedbackToyDem) { // 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 = 4; 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=*/4, noise); - EXPECT_GT(dem.num_detectors(), 0); + 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::apply_inlined_feedback(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::apply_inlined_feedback(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::apply_inlined_feedback(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::apply_inlined_feedback(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::apply_inlined_feedback(bad, empty, 2, 0), + std::runtime_error); + cudaqx::tensor bad_obs({2, 2}); // expected [1 x 2] + EXPECT_THROW(cudaq::qec::apply_inlined_feedback(empty, bad_obs, 2, 1), + std::runtime_error); +} From 8dbe2daa183e7efb99650e7533c5f96ea11545c0 Mon Sep 17 00:00:00 2001 From: kvmto Date: Sat, 11 Jul 2026 02:20:43 +0200 Subject: [PATCH 3/4] Restore data-shaped memory_circuit API; extract feedback logic into named device helpers The kernel keeps its flat-matrix signature; all record-composition logic lives in reusable __qpu__ functions in device/inlined_feedback.h. The host CSR layout builder (build_inlined_feedback_layout) is retained for record-stream consumers and pinned to the kernel by an M2D conformance test Signed-off-by: kvmto --- .../include/cudaq/qec/detector_error_model.h | 8 +- libs/qec/lib/detector_error_model.cpp | 8 +- libs/qec/lib/device/inlined_feedback.h | 161 ++++++++++++++++++ libs/qec/lib/device/memory_circuit.cpp | 143 ++++------------ libs/qec/lib/device/memory_circuit.h | 30 ++-- libs/qec/lib/experiments.cpp | 60 +++++-- libs/qec/unittests/test_qec.cpp | 131 +++++++++++++- 7 files changed, 385 insertions(+), 156 deletions(-) create mode 100644 libs/qec/lib/device/inlined_feedback.h diff --git a/libs/qec/include/cudaq/qec/detector_error_model.h b/libs/qec/include/cudaq/qec/detector_error_model.h index 3899dd2f7..88b053303 100644 --- a/libs/qec/include/cudaq/qec/detector_error_model.h +++ b/libs/qec/include/cudaq/qec/detector_error_model.h @@ -120,10 +120,10 @@ struct inlined_feedback_layout { /// @throws std::runtime_error if a non-empty tensor's shape does not match /// the expected dimensions. inlined_feedback_layout -apply_inlined_feedback(const cudaqx::tensor &feedback, - const cudaqx::tensor &obs_feedback, - std::size_t num_syndromes_per_round, - std::size_t num_observables); +build_inlined_feedback_layout(const cudaqx::tensor &feedback, + const cudaqx::tensor &obs_feedback, + std::size_t num_syndromes_per_round, + std::size_t num_observables); /// Parse the Stim DEM string @p dem_text into detector/observable flip /// matrices and error rates. DEM-native decoders should consume raw DEM text diff --git a/libs/qec/lib/detector_error_model.cpp b/libs/qec/lib/detector_error_model.cpp index aa1c683fc..338f455c3 100644 --- a/libs/qec/lib/detector_error_model.cpp +++ b/libs/qec/lib/detector_error_model.cpp @@ -327,10 +327,10 @@ static void feedback_rows_to_csr(const cudaqx::tensor &matrix, } inlined_feedback_layout -apply_inlined_feedback(const cudaqx::tensor &feedback, - const cudaqx::tensor &obs_feedback, - std::size_t num_syndromes_per_round, - std::size_t num_observables) { +build_inlined_feedback_layout(const cudaqx::tensor &feedback, + const cudaqx::tensor &obs_feedback, + std::size_t num_syndromes_per_round, + std::size_t num_observables) { inlined_feedback_layout layout; feedback_rows_to_csr(feedback, num_syndromes_per_round, num_syndromes_per_round, "inlined feedback", diff --git a/libs/qec/lib/device/inlined_feedback.h b/libs/qec/lib/device/inlined_feedback.h new file mode 100644 index 000000000..9884ba455 --- /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 `prev[k]` for each column k (ascending) with +/// `feedback_flat(j, k) != 0`. +inline __qpu__ std::vector cross_round_detector_records( + const std::vector &prev, + const std::vector &curr, std::size_t j, + const std::vector &feedback_flat, std::size_t num_cols) { + std::size_t weight = row_support_weight(feedback_flat, j, num_cols); + std::vector det(2 + weight); + det[0] = prev[j]; + det[1] = curr[j]; + std::size_t idx = 2; + for (std::size_t k = 0; k < num_cols; ++k) + if (feedback_flat[j * num_cols + k] != 0) + det[idx++] = prev[k]; + 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 +/// `last_syndrome[k]` for each column k (ascending) with +/// `feedback_flat(record_row, k) != 0`. When @p feedback_flat is empty the +/// herald part contributes nothing, so one helper serves both the +/// 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_flat, + std::size_t num_cols) { + 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_flat.size() > 0 + ? row_support_weight(feedback_flat, record_row, num_cols) + : 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_flat.size() > 0) + for (std::size_t k = 0; k < num_cols; ++k) + if (feedback_flat[record_row * num_cols + k] != 0) + support[idx++] = last_syndrome[k]; + 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 * (weight of obs_feedback_flat row m)` records. The returned +/// vector has `num_observables + 1` entries; the last entry is the total +/// buffer size. Returns an empty vector when @p obs_feedback_flat is empty +/// (no observable feedback declared). +inline __qpu__ std::vector +observable_feedback_offsets(const std::vector &obs_feedback_flat, + std::size_t num_observables, std::size_t num_cols, + std::size_t num_rounds) { + // Sized (not default) construction: the __qpu__ dialect forbids the default + // std::vector constructor. An empty obs_feedback_flat yields a size-0 vector. + std::size_t num_offsets = + obs_feedback_flat.size() > 0 ? num_observables + 1 : 0; + std::vector offsets(num_offsets); + if (obs_feedback_flat.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 * row_support_weight(obs_feedback_flat, m, num_cols); + } + 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 `offsets[m] + round * (weight of obs_feedback_flat +/// row m)` and holds `syndrome[k]` for each column k (ascending) with +/// `obs_feedback_flat(m, k) != 0`. 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 &offsets, + const std::vector &syndrome, std::size_t round, + const std::vector &obs_feedback_flat, + std::size_t num_observables, std::size_t num_cols) { + for (std::size_t m = 0; m < num_observables; ++m) { + std::size_t weight = row_support_weight(obs_feedback_flat, m, num_cols); + std::size_t idx = offsets[m] + round * weight; + for (std::size_t k = 0; k < num_cols; ++k) + if (obs_feedback_flat[m * num_cols + k] != 0) + obs_fb_records[idx++] = syndrome[k]; + } +} + +/// @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 0fa0db3c2..ad1a88117 100644 --- a/libs/qec/lib/device/memory_circuit.cpp +++ b/libs/qec/lib/device/memory_circuit.cpp @@ -6,7 +6,7 @@ * the terms of the Apache License 2.0 which accompanies this distribution. * ******************************************************************************/ #include "memory_circuit.h" -#include +#include "inlined_feedback.h" namespace cudaq::qec { @@ -19,10 +19,8 @@ __qpu__ void memory_circuit(const code::stabilizer_round &stabilizer_round, const std::vector &obs_matrix_flat, std::size_t num_observables, bool measure_in_x_basis, - const std::vector &fb_indices, - const std::vector &fb_offsets, - const std::vector &obs_fb_indices, - const std::vector &obs_fb_offsets) { + const std::vector &feedback_flat, + const std::vector &obs_feedback_flat) { // Allocate the data and ancilla qubits cudaq::qvector data(num_data), xstab_anc(numAncx), zstab_anc(numAncz); @@ -45,53 +43,39 @@ __qpu__ void memory_circuit(const code::stabilizer_round &stabilizer_round, std::size_t numCols = numAncx + numAncz; - // Observable inlined-feedback accumulation. Nested std::vector is not - // supported in __qpu__ code, so observable m's per-round feedback records - // occupy the round-major slice starting at num_rounds * obs_fb_offsets[m] - // of a single flat buffer; the round-r block of observable m starts at - // num_rounds * obs_fb_offsets[m] + r * row_weight(m) with - // row_weight(m) = obs_fb_offsets[m + 1] - obs_fb_offsets[m]. - bool has_obs_fb = obs_fb_offsets.size() > 0; + // 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_offsets / collect_observable_feedback_round in + // inlined_feedback.h. + std::vector obs_fb_offsets = observable_feedback_offsets( + obs_feedback_flat, num_observables, numCols, num_rounds); std::size_t obs_fb_total = - has_obs_fb ? num_rounds * obs_fb_offsets[num_observables] : 0; + obs_feedback_flat.size() > 0 ? obs_fb_offsets[num_observables] : 0; std::vector obs_fb_records(obs_fb_total); // Collect the first-round feedback records for each observable. - if (has_obs_fb) { - for (std::size_t m = 0; m < num_observables; ++m) { - std::size_t idx = num_rounds * obs_fb_offsets[m]; - for (std::size_t i = obs_fb_offsets[m]; i < obs_fb_offsets[m + 1]; ++i) - obs_fb_records[idx++] = final_syndrome[obs_fb_indices[i]]; - } - } + if (obs_feedback_flat.size() > 0) + collect_observable_feedback_round(obs_fb_records, obs_fb_offsets, + final_syndrome, 0, obs_feedback_flat, + num_observables, numCols); // Generate syndrome data for (std::size_t round = 1; round < num_rounds; ++round) { auto syndrome = stabilizer_round(logical, x_stabilizers, z_stabilizers); - if (has_obs_fb) { - for (std::size_t m = 0; m < num_observables; ++m) { - std::size_t row_weight = obs_fb_offsets[m + 1] - obs_fb_offsets[m]; - std::size_t idx = num_rounds * obs_fb_offsets[m] + round * row_weight; - for (std::size_t i = obs_fb_offsets[m]; i < obs_fb_offsets[m + 1]; ++i) - obs_fb_records[idx++] = syndrome[obs_fb_indices[i]]; - } - } - if (fb_offsets.size() == 0) { + if (obs_feedback_flat.size() > 0) + collect_observable_feedback_round(obs_fb_records, obs_fb_offsets, + syndrome, round, obs_feedback_flat, + num_observables, numCols); + if (feedback_flat.size() == 0) { cudaq::detectors(final_syndrome, syndrome); } else { // Cross-round detector for record j: earlier vs current round record, - // augmented with the earlier-round herald records in row j of the - // feedback layout. - for (std::size_t j = 0; j < numCols; ++j) { - std::size_t begin = fb_offsets[j]; - std::size_t end = fb_offsets[j + 1]; - std::vector det(2 + end - begin); - det[0] = final_syndrome[j]; - det[1] = syndrome[j]; - for (std::size_t i = begin; i < end; ++i) - det[2 + (i - begin)] = final_syndrome[fb_indices[i]]; - cudaq::detector(det); - } + // 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_flat, numCols)); } final_syndrome = syndrome; } @@ -101,76 +85,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++; - } - if (!has_obs_fb) { - 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); - } else { - // Feedback records from every round, then the data-qubit support. - std::size_t row_weight = obs_fb_offsets[obs + 1] - obs_fb_offsets[obs]; - std::size_t fb_count = num_rounds * row_weight; - std::size_t base = num_rounds * obs_fb_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]; - } - 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_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++; - } - } - - // With inlined feedback, the boundary detector for record - // (fixed_offset + x) is extended with its declared last-round records. - std::size_t fb_begin = 0; - std::size_t fb_end = 0; - if (fb_offsets.size() > 0) { - fb_begin = fb_offsets[fixed_offset + x]; - fb_end = fb_offsets[fixed_offset + x + 1]; - } - - std::vector support(support_weight + 1 + - (fb_end - fb_begin)); - 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]; - } - } - for (std::size_t i = fb_begin; i < fb_end; ++i) - support[support_idx++] = final_syndrome[fb_indices[i]]; - - 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_flat, numCols)); } } // namespace cudaq::qec diff --git a/libs/qec/lib/device/memory_circuit.h b/libs/qec/lib/device/memory_circuit.h index e21ac943b..a15b25842 100644 --- a/libs/qec/lib/device/memory_circuit.h +++ b/libs/qec/lib/device/memory_circuit.h @@ -29,20 +29,16 @@ 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 -/// @param fb_indices CSR herald-column indices for the detector inlined -/// feedback layout. `fb_indices[fb_offsets[j] .. fb_offsets[j+1])` -/// lists the herald record columns XORed into record j's cross-round -/// detector (records of the earlier round) and into its final boundary -/// detector (records of the last round). -/// @param fb_offsets CSR row offsets for the detector inlined feedback layout. -/// Size is 0 (no feedback; legacy `cudaq::detectors` structure) or -/// numCols + 1 with numCols = numAncx + numAncz. -/// @param obs_fb_indices CSR record-column indices for the observable inlined -/// feedback layout. `obs_fb_indices[obs_fb_offsets[m] .. -/// obs_fb_offsets[m+1])` lists the record columns XORed into logical -/// observable m on every round. -/// @param obs_fb_offsets CSR row offsets for the observable inlined feedback -/// layout. Size is 0 (no feedback) or num_observables + 1. +/// @param feedback_flat Row-major flattened inlined feedback matrix for +/// detectors. Size is 0 (no feedback; legacy detector structure) or +/// numCols*numCols with numCols = numAncx + numAncz. Entry (j, k) = 1 +/// means the cross-round detector for record j additionally XORs +/// record k of the earlier round, and the final boundary detector for +/// record j additionally XORs record k of the last round. +/// @param obs_feedback_flat Row-major flattened inlined feedback matrix for +/// logical observables. Size is 0 (no feedback) or +/// num_observables*numCols. Entry (m, k) = 1 means logical observable +/// m additionally XORs record k of every round. __qpu__ void memory_circuit(const code::stabilizer_round &stabilizer_round, const code::one_qubit_encoding &statePrep, std::size_t numData, std::size_t numAncx, @@ -52,8 +48,6 @@ __qpu__ void memory_circuit(const code::stabilizer_round &stabilizer_round, const std::vector &obs_matrix_flat, std::size_t num_observables, bool measure_in_x_basis, - const std::vector &fb_indices, - const std::vector &fb_offsets, - const std::vector &obs_fb_indices, - const std::vector &obs_fb_offsets); + const std::vector &feedback_flat, + const std::vector &obs_feedback_flat); } // namespace cudaq::qec diff --git a/libs/qec/lib/experiments.cpp b/libs/qec/lib/experiments.cpp index 39a465475..efb9add83 100644 --- a/libs/qec/lib/experiments.cpp +++ b/libs/qec/lib/experiments.cpp @@ -254,6 +254,29 @@ stabilizer_detector_ranges(std::size_t numRounds, std::size_t numZStabs, return ranges; } +/// @brief Validate and flatten (row-major) an inlined feedback tensor +/// declared by a code. An empty tensor (the identity default) flattens to an +/// empty vector; a non-empty tensor must have shape +/// [expected_rows, expected_cols]. +static std::vector +flatten_feedback_tensor(const cudaqx::tensor &t, + std::size_t expected_rows, std::size_t expected_cols, + const std::string &name) { + if (t.size() == 0) + return {}; + if (t.rank() != 2 || t.shape()[0] != expected_rows || + t.shape()[1] != expected_cols) { + std::string actual; + for (std::size_t d = 0; d < t.rank(); ++d) + actual += (d ? ", " : "") + std::to_string(t.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."); + } + return std::vector(t.data(), t.data() + t.size()); +} + /// @brief Given a memory circuit setup, sample syndrome and data qubit /// measurements. This is the main driver function that all of the /// sample_memory_circuit overloads invoke. @@ -309,29 +332,29 @@ sample_memory_circuit(const code &code, operation statePrep, const std::size_t numCols = numAncx + numAncz; - auto fb_layout = apply_inlined_feedback( - code.get_inlined_feedback(), code.get_observable_inlined_feedback(), - numCols, num_obs); + auto feedback_flat = flatten_feedback_tensor( + code.get_inlined_feedback(), numCols, numCols, "get_inlined_feedback()"); + auto obs_feedback_flat = + flatten_feedback_tensor(code.get_observable_inlined_feedback(), num_obs, + numCols, "get_observable_inlined_feedback()"); // 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, fb_layout.detector_indices, fb_layout.detector_offsets, - fb_layout.observable_indices, fb_layout.observable_offsets); + 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_flat, obs_feedback_flat); // 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, - fb_layout.detector_indices, fb_layout.detector_offsets, - fb_layout.observable_indices, fb_layout.observable_offsets); + auto result = + cudaq::sample(opts, memory_circuit, stabRound, prep, numData, numAncx, + numAncz, numRounds, xVec, zVec, obs_flat, num_obs, + !is_z_prep, feedback_flat, obs_feedback_flat); // mzTable[shot, meas_idx] = raw 0/1 outcome; shape (numShots, // numMeasPerShot). Measurement layout per shot: numRounds*numCols ancilla, @@ -488,17 +511,18 @@ dem_from_memory_circuit(const code &code, operation statePrep, logical_obs.data() + logical_obs.size()); const std::size_t numAnc = numAncx + numAncz; - auto fb_layout = apply_inlined_feedback( - code.get_inlined_feedback(), code.get_observable_inlined_feedback(), - numAnc, num_obs); + auto feedback_flat = flatten_feedback_tensor( + code.get_inlined_feedback(), numAnc, numAnc, "get_inlined_feedback()"); + auto obs_feedback_flat = + flatten_feedback_tensor(code.get_observable_inlined_feedback(), num_obs, + numAnc, "get_observable_inlined_feedback()"); 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, - fb_layout.detector_indices, fb_layout.detector_offsets, - fb_layout.observable_indices, fb_layout.observable_offsets); + feedback_flat, obs_feedback_flat); 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/unittests/test_qec.cpp b/libs/qec/unittests/test_qec.cpp index 429f960b3..7bba70357 100644 --- a/libs/qec/unittests/test_qec.cpp +++ b/libs/qec/unittests/test_qec.cpp @@ -11,9 +11,11 @@ #include #include #include +#include #include #include "cuda-qx/core/library_utils.h" +#include "cudaq/algorithms/dem.h" #include "cudaq/qec/codes/surface_code.h" #include "cudaq/qec/experiments.h" #include "cudaq/qec/pcm_utils.h" @@ -2362,7 +2364,7 @@ TEST(QECCodeTester, checkInlinedFeedbackToyDem) { TEST(InlinedFeedbackLayout, EmptyTensorsGiveEmptyLayout) { cudaqx::tensor empty; - auto layout = cudaq::qec::apply_inlined_feedback(empty, empty, 2, 1); + auto layout = cudaq::qec::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()); @@ -2376,7 +2378,7 @@ TEST(InlinedFeedbackLayout, ToyMatricesCsr) { // 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::apply_inlined_feedback(fb, obs_fb, 2, 1); + auto layout = cudaq::qec::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})); @@ -2390,7 +2392,7 @@ TEST(InlinedFeedbackLayout, MultiRowWeights) { fb.at({0, 2}) = 1; fb.at({2, 1}) = 1; cudaqx::tensor empty; - auto layout = cudaq::qec::apply_inlined_feedback(fb, empty, 3, 0); + auto layout = cudaq::qec::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()); @@ -2399,7 +2401,7 @@ TEST(InlinedFeedbackLayout, MultiRowWeights) { TEST(InlinedFeedbackLayout, AllZeroMatrixGivesZeroWeightRows) { cudaqx::tensor fb({2, 2}); cudaqx::tensor empty; - auto layout = cudaq::qec::apply_inlined_feedback(fb, empty, 2, 0); + auto layout = cudaq::qec::build_inlined_feedback_layout(fb, empty, 2, 0); EXPECT_TRUE(layout.detector_indices.empty()); EXPECT_EQ(layout.detector_offsets, (std::vector{0, 0, 0})); } @@ -2407,9 +2409,126 @@ TEST(InlinedFeedbackLayout, AllZeroMatrixGivesZeroWeightRows) { TEST(InlinedFeedbackLayout, ThrowsOnWrongShape) { cudaqx::tensor bad({2, 3}); // expected [2 x 2] cudaqx::tensor empty; - EXPECT_THROW(cudaq::qec::apply_inlined_feedback(bad, empty, 2, 0), + EXPECT_THROW(cudaq::qec::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::apply_inlined_feedback(empty, bad_obs, 2, 1), + EXPECT_THROW(cudaq::qec::build_inlined_feedback_layout(empty, bad_obs, 2, 1), std::runtime_error); } + +// 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_flat, + const std::vector &obs_feedback_flat); +} // 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. Build the measurements-to-detectors matrix (M2D) from the kernel +// and compare it, row by row, against the M2D reconstructed from the host +// layout plus the known emission order. If either the device fold or the host +// layout drifts, this test fails. +TEST(InlinedFeedbackLayout, HostLayoutMatchesKernelM2D) { + feedback_toy_code toy(/*declare_feedback=*/true); + const std::size_t numRounds = 4; + 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. prep0 is a + // Z-basis memory, so measure_in_x_basis = false and the fixed records are the + // Z ancillas (fixed_offset = 0). --- + auto &prep = toy.get_operation( + cudaq::qec::operation::prep0); + 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 = 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 = toy.get_observable_inlined_feedback(); + std::vector feedback_flat(feedback.data(), + feedback.data() + feedback.size()); + std::vector obs_feedback_flat( + obs_feedback.data(), obs_feedback.data() + obs_feedback.size()); + + 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, + /*measure_in_x_basis=*/false, feedback_flat, obs_feedback_flat); + + // --- 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 layout = cudaq::qec::build_inlined_feedback_layout( + feedback, obs_feedback, numCols, num_obs); + 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; + // Round-0 fixed detector: the fixed record (fixed_offset + 0) of round 0. + expected.push_back({rec(0, 0)}); + // 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 (row 0): the last-round fixed + // record, the data readouts in the stabilizer support, and the layout-row-0 + // heralds taken from the last round. + { + std::set row{rec(numRounds - 1, 0)}; + for (std::size_t q = 0; q < numData; ++q) + if (zVec[q] != 0) + row.insert(numRounds * numCols + q); + for (auto k : + herald_cols(layout.detector_indices, layout.detector_offsets, 0)) + 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; + } +} From 43b6bd18464b7c8aaaa69c0e216fca52b9a464c7 Mon Sep 17 00:00:00 2001 From: kvmto Date: Sat, 11 Jul 2026 03:07:04 +0200 Subject: [PATCH 4/4] split between two types of stabiizer Signed-off-by: kvmto --- libs/qec/include/cudaq/qec/code.h | 29 +++++++-- libs/qec/lib/code.cpp | 6 +- libs/qec/lib/experiments.cpp | 18 ++++-- libs/qec/python/bindings/py_code.cpp | 37 ++++++++--- libs/qec/python/tests/test_code.py | 67 +++++++++++++++++++- libs/qec/unittests/test_qec.cpp | 94 ++++++++++++++++++++++++++-- 6 files changed, 221 insertions(+), 30 deletions(-) diff --git a/libs/qec/include/cudaq/qec/code.h b/libs/qec/include/cudaq/qec/code.h index 3a9eefce7..e41177a76 100644 --- a/libs/qec/include/cudaq/qec/code.h +++ b/libs/qec/include/cudaq/qec/code.h @@ -214,16 +214,33 @@ class code : public cudaqx::extension_point { /// formed from same-record comparisons only. virtual cudaqx::tensor get_inlined_feedback() const; - /// @brief Get the inlined feedback matrix for logical observables + /// @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). Entry (m, k) = 1 means logical observable m - /// additionally XORs record k of every round. - /// @return Tensor representing the observable inlined feedback matrix. - /// An empty tensor (the default) means no feedback is applied to the - /// logical observables. - virtual cudaqx::tensor get_observable_inlined_feedback() const; + /// 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). 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 diff --git a/libs/qec/lib/code.cpp b/libs/qec/lib/code.cpp index bb6b720a6..39be5efed 100644 --- a/libs/qec/lib/code.cpp +++ b/libs/qec/lib/code.cpp @@ -65,7 +65,11 @@ cudaqx::tensor code::get_inlined_feedback() const { return cudaqx::tensor(); } -cudaqx::tensor code::get_observable_inlined_feedback() const { +cudaqx::tensor code::get_observable_inlined_feedback_z() const { + return cudaqx::tensor(); +} + +cudaqx::tensor code::get_observable_inlined_feedback_x() const { return cudaqx::tensor(); } diff --git a/libs/qec/lib/experiments.cpp b/libs/qec/lib/experiments.cpp index efb9add83..d4a9f1780 100644 --- a/libs/qec/lib/experiments.cpp +++ b/libs/qec/lib/experiments.cpp @@ -334,9 +334,12 @@ sample_memory_circuit(const code &code, operation statePrep, auto feedback_flat = flatten_feedback_tensor( code.get_inlined_feedback(), numCols, numCols, "get_inlined_feedback()"); - auto obs_feedback_flat = - flatten_feedback_tensor(code.get_observable_inlined_feedback(), num_obs, - numCols, "get_observable_inlined_feedback()"); + auto obs_feedback_flat = flatten_feedback_tensor( + is_z_prep ? code.get_observable_inlined_feedback_z() + : code.get_observable_inlined_feedback_x(), + num_obs, numCols, + 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 @@ -513,9 +516,12 @@ dem_from_memory_circuit(const code &code, operation statePrep, const std::size_t numAnc = numAncx + numAncz; auto feedback_flat = flatten_feedback_tensor( code.get_inlined_feedback(), numAnc, numAnc, "get_inlined_feedback()"); - auto obs_feedback_flat = - flatten_feedback_tensor(code.get_observable_inlined_feedback(), num_obs, - numAnc, "get_observable_inlined_feedback()"); + auto obs_feedback_flat = flatten_feedback_tensor( + is_z_prep ? code.get_observable_inlined_feedback_z() + : code.get_observable_inlined_feedback_x(), + num_obs, numAnc, + is_z_prep ? "get_observable_inlined_feedback_z()" + : "get_observable_inlined_feedback_x()"); cudaq::dem_options dem_opts; dem_opts.decompose_errors = decompose_errors; diff --git a/libs/qec/python/bindings/py_code.cpp b/libs/qec/python/bindings/py_code.cpp index 3c63816b6..52a59b43b 100644 --- a/libs/qec/python/bindings/py_code.cpp +++ b/libs/qec/python/bindings/py_code.cpp @@ -51,7 +51,7 @@ static cudaqx::tensor feedbackArrayToTensor(nb::object array, class PyCode : public qec::code { public: - NB_TRAMPOLINE(qec::code, 8); + NB_TRAMPOLINE(qec::code, 9); protected: // Trampoline methods for the optional inlined-feedback virtuals: forward to @@ -66,13 +66,22 @@ class PyCode : public qec::code { return NBBase::get_inlined_feedback(); } - cudaqx::tensor get_observable_inlined_feedback() const override { + cudaqx::tensor get_observable_inlined_feedback_z() const override { nb::detail::ticket nb_ticket(nb_trampoline, - "get_observable_inlined_feedback", false); + "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"); - return NBBase::get_observable_inlined_feedback(); + "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 @@ -237,12 +246,20 @@ class PyCodeHandle : public qec::code { "get_inlined_feedback"); } - cudaqx::tensor get_observable_inlined_feedback() const override { - if (!nb::hasattr(pyCode, "get_observable_inlined_feedback")) - return code::get_observable_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")(), - "get_observable_inlined_feedback"); + pyCode.attr("get_observable_inlined_feedback_x")(), + "get_observable_inlined_feedback_x"); } protected: diff --git a/libs/qec/python/tests/test_code.py b/libs/qec/python/tests/test_code.py index 0ba39a425..0a08c3bc4 100644 --- a/libs/qec/python/tests/test_code.py +++ b/libs/qec/python/tests/test_code.py @@ -392,7 +392,7 @@ def get_num_z_stabilizers(self): def get_inlined_feedback(self): return np.array([[0, 1], [0, 0]], dtype=np.uint8) - def get_observable_inlined_feedback(self): + def get_observable_inlined_feedback_z(self): # Exercise the bridge's dtype coercion: int64 instead of uint8. return np.array([[0, 1]]) @@ -473,5 +473,70 @@ def get_num_z_stabilizers(self): numRounds=4) +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). + @qec.code('py-feedback-toy-basis-select') + class FeedbackToyBasisSelect: + + 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 + + def get_inlined_feedback(self): + return np.array([[0, 1], [0, 0]], dtype=np.uint8) + + def get_observable_inlined_feedback_z(self): + return np.array([[0, 1]], dtype=np.uint8) + + def get_observable_inlined_feedback_x(self): + # 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) + + 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) + assert not np.any(data[:, 0] ^ data[:, 1]) + finally: + cudaq.reset_target() + + if __name__ == "__main__": pytest.main() diff --git a/libs/qec/unittests/test_qec.cpp b/libs/qec/unittests/test_qec.cpp index 7bba70357..4b44344d8 100644 --- a/libs/qec/unittests/test_qec.cpp +++ b/libs/qec/unittests/test_qec.cpp @@ -2200,9 +2200,12 @@ TEST(QECCodeTester, checkInlinedFeedbackDefaults) { auto feedback = steane->get_inlined_feedback(); EXPECT_EQ(feedback.rank(), 0); EXPECT_EQ(feedback.size(), 0); - auto obs_feedback = steane->get_observable_inlined_feedback(); - EXPECT_EQ(obs_feedback.rank(), 0); - EXPECT_EQ(obs_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 @@ -2287,15 +2290,38 @@ class feedback_toy_code : public cudaq::qec::code { return feedback; } - cudaqx::tensor get_observable_inlined_feedback() const override { + cudaqx::tensor get_observable_inlined_feedback_z() const override { if (!declare_feedback) - return code::get_observable_inlined_feedback(); + return code::get_observable_inlined_feedback_z(); cudaqx::tensor feedback({1, 2}); feedback.at({0, 1}) = 1; 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) { @@ -2324,6 +2350,62 @@ TEST(QECCodeTester, checkInlinedFeedbackToyMemoryCircuit) { EXPECT_EQ(data.at({shot, 0}) ^ data.at({shot, 1}), 0) << "shot " << shot; } +TEST(QECCodeTester, checkInlinedFeedbackObservableBasisSelection) { + // The observable-feedback getter must be selected by the memory basis: a + // Z-basis (prep0) memory queries only get_observable_inlined_feedback_z() + // and never get_observable_inlined_feedback_x(). The toy device kernels + // register only prep0, so selection is proved from the Z path in both + // directions. + auto valid_obs = [] { + cudaqx::tensor t({1, 2}); + t.at({0, 1}) = 1; + return t; + }; + // Non-empty and the wrong shape for [num_obs=1 x numCols=2], so any code that + // flattens it trips flatten_feedback_tensor's shape check. + auto wrong_shape_obs = [] { return cudaqx::tensor({3, 5}); }; + + // (a) Valid Z matrix, deliberately wrong-shape X matrix: the Z-memory run + // must succeed exactly as the positive toy test, proving the X getter is + // never queried on the Z path. + { + feedback_toy_basis_select_code toy(valid_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); + for (std::size_t shot = 0; shot < numShots; ++shot) + EXPECT_EQ(data.at({shot, 0}) ^ data.at({shot, 1}), 0) << "shot " << shot; + } + + // (converse) Wrong-shape Z matrix, valid X matrix: the same Z-memory run + // must now throw the shape-validation error, proving the Z getter IS the one + // queried on the Z path (the valid X matrix is irrelevant here). + { + feedback_toy_basis_select_code toy(wrong_shape_obs(), valid_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(); + } + } +} + TEST(QECCodeTester, checkInlinedFeedbackToyNegativeControl) { // The identical toy without the feedback declaration: the uncorrected // byproduct makes the record-0 detectors non-deterministic and stim must @@ -2464,7 +2546,7 @@ TEST(InlinedFeedbackLayout, HostLayoutMatchesKernelM2D) { logical_obs.data() + logical_obs.size()); auto feedback = toy.get_inlined_feedback(); - auto obs_feedback = toy.get_observable_inlined_feedback(); + auto obs_feedback = toy.get_observable_inlined_feedback_z(); std::vector feedback_flat(feedback.data(), feedback.data() + feedback.size()); std::vector obs_feedback_flat(