From 31c4eda57f2b22bcf4be7493f5bcf323d9755102 Mon Sep 17 00:00:00 2001 From: Ben Howe Date: Wed, 15 Jul 2026 03:03:33 +0000 Subject: [PATCH] Replace msm execution contexts with cudaq::dem_from_kernel in app examples The surface_code-1/2/3 C++ app examples and the surface_code_1.py Python example built their detector error models by running the kernels under the "msm_size"/"msm" execution contexts and hand-assembling the detector and observable matrices from the raw measurement syndrome matrix. The new cudaq::dem_from_kernel API supersedes that: the demo kernels now annotate themselves with cudaq::detector/detectors and cudaq::logical_observable when a declare_detectors flag is set (detectors pair each round against the previous one, with the lock-in round as the reference; prep0's deterministic first-round Z stabilizers are declared as singles where the old code used them), and the hosts parse the returned Stim DEM text with dem_from_stim_text. The runtime D_sparse now comes straight from the analysis' measurements-to-detectors (m2d) map instead of generate_timelike_sparse_detector_matrix, remapped into each decoder's own enqueue stream for surface_code-3's split Z/X decoders (which run one dem_from_kernel pass per decoder via declare_detectors_z/_x flags). Notable constraints baked into the kernels: - dem_from_kernel rejects kernels tagged qubitMeasurementFeedback, so the final data measurements are packed branch-free instead of via to_integer(to_bools(...)), whose discriminate-into-call flow sets the tag. - The Python kernels avoid measure-handle lists crossing kernel-call boundaries and avoid iterating returned handle lists: both leave loops and heap traffic that do not survive the full loop unrolling the adaptive QIR profile (quantinuum target) requires. The DEM slice of demo_circuit_qpu is instead a single custom_memory_circuit_stabs call covering lock-in plus one decoder window, which is gate-for-gate identical to the live path. The now-unused compute_msm/construct_mz_table Python bindings are removed. Testing: all 54 app_examples ctest cases pass (local, sliding-window, cqr in-process and two-process, quantinuum-emulate at d=3/5, prep0/prepp), and the 12 surface_code-1-test.py pytest cases pass in ~10s (with /usr/local/cudaq and the built python dir on PYTHONPATH). Co-Authored-By: Claude Fable 5 Signed-off-by: Ben Howe --- libs/qec/python/bindings/py_decoder.cpp | 45 -- libs/qec/python/cudaq_qec/__init__.py | 2 - .../realtime/app_examples/surface_code-1.cpp | 254 +++++------ .../realtime/app_examples/surface_code-2.cpp | 283 ++++++------ .../realtime/app_examples/surface_code-3.cpp | 405 ++++++++---------- .../realtime/app_examples/surface_code_1.py | 362 +++++++++------- 6 files changed, 617 insertions(+), 734 deletions(-) diff --git a/libs/qec/python/bindings/py_decoder.cpp b/libs/qec/python/bindings/py_decoder.cpp index b646645b8..12a03d1c1 100644 --- a/libs/qec/python/bindings/py_decoder.cpp +++ b/libs/qec/python/bindings/py_decoder.cpp @@ -5,7 +5,6 @@ * This source code and the accompanying materials are made available under * * the terms of the Apache License 2.0 which accompanies this distribution. * ******************************************************************************/ -#include "common/ExecutionContext.h" #include "cuda-qx/core/kwargs_utils.h" #include "cuda-qx/core/library_utils.h" #include "type_casters.h" @@ -1404,50 +1403,6 @@ void bindDecoder(nb::module_ &mod) { [Internal] Load local simulation realtime decoder library. )pbdoc"); - qecmod.def( - "compute_msm", - [](std::function kernel, bool verbose = false) { - cudaq::ExecutionContext ctx_msm_size("msm_size"); - auto &platform = cudaq::get_platform(); - platform.with_execution_context(ctx_msm_size, kernel); - if (!ctx_msm_size.msm_dimensions.has_value()) { - throw std::runtime_error("No MSM dimensions found"); - } - if (ctx_msm_size.msm_dimensions.value().second == 0) { - throw std::runtime_error("No MSM dimensions found"); - } - cudaq::ExecutionContext ctx_msm("msm"); - ctx_msm.msm_dimensions = ctx_msm_size.msm_dimensions; - platform.with_execution_context(ctx_msm, kernel); - - auto msm_as_strings = ctx_msm.result.sequential_data(); - if (verbose) { - printf("MSM Dimensions: %ld measurements x %ld error mechanisms\n", - ctx_msm.msm_dimensions.value().first, - ctx_msm.msm_dimensions.value().second); - for (std::size_t i = 0; i < ctx_msm.msm_dimensions.value().first; - i++) { - for (std::size_t j = 0; j < ctx_msm.msm_dimensions.value().second; - j++) { - printf("%c", msm_as_strings[j][i] == '1' ? '1' : '.'); - } - printf("\n"); - } - } - return std::make_tuple(msm_as_strings, ctx_msm.msm_dimensions.value(), - ctx_msm.msm_probabilities.value(), - ctx_msm.msm_prob_err_id.value()); - }, - ""); - qecmod.def( - "construct_mz_table", - [](const std::vector &msm_as_strings) { - cudaqx::tensor mzTable(msm_as_strings); - mzTable = mzTable.transpose(); - return cudaq::python::copyCUDAQXTensorToPyArray(mzTable); - }, - ""); - qecmod.def( "generate_timelike_sparse_detector_matrix", [](std::uint32_t num_syndromes_per_round, std::uint32_t num_rounds, diff --git a/libs/qec/python/cudaq_qec/__init__.py b/libs/qec/python/cudaq_qec/__init__.py index 310785e7d..144886ad0 100644 --- a/libs/qec/python/cudaq_qec/__init__.py +++ b/libs/qec/python/cudaq_qec/__init__.py @@ -107,8 +107,6 @@ def checked_decode_batch(self, *args, **kwargs): simplify_pcm = qecrt.simplify_pcm sort_pcm_columns = qecrt.sort_pcm_columns pcm_extend_to_n_rounds = qecrt.pcm_extend_to_n_rounds -compute_msm = qecrt.compute_msm -construct_mz_table = qecrt.construct_mz_table generate_timelike_sparse_detector_matrix = qecrt.generate_timelike_sparse_detector_matrix pcm_to_sparse_vec = qecrt.pcm_to_sparse_vec diff --git a/libs/qec/unittests/realtime/app_examples/surface_code-1.cpp b/libs/qec/unittests/realtime/app_examples/surface_code-1.cpp index a5ab846d2..5a688da7e 100644 --- a/libs/qec/unittests/realtime/app_examples/surface_code-1.cpp +++ b/libs/qec/unittests/realtime/app_examples/surface_code-1.cpp @@ -17,9 +17,9 @@ #include "cudaq/qec/realtime/decoding.h" #include "cudaq/qec/realtime/decoding_config.h" #include -#include #include #include +#include #include #include #include @@ -106,13 +106,12 @@ static int g_syndromes_per_shot = 0; // #define MANUALLY_INJECT_ERRORS void save_dem_to_file(const cudaq::qec::detector_error_model &dem, + const std::vector &d_sparse, std::string dem_filename, uint64_t numSyndromesPerRound, uint64_t numLogical, const std::string &decoder_type, int sw_window_size, int sw_step_size, bool use_relay_bp) { cudaq::qec::decoding::config::multi_decoder_config multi_config; for (uint64_t i = 0; i < numLogical; i++) { - // We actually send 1 additional round in this example, so add 1. - auto numRounds = dem.num_detectors() / numSyndromesPerRound + 1; cudaq::qec::decoding::config::decoder_config config; config.id = i; config.type = decoder_type; // Use parameter instead of hardcoded @@ -121,8 +120,7 @@ void save_dem_to_file(const cudaq::qec::detector_error_model &dem, config.H_sparse = cudaq::qec::pcm_to_sparse_vec(dem.detector_error_matrix); config.O_sparse = cudaq::qec::pcm_to_sparse_vec(dem.observables_flips_matrix); - config.D_sparse = cudaq::qec::generate_timelike_sparse_detector_matrix( - numSyndromesPerRound, numRounds, /*include_first_round=*/false); + config.D_sparse = d_sparse; if (decoder_type == "nv-qldpc-decoder") { config.decoder_custom_args = @@ -323,19 +321,38 @@ se_x_ft(cudaq::qec::patch logicalQubit, return results; } -__qpu__ void custom_memory_circuit_stabs( +// Runs the stabilizer measurement rounds and returns the last round's +// combined syndrome so the caller can chain it into a later call. When +// declare_detectors is set (DEM generation), every round declares one +// cross-round detector per syndrome bit against the previous round, in the +// same [Z..., X...] order the syndromes are enqueued. prev_syndrome supplies +// the reference round for the first transition (pass an empty vector for +// none, e.g. for the lock-in round itself). +__qpu__ std::vector custom_memory_circuit_stabs( cudaq::qview<> data, cudaq::qview<> xstab_anc, cudaq::qview<> zstab_anc, std::size_t numRounds, const std::vector &cnot_schedX_flat, const std::vector &cnot_schedZ_flat, bool enqueue_syndromes, bool do_errors_after_non_last_rounds, double p_spam, int logical_qubit_idx, - int decoder_window) { + int decoder_window, bool declare_detectors, + const std::vector &prev_syndrome) { // Create the logical patch patch logical(data, xstab_anc, zstab_anc); - std::vector combined_syndrome(xstab_anc.size() + - zstab_anc.size()); + + // Local copy of the reference syndrome (kernel vector parameters are + // read-only, and the round loop below reassigns this as it advances). The + // local is always full-size so the reassignment never changes its length; + // have_prev says whether it currently holds a valid reference round. + bool have_prev = prev_syndrome.size() == xstab_anc.size() + zstab_anc.size(); + std::vector prev(xstab_anc.size() + zstab_anc.size()); + if (have_prev) { + for (std::size_t k = 0; k < prev.size(); ++k) + prev[k] = prev_syndrome[k]; + } // Handle the stabilizer lock-in round (numRounds == 1) if (numRounds == 1) { + std::vector combined_syndrome(xstab_anc.size() + + zstab_anc.size()); auto syndrome_z = se_z_ft(logical, cnot_schedZ_flat); auto syndrome_x = se_x_ft(logical, cnot_schedX_flat); int i = 0; @@ -347,7 +364,10 @@ __qpu__ void custom_memory_circuit_stabs( cudaq::qec::decoding::enqueue_syndromes( /*decoder_id=*/logical_qubit_idx, combined_syndrome); } - return; + if (declare_detectors && have_prev) { + cudaq::detectors(prev, combined_syndrome); + } + return combined_syndrome; } // Process rounds window by window for the main measurement rounds @@ -358,12 +378,14 @@ __qpu__ void custom_memory_circuit_stabs( // For window_idx > 0, enqueue the last syndrome from previous window first if (window_idx > 0 && enqueue_syndromes) { cudaq::qec::decoding::enqueue_syndromes( - /*decoder_id=*/logical_qubit_idx, combined_syndrome); + /*decoder_id=*/logical_qubit_idx, prev); } // Process the current window rounds for (std::size_t round = window_idx * decoder_window; round < (window_idx + 1) * decoder_window; round++) { + std::vector combined_syndrome(xstab_anc.size() + + zstab_anc.size()); auto syndrome_z = se_z_ft(logical, cnot_schedZ_flat); auto syndrome_x = se_x_ft(logical, cnot_schedX_flat); int i = 0; @@ -375,6 +397,11 @@ __qpu__ void custom_memory_circuit_stabs( cudaq::qec::decoding::enqueue_syndromes( /*decoder_id=*/logical_qubit_idx, combined_syndrome); } + if (declare_detectors && have_prev) { + cudaq::detectors(prev, combined_syndrome); + } + prev = combined_syndrome; + have_prev = true; #if PER_SHOT_DEBUG debug_print_syndromes(syndrome_x_int, syndrome_z_int); #endif @@ -393,16 +420,24 @@ __qpu__ void custom_memory_circuit_stabs( } } } + return prev; } +// When declare_detectors is set (only meaningful with numLogical = 1 and +// allow_device_calls = false), the kernel annotates itself for DEM +// generation via cudaq::dem_from_kernel: every stabilizer round declares +// cross-round detectors against the previous round (the lock-in round is the +// first reference), and the Z logical observable is declared over the final +// data measurements at z_obs_indices. __qpu__ std::int64_t -demo_circuit_qpu(bool allow_device_calls, +demo_circuit_qpu(bool allow_device_calls, bool declare_detectors, const cudaq::qec::code::one_qubit_encoding &statePrep, std::size_t numData, std::size_t numAncx, std::size_t numAncz, std::size_t numRounds, std::size_t numLogical, const std::vector &cnot_schedX_flat, const std::vector &cnot_schedZ_flat, - double p_spam, bool apply_corrections, int decoder_window) { + double p_spam, bool apply_corrections, int decoder_window, + const std::vector &z_obs_indices) { #if PER_SHOT_DEBUG debug_start_shot(); #endif @@ -428,18 +463,25 @@ demo_circuit_qpu(bool allow_device_calls, statePrep(logical); } - // Do 1 stabilizer round to lock in the stabilizers + // Do 1 stabilizer round to lock in the stabilizers. Its syndrome is the + // reference round for the first cross-round detectors when + // declare_detectors is set (DEM generation always uses numLogical = 1). + std::vector lockin_syndrome(numAncx + numAncz); { for (int i = 0; i < numLogical; i++) { auto subData = data.slice(i * numData, numData); auto subXstab_anc = xstab_anc.slice(i * numAncx, numAncx); auto subZstab_anc = zstab_anc.slice(i * numAncz, numAncz); - custom_memory_circuit_stabs( + std::vector no_prev(0); + auto syndrome = custom_memory_circuit_stabs( subData, subXstab_anc, subZstab_anc, /*numRounds=*/1, cnot_schedX_flat, cnot_schedZ_flat, /*enqueue_syndromes=*/allow_device_calls, - /*do_errors_after_non_last_rounds=*/false, p_spam, i, decoder_window); + /*do_errors_after_non_last_rounds=*/false, p_spam, i, decoder_window, + /*declare_detectors=*/false, no_prev); + if (i == 0) + lockin_syndrome = syndrome; } } @@ -459,10 +501,12 @@ demo_circuit_qpu(bool allow_device_calls, auto subXstab_anc = xstab_anc.slice(i * numAncx, numAncx); auto subZstab_anc = zstab_anc.slice(i * numAncz, numAncz); - custom_memory_circuit_stabs( - subData, subXstab_anc, subZstab_anc, numRounds, cnot_schedX_flat, - cnot_schedZ_flat, /*enqueue_syndromes=*/allow_device_calls, - /*do_errors_after_non_last_rounds=*/true, p_spam, i, decoder_window); + custom_memory_circuit_stabs(subData, subXstab_anc, subZstab_anc, numRounds, + cnot_schedX_flat, cnot_schedZ_flat, + /*enqueue_syndromes=*/allow_device_calls, + /*do_errors_after_non_last_rounds=*/true, + p_spam, i, decoder_window, declare_detectors, + lockin_syndrome); } // Only apply corrections after processing all windows @@ -492,7 +536,20 @@ demo_circuit_qpu(bool allow_device_calls, ret <<= numData; auto subData = data.slice(i * numData, numData); auto subMeas = mz(subData); - ret |= cudaq::to_integer(cudaq::to_bools(subMeas)); + if (declare_detectors && i == 0) { + std::vector zlog(z_obs_indices.size()); + for (std::size_t k = 0; k < z_obs_indices.size(); ++k) + zlog[k] = subMeas[z_obs_indices[k]]; + cudaq::logical_observable(zlog, /*observable_index=*/0); + } + // Pack the measured bits branch-free (bit j = data qubit j, the same + // LSB-first order as cudaq::to_integer). Routing the measurement results + // through a call (to_bools/to_integer) or a branch would tag this kernel + // with qubitMeasurementFeedback, which cudaq::dem_from_kernel rejects. + for (std::size_t j = 0; j < subMeas.size(); j++) { + std::uint64_t bitval = subMeas[j]; + ret |= bitval << j; + } } // The remaining bits are allocated to the number of corrections. ret |= num_corrections << (numData * numLogical); @@ -543,131 +600,57 @@ void demo_circuit_host(const cudaq::qec::code &code, int distance, cudaq::noise_model noise; - // First get the MSM + // The Z logical observable's data-qubit support (row 0 of the Z observables + // matrix); demo_circuit_qpu declares the matching logical_observable over + // the final data measurements when generating the DEM. + auto obs_matrix = code.get_observables_z(); + std::vector z_obs_indices; + for (std::size_t col = 0; col < obs_matrix.shape()[1]; ++col) + if (obs_matrix.at({0, col})) + z_obs_indices.push_back(col); + + // First generate (or load) the DEM cudaq::qec::detector_error_model dem; if (load_dem) { load_dem_from_file(dem_filename, dem, numLogical); } else { if (p_spam == 0.0) { - printf("p_spam is 0.0, cannot get the MSM\n"); + printf("p_spam is 0.0, cannot generate the DEM\n"); exit(0); } - cudaq::ExecutionContext ctx_msm_size("msm_size"); - ctx_msm_size.noiseModel = &noise; - auto &platform = cudaq::get_platform(); - platform.with_execution_context(ctx_msm_size, [&] { - // Always use numLogical = 1 for the MSM - cudaq::qec::qpu::demo_circuit_qpu( - /*allow_device_calls=*/false, prep, numData, numAncx, numAncz, - decoder_window, // Use decoder_window instead of numRounds for DEM - // generation - /*numLogical=*/1, cnot_schedX_flat, cnot_schedZ_flat, p_spam, - /*apply_corrections=*/false, decoder_window); - }); - if (!ctx_msm_size.msm_dimensions.has_value()) { - throw std::runtime_error("No MSM dimensions found"); - } - if (ctx_msm_size.msm_dimensions.value().second == 0) { - throw std::runtime_error("No MSM dimensions found"); - } - cudaq::ExecutionContext ctx_msm("msm"); - ctx_msm.noiseModel = &noise; - ctx_msm.msm_dimensions = ctx_msm_size.msm_dimensions; - platform.with_execution_context(ctx_msm, [&] { - // Always use numLogical = 1 for the MSM - cudaq::qec::qpu::demo_circuit_qpu( - /*allow_device_calls=*/false, prep, numData, numAncx, numAncz, - decoder_window, // Use decoder_window instead of numRounds for DEM - // generation - /*numLogical=*/1, cnot_schedX_flat, cnot_schedZ_flat, p_spam, - /*apply_corrections=*/false, decoder_window); - }); - - auto msm_as_strings = ctx_msm.result.sequential_data(); - printf("MSM Dimensions: %ld measurements x %ld error mechanisms\n", - ctx_msm.msm_dimensions.value().first, - ctx_msm.msm_dimensions.value().second); - for (std::size_t i = 0; i < ctx_msm.msm_dimensions.value().first; i++) { - for (std::size_t j = 0; j < ctx_msm.msm_dimensions.value().second; j++) { - printf("%c", msm_as_strings[j][i] == '1' ? '1' : '.'); - } - printf("\n"); - } - // Populate error rates and error IDs - dem.error_rates = std::move(ctx_msm.msm_probabilities.value()); - dem.error_ids = std::move(ctx_msm.msm_prob_err_id.value()); - - cudaqx::tensor mzTable(msm_as_strings); - mzTable = mzTable.transpose(); - printf("mzTable:\n"); - mzTable.dump_bits(); - // Subtract the number of data qubits to get the number of syndrome - // measurements. - std::size_t totalNumSyndromes = mzTable.shape()[0] - distance * distance; - std::size_t numNoiseMechs = mzTable.shape()[1]; - std::size_t numSyndromesPerRound = distance * distance - 1; - if (totalNumSyndromes % numSyndromesPerRound != 0) { - throw std::runtime_error("Num syndromes per round is not a divisor of " - "the number of syndrome measurements"); - } - std::size_t numRoundsOfSyndromData = - totalNumSyndromes / numSyndromesPerRound; - if (numRoundsOfSyndromData != - decoder_window + 1) { // Use decoder_window instead of numRounds - throw std::runtime_error("Num rounds of syndrome data [" + - std::to_string(numRoundsOfSyndromData) + - "] is not equal to the decoder_window + 1[" + - std::to_string(decoder_window + 1) + "]"); - } - dem.detector_error_matrix = cudaqx::tensor( - {decoder_window * numSyndromesPerRound, - numNoiseMechs}); // Use decoder_window instead of numRounds - // There should be (decoder_window + 1) rounds of data in MSM. - // TODO: [feature] Good candidate. Auto-generating the detector error - // matrix. Currently, we need to manually construct the detector error - // matrix by copying the measurements from the MSM. - for (std::size_t round = 0; round < decoder_window; - round++) { // Use decoder_window instead of numRounds - for (std::size_t syndrome = 0; syndrome < numSyndromesPerRound; - syndrome++) { - for (std::size_t noise_mech = 0; noise_mech < numNoiseMechs; - noise_mech++) { - dem.detector_error_matrix.at( - {round * numSyndromesPerRound + syndrome, noise_mech}) = - mzTable.at( - {(round + 0) * numSyndromesPerRound + syndrome, noise_mech}) ^ - mzTable.at( - {(round + 1) * numSyndromesPerRound + syndrome, noise_mech}); - } - } - } - auto first_data_row = - (decoder_window + 1) * - numSyndromesPerRound; // Use decoder_window instead of numRounds - cudaqx::tensor msm_obs( - {mzTable.shape()[0] - first_data_row, numNoiseMechs}); - for (std::size_t row = first_data_row; row < mzTable.shape()[0]; row++) - for (std::size_t col = 0; col < numNoiseMechs; col++) - msm_obs.at({row - first_data_row, col}) = mzTable.at({row, col}); - - // Populate dem.observables_flips_matrix by converting the physical data - // qubit measurements to logical observables. - auto obs_matrix = code.get_observables_z(); - printf("obs_matrix:\n"); - obs_matrix.dump_bits(); - dem.observables_flips_matrix = obs_matrix.dot(msm_obs) % 2; - printf("numSyndromesPerRound: %ld\n", numSyndromesPerRound); + cudaq::M2DSparseMatrix m2d; + cudaq::M2OSparseMatrix m2o; + std::string dem_text = cudaq::dem_from_kernel( + cudaq::qec::qpu::demo_circuit_qpu, &noise, m2d, m2o, + /*allow_device_calls=*/false, + /*declare_detectors=*/true, prep, numData, numAncx, numAncz, + decoder_window, // Use decoder_window instead of numRounds for DEM + // generation + /*numLogical=*/1, cnot_schedX_flat, cnot_schedZ_flat, p_spam, + /*apply_corrections=*/false, decoder_window, z_obs_indices); + dem = cudaq::qec::dem_from_stim_text(dem_text); + auto numSyndromesPerRound = distance * distance - 1; dem.canonicalize_for_rounds(numSyndromesPerRound, /*remove_zero_syndrome_errors=*/true); + // The runtime detector matrix comes straight from the analysis' + // measurements-to-detectors map: row d lists the (chronological, and thus + // enqueue-ordered) measurement indices whose XOR forms detector d. + std::vector d_sparse; + for (const auto &row : m2d.rows) { + for (auto m : row) + d_sparse.push_back(static_cast(m)); + d_sparse.push_back(-1); + } + printf("dem.detector_error_matrix:\n"); dem.detector_error_matrix.dump_bits(); printf("dem.observables_flips_matrix:\n"); dem.observables_flips_matrix.dump_bits(); if (save_dem) { - save_dem_to_file(dem, dem_filename, numSyndromesPerRound, numLogical, - decoder_type, sw_window_size, sw_step_size, + save_dem_to_file(dem, d_sparse, dem_filename, numSyndromesPerRound, + numLogical, decoder_type, sw_window_size, sw_step_size, use_relay_bp); return; } @@ -907,19 +890,20 @@ void demo_circuit_host(const cudaq::qec::code &code, int distance, run_result = cudaq::get_platform().is_remote() ? cudaq::run(numShots, cudaq::qec::qpu::demo_circuit_qpu, - /*allow_device_calls=*/true, prep, numData, numAncx, + /*allow_device_calls=*/true, + /*declare_detectors=*/false, prep, numData, numAncx, numAncz, numRounds, numLogical, cnot_schedX_flat, cnot_schedZ_flat, p_spam, /*apply_corrections=*/true, - decoder_window) + decoder_window, z_obs_indices) : cudaq::run(numShots, noise, cudaq::qec::qpu::demo_circuit_qpu, - /*allow_device_calls=*/true, prep, numData, numAncx, + /*allow_device_calls=*/true, + /*declare_detectors=*/false, prep, numData, numAncx, numAncz, numRounds, numLogical, cnot_schedX_flat, cnot_schedZ_flat, p_spam, /*apply_corrections=*/true, - decoder_window); + decoder_window, z_obs_indices); } printf("Result size: %ld\n", run_result.size()); std::vector> logical_results; - auto obs_matrix = code.get_observables_z(); int num_non_zero_values = 0; std::int64_t num_corrections = 0; for (int i = 0; i < run_result.size(); i++) { diff --git a/libs/qec/unittests/realtime/app_examples/surface_code-2.cpp b/libs/qec/unittests/realtime/app_examples/surface_code-2.cpp index b61e174f2..57632f93f 100644 --- a/libs/qec/unittests/realtime/app_examples/surface_code-2.cpp +++ b/libs/qec/unittests/realtime/app_examples/surface_code-2.cpp @@ -17,8 +17,8 @@ #include "cudaq/qec/realtime/decoding.h" #include "cudaq/qec/realtime/decoding_config.h" #include -#include #include +#include #include // Whether or not to put calls to debug functions in the QIR program. You cannot @@ -183,19 +183,39 @@ se_x_ft(cudaq::qec::patch logicalQubit, return results; } -__qpu__ void custom_memory_circuit_stabs( +// Runs the stabilizer measurement rounds and returns the last round's +// combined syndrome so the caller can chain it into a later call. When +// declare_detectors is set (DEM generation), every round declares detectors: +// cross-round pairs against the previous round in the same [Z..., X...] order +// the syndromes are enqueued or, when there is no reference round yet (the +// lock-in round after prep0), one single-measurement detector per +// Z-stabilizer, which is deterministic in the prepared basis. prev_syndrome +// supplies the reference round (pass an empty vector for none). +__qpu__ std::vector custom_memory_circuit_stabs( cudaq::qview<> data, cudaq::qview<> xstab_anc, cudaq::qview<> zstab_anc, std::size_t numRounds, const std::vector &cnot_schedX_flat, const std::vector &cnot_schedZ_flat, bool enqueue_syndromes, - bool do_errors_after_non_last_rounds, double p_spam, - int logical_qubit_idx) { + bool do_errors_after_non_last_rounds, double p_spam, int logical_qubit_idx, + bool declare_detectors, + const std::vector &prev_syndrome) { // Create the logical patch patch logical(data, xstab_anc, zstab_anc); - std::vector combined_syndrome(xstab_anc.size() + - zstab_anc.size()); + + // Local copy of the reference syndrome (kernel vector parameters are + // read-only, and the round loop below reassigns this as it advances). The + // local is always full-size so the reassignment never changes its length; + // have_prev says whether it currently holds a valid reference round. + bool have_prev = prev_syndrome.size() == xstab_anc.size() + zstab_anc.size(); + std::vector prev(xstab_anc.size() + zstab_anc.size()); + if (have_prev) { + for (std::size_t k = 0; k < prev.size(); ++k) + prev[k] = prev_syndrome[k]; + } // Generate syndrome data for (std::size_t round = 0; round < numRounds; round++) { + std::vector combined_syndrome(xstab_anc.size() + + zstab_anc.size()); auto syndrome_z = se_z_ft(logical, cnot_schedZ_flat); auto syndrome_x = se_x_ft(logical, cnot_schedX_flat); int i = 0; @@ -207,6 +227,19 @@ __qpu__ void custom_memory_circuit_stabs( cudaq::qec::decoding::enqueue_syndromes( /*decoder_id=*/logical_qubit_idx, combined_syndrome); } + if (declare_detectors) { + if (have_prev) { + cudaq::detectors(prev, combined_syndrome); + } else { + // In this Z-basis (prep0) example the Z stabilizers measured by the + // very first round are deterministic, so they are detectors on their + // own; the X stabilizers only project a random frame. + for (std::size_t k = 0; k < zstab_anc.size(); ++k) + cudaq::detector(combined_syndrome[k]); + } + } + prev = combined_syndrome; + have_prev = true; #if PER_SHOT_DEBUG debug_print_syndromes(syndrome_x_int, syndrome_z_int); #endif @@ -223,16 +256,22 @@ __qpu__ void custom_memory_circuit_stabs( #endif } } + return prev; } -__qpu__ std::int64_t -demo_circuit_qpu(bool allow_device_calls, - const cudaq::qec::code::one_qubit_encoding &statePrep, - std::size_t numData, std::size_t numAncx, std::size_t numAncz, - std::size_t numRounds, std::size_t numLogical, - const std::vector &cnot_schedX_flat, - const std::vector &cnot_schedZ_flat, - double p_spam, bool apply_corrections) { +// When declare_detectors is set (only meaningful with numLogical = 1 and +// allow_device_calls = false), the kernel annotates itself for DEM +// generation via cudaq::dem_from_kernel: the lock-in round declares the +// deterministic Z-stabilizer singles, every later round declares cross-round +// detectors against the previous round, and the Z logical observable is +// declared over the final data measurements at z_obs_indices. +__qpu__ std::int64_t demo_circuit_qpu( + bool allow_device_calls, bool declare_detectors, + const cudaq::qec::code::one_qubit_encoding &statePrep, std::size_t numData, + std::size_t numAncx, std::size_t numAncz, std::size_t numRounds, + std::size_t numLogical, const std::vector &cnot_schedX_flat, + const std::vector &cnot_schedZ_flat, double p_spam, + bool apply_corrections, const std::vector &z_obs_indices) { #if PER_SHOT_DEBUG debug_start_shot(); #endif @@ -258,18 +297,25 @@ demo_circuit_qpu(bool allow_device_calls, statePrep(logical); } - // Do 1 stabilizer round to lock in the stabilizers + // Do 1 stabilizer round to lock in the stabilizers. Its syndrome is the + // reference round for the first cross-round detectors when + // declare_detectors is set (DEM generation always uses numLogical = 1). + std::vector lockin_syndrome(numAncx + numAncz); { for (int i = 0; i < numLogical; i++) { auto subData = data.slice(i * numData, numData); auto subXstab_anc = xstab_anc.slice(i * numAncx, numAncx); auto subZstab_anc = zstab_anc.slice(i * numAncz, numAncz); - custom_memory_circuit_stabs( + std::vector no_prev(0); + auto syndrome = custom_memory_circuit_stabs( subData, subXstab_anc, subZstab_anc, /*numRounds=*/1, cnot_schedX_flat, cnot_schedZ_flat, /*enqueue_syndromes=*/allow_device_calls, - /*do_errors_after_non_last_rounds=*/false, p_spam, i); + /*do_errors_after_non_last_rounds=*/false, p_spam, i, + declare_detectors, no_prev); + if (i == 0) + lockin_syndrome = syndrome; } } @@ -289,10 +335,11 @@ demo_circuit_qpu(bool allow_device_calls, auto subXstab_anc = xstab_anc.slice(i * numAncx, numAncx); auto subZstab_anc = zstab_anc.slice(i * numAncz, numAncz); - custom_memory_circuit_stabs( - subData, subXstab_anc, subZstab_anc, numRounds, cnot_schedX_flat, - cnot_schedZ_flat, /*enqueue_syndromes=*/allow_device_calls, - /*do_errors_after_non_last_rounds=*/true, p_spam, i); + custom_memory_circuit_stabs(subData, subXstab_anc, subZstab_anc, numRounds, + cnot_schedX_flat, cnot_schedZ_flat, + /*enqueue_syndromes=*/allow_device_calls, + /*do_errors_after_non_last_rounds=*/true, + p_spam, i, declare_detectors, lockin_syndrome); } if (allow_device_calls && apply_corrections) { @@ -321,7 +368,20 @@ demo_circuit_qpu(bool allow_device_calls, ret <<= numData; auto subData = data.slice(i * numData, numData); auto subMeas = mz(subData); - ret |= cudaq::to_integer(cudaq::to_bools(subMeas)); + if (declare_detectors && i == 0) { + std::vector zlog(z_obs_indices.size()); + for (std::size_t k = 0; k < z_obs_indices.size(); ++k) + zlog[k] = subMeas[z_obs_indices[k]]; + cudaq::logical_observable(zlog, /*observable_index=*/0); + } + // Pack the measured bits branch-free (bit j = data qubit j, the same + // LSB-first order as cudaq::to_integer). Routing the measurement results + // through a call (to_bools/to_integer) or a branch would tag this kernel + // with qubitMeasurementFeedback, which cudaq::dem_from_kernel rejects. + for (std::size_t j = 0; j < subMeas.size(); j++) { + std::uint64_t bitval = subMeas[j]; + ret |= bitval << j; + } } // The remaining bits are allocated to the number of corrections. ret |= num_corrections << (numData * numLogical); @@ -367,132 +427,34 @@ void demo_circuit_host(const cudaq::qec::code &code, int distance, cudaq::noise_model noise; - // First get the MSM + // The Z logical observable's data-qubit support (row 0 of the Z observables + // matrix); demo_circuit_qpu declares the matching logical_observable over + // the final data measurements when generating the DEM. + auto obs_matrix = code.get_observables_z(); + std::vector z_obs_indices; + for (std::size_t col = 0; col < obs_matrix.shape()[1]; ++col) + if (obs_matrix.at({0, col})) + z_obs_indices.push_back(col); + + // First generate (or load) the DEM cudaq::qec::detector_error_model dem; if (load_dem) { load_dem_from_file(dem_filename, dem, numLogical); } else { if (p_spam == 0.0) { - printf("p_spam is 0.0, cannot get the MSM\n"); + printf("p_spam is 0.0, cannot generate the DEM\n"); exit(0); } - cudaq::ExecutionContext ctx_msm_size("msm_size"); - ctx_msm_size.noiseModel = &noise; - auto &platform = cudaq::get_platform(); - platform.with_execution_context(ctx_msm_size, [&] { - // Always use numLogical = 1 for the MSM - cudaq::qec::qpu::demo_circuit_qpu( - /*allow_device_calls=*/false, prep, numData, numAncx, numAncz, - numRounds, - /*numLogical=*/1, cnot_schedX_flat, cnot_schedZ_flat, p_spam, - /*apply_corrections=*/false); - }); - if (!ctx_msm_size.msm_dimensions.has_value()) { - throw std::runtime_error("No MSM dimensions found"); - } - if (ctx_msm_size.msm_dimensions.value().second == 0) { - throw std::runtime_error("No MSM dimensions found"); - } - cudaq::ExecutionContext ctx_msm("msm"); - ctx_msm.noiseModel = &noise; - ctx_msm.msm_dimensions = ctx_msm_size.msm_dimensions; - platform.with_execution_context(ctx_msm, [&] { - // Always use numLogical = 1 for the MSM - cudaq::qec::qpu::demo_circuit_qpu( - /*allow_device_calls=*/false, prep, numData, numAncx, numAncz, - numRounds, - /*numLogical=*/1, cnot_schedX_flat, cnot_schedZ_flat, p_spam, - /*apply_corrections=*/false); - }); - - auto msm_as_strings = ctx_msm.result.sequential_data(); - printf("MSM Dimensions: %ld measurements x %ld error mechanisms\n", - ctx_msm.msm_dimensions.value().first, - ctx_msm.msm_dimensions.value().second); - for (std::size_t i = 0; i < ctx_msm.msm_dimensions.value().first; i++) { - for (std::size_t j = 0; j < ctx_msm.msm_dimensions.value().second; j++) { - printf("%c", msm_as_strings[j][i] == '1' ? '1' : '.'); - } - printf("\n"); - } - // Populate error rates and error IDs - dem.error_rates = std::move(ctx_msm.msm_probabilities.value()); - dem.error_ids = std::move(ctx_msm.msm_prob_err_id.value()); - - cudaqx::tensor mzTable(msm_as_strings); - mzTable = mzTable.transpose(); - printf("mzTable:\n"); - mzTable.dump_bits(); - // Subtract the number of data qubits to get the number of syndrome - // measurements. - std::size_t totalNumSyndromes = mzTable.shape()[0] - distance * distance; - std::size_t numNoiseMechs = mzTable.shape()[1]; + cudaq::M2DSparseMatrix m2d; + cudaq::M2OSparseMatrix m2o; + std::string dem_text = cudaq::dem_from_kernel( + cudaq::qec::qpu::demo_circuit_qpu, &noise, m2d, m2o, + /*allow_device_calls=*/false, + /*declare_detectors=*/true, prep, numData, numAncx, numAncz, numRounds, + /*numLogical=*/1, cnot_schedX_flat, cnot_schedZ_flat, p_spam, + /*apply_corrections=*/false, z_obs_indices); + dem = cudaq::qec::dem_from_stim_text(dem_text); std::size_t numSyndromesPerRound = distance * distance - 1; - - std::size_t numRoundsOfSyndromData = - totalNumSyndromes / numSyndromesPerRound; - if (numRoundsOfSyndromData != numRounds + 1) { - throw std::runtime_error("Num rounds of syndrome data [" + - std::to_string(numRoundsOfSyndromData) + - "] is not equal to the number of rounds + 1[" + - std::to_string(numRounds + 1) + "]"); - } - - // There should be (numRounds + 1) rounds of data in MSM. - // This corresponds to numRounds + measurements during state prep - // Not every measurement during stateprep is a detector, but some - // may be. - // In this Z-basis surface code case, the Z stabs during state prep - // are detectors. - // Skip the X stabs during the first round. - std::size_t numDetectors = - numSyndromesPerRound * numRounds + numSyndromesPerRound / 2; - - dem.detector_error_matrix = - cudaqx::tensor({numDetectors, numNoiseMechs}); - // Grab first half of first "round" - std::size_t r0_offset = 0; - for (std::size_t syndrome = 0; syndrome < numSyndromesPerRound / 2; - syndrome++) { - for (std::size_t noise_mech = 0; noise_mech < numNoiseMechs; - noise_mech++) { - // round 0 - dem.detector_error_matrix.at({r0_offset, noise_mech}) = - mzTable.at({syndrome, noise_mech}); - } - r0_offset += 1; - } - - // Grab all of rounds >=1. - for (std::size_t round = 0; round < numRounds; round++) { - for (std::size_t syndrome = 0; syndrome < numSyndromesPerRound; - syndrome++) { - for (std::size_t noise_mech = 0; noise_mech < numNoiseMechs; - noise_mech++) { - dem.detector_error_matrix.at( - {round * numSyndromesPerRound + syndrome + r0_offset, - noise_mech}) = - mzTable.at( - {(round + 0) * numSyndromesPerRound + syndrome, noise_mech}) ^ - mzTable.at( - {(round + 1) * numSyndromesPerRound + syndrome, noise_mech}); - } - } - } - auto first_data_row = (numRounds + 1) * numSyndromesPerRound; - cudaqx::tensor msm_obs( - {mzTable.shape()[0] - first_data_row, numNoiseMechs}); - for (std::size_t row = first_data_row; row < mzTable.shape()[0]; row++) - for (std::size_t col = 0; col < numNoiseMechs; col++) - msm_obs.at({row - first_data_row, col}) = mzTable.at({row, col}); - - // Populate dem.observables_flips_matrix by converting the physical data - // qubit measurements to logical observables. - auto obs_matrix = code.get_observables_z(); - printf("obs_matrix:\n"); - obs_matrix.dump_bits(); - dem.observables_flips_matrix = obs_matrix.dot(msm_obs) % 2; - printf("numSyndromesPerRound: %ld\n", numSyndromesPerRound); dem.canonicalize_for_rounds(numSyndromesPerRound, /*remove_zero_syndrome_errors=*/true); @@ -501,25 +463,17 @@ void demo_circuit_host(const cudaq::qec::code &code, int distance, printf("dem.observables_flips_matrix:\n"); dem.observables_flips_matrix.dump_bits(); - // Prep0 means that first round of Z stabs should be deterministic - // These are measured first. - std::vector first_round; - for (int i = 0; i < numSyndromesPerRound / 2; i++) { - first_round.push_back(i); - first_round.push_back(-1); - } - - // TO DO: - // Does numRounds include first round? - std::vector det_mat = - cudaq::qec::generate_timelike_sparse_detector_matrix( - numSyndromesPerRound, numRounds + 1, first_round); - - printf("detector_matrix with first round:\n"); - for (int i = 0; i < det_mat.size(); i++) { - printf("%ld ", det_mat[i]); + // The runtime detector matrix comes straight from the analysis' + // measurements-to-detectors map: row d lists the (chronological, and thus + // enqueue-ordered) measurement indices whose XOR forms detector d. This + // reproduces the round-0 Z-stabilizer singles followed by the timelike + // pairs of the declared detectors. + std::vector det_mat; + for (const auto &row : m2d.rows) { + for (auto m : row) + det_mat.push_back(static_cast(m)); + det_mat.push_back(-1); } - printf("\n"); if (save_dem) { save_dem_to_file(dem, det_mat, dem_filename, numSyndromesPerRound, @@ -544,16 +498,19 @@ void demo_circuit_host(const cudaq::qec::code &code, int distance, auto run_result = cudaq::get_platform().is_remote() ? cudaq::run(numShots, cudaq::qec::qpu::demo_circuit_qpu, - /*allow_device_calls=*/true, prep, numData, numAncx, + /*allow_device_calls=*/true, + /*declare_detectors=*/false, prep, numData, numAncx, numAncz, numRounds, numLogical, cnot_schedX_flat, - cnot_schedZ_flat, p_spam, /*apply_corrections=*/true) + cnot_schedZ_flat, p_spam, /*apply_corrections=*/true, + z_obs_indices) : cudaq::run(numShots, noise, cudaq::qec::qpu::demo_circuit_qpu, - /*allow_device_calls=*/true, prep, numData, numAncx, + /*allow_device_calls=*/true, + /*declare_detectors=*/false, prep, numData, numAncx, numAncz, numRounds, numLogical, cnot_schedX_flat, - cnot_schedZ_flat, p_spam, /*apply_corrections=*/true); + cnot_schedZ_flat, p_spam, /*apply_corrections=*/true, + z_obs_indices); printf("Result size: %ld\n", run_result.size()); std::vector> logical_results; - auto obs_matrix = code.get_observables_z(); int num_non_zero_values = 0; std::int64_t num_corrections = 0; for (int i = 0; i < run_result.size(); i++) { diff --git a/libs/qec/unittests/realtime/app_examples/surface_code-3.cpp b/libs/qec/unittests/realtime/app_examples/surface_code-3.cpp index 655049e7f..c72a00dda 100644 --- a/libs/qec/unittests/realtime/app_examples/surface_code-3.cpp +++ b/libs/qec/unittests/realtime/app_examples/surface_code-3.cpp @@ -17,8 +17,8 @@ #include "cudaq/qec/realtime/decoding.h" #include "cudaq/qec/realtime/decoding_config.h" #include -#include #include +#include #include // Whether or not to put calls to debug functions in the QIR program. You cannot @@ -230,14 +230,35 @@ se_x_ft(cudaq::qec::patch logicalQubit, return results; } -__qpu__ void custom_memory_circuit_stabs( +// Runs the stabilizer measurement rounds and returns the last round's +// combined [Z..., X...] syndrome so the caller can chain it into a later +// call. When declare_detectors_z (or _x) is set (DEM generation of that +// stabilizer type's decoder model), every round declares that type's +// detectors: cross-round pairs against the previous round or, when there is +// no reference round yet (the lock-in round), one single-measurement +// detector per stabilizer if the type matches the prepared basis - the only +// type that is deterministic right after state prep. prev_syndrome supplies +// the reference round (pass an empty vector for none). +__qpu__ std::vector custom_memory_circuit_stabs( cudaq::qview<> data, cudaq::qview<> xstab_anc, cudaq::qview<> zstab_anc, std::size_t numRounds, const std::vector &cnot_schedX_flat, const std::vector &cnot_schedZ_flat, bool enqueue_syndromes, bool do_errors_after_non_last_rounds, double p_spam, int logical_qubit_idx, - bool is_on_Z_basis) { + bool is_on_Z_basis, bool declare_detectors_z, bool declare_detectors_x, + const std::vector &prev_syndrome) { patch logical(data, xstab_anc, zstab_anc); + // Local copy of the reference syndrome (kernel vector parameters are + // read-only, and the round loop below reassigns this as it advances). The + // local is always full-size so the reassignment never changes its length; + // have_prev says whether it currently holds a valid reference round. + bool have_prev = prev_syndrome.size() == xstab_anc.size() + zstab_anc.size(); + std::vector prev(xstab_anc.size() + zstab_anc.size()); + if (have_prev) { + for (std::size_t k = 0; k < prev.size(); ++k) + prev[k] = prev_syndrome[k]; + } + for (std::size_t round = 0; round < numRounds; round++) { auto syndrome_z = se_z_ft(logical, cnot_schedZ_flat); auto syndrome_x = se_x_ft(logical, cnot_schedX_flat); @@ -253,6 +274,35 @@ __qpu__ void custom_memory_circuit_stabs( /*decoder_id=*/2 * logical_qubit_idx + 1, syndrome_x); } + if (have_prev) { + if (declare_detectors_z) { + for (std::size_t k = 0; k < zstab_anc.size(); ++k) + cudaq::detector(prev[k], syndrome_z[k]); + } + if (declare_detectors_x) { + for (std::size_t k = 0; k < xstab_anc.size(); ++k) + cudaq::detector(prev[zstab_anc.size() + k], syndrome_x[k]); + } + } else { + // Lock-in round: only the prepared-basis stabilizers are deterministic. + if (declare_detectors_z && is_on_Z_basis) { + for (std::size_t k = 0; k < zstab_anc.size(); ++k) + cudaq::detector(syndrome_z[k]); + } + if (declare_detectors_x && !is_on_Z_basis) { + for (std::size_t k = 0; k < xstab_anc.size(); ++k) + cudaq::detector(syndrome_x[k]); + } + } + std::vector combined_syndrome(xstab_anc.size() + + zstab_anc.size()); + for (std::size_t k = 0; k < zstab_anc.size(); ++k) + combined_syndrome[k] = syndrome_z[k]; + for (std::size_t k = 0; k < xstab_anc.size(); ++k) + combined_syndrome[zstab_anc.size() + k] = syndrome_x[k]; + prev = combined_syndrome; + have_prev = true; + if (do_errors_after_non_last_rounds && round < numRounds - 1) { // spam_error(logical, p_spam, p_spam, p_spam); spam_error(logical, /*p_spam_data=*/p_spam, /*p_spam_ancx=*/0.001, @@ -270,17 +320,27 @@ __qpu__ void custom_memory_circuit_stabs( #endif } } + return prev; } -__qpu__ std::int64_t -demo_circuit_qpu(bool allow_device_calls, - const cudaq::qec::code::one_qubit_encoding &statePrep, - bool is_on_Z_basis, std::size_t numData, std::size_t numAncx, - std::size_t numAncz, std::size_t numRounds, - std::size_t numLogical, - const std::vector &cnot_schedX_flat, - const std::vector &cnot_schedZ_flat, - double p_spam, bool apply_corrections) { +// When declare_detectors_z (or _x) is set (only meaningful with +// numLogical = 1 and allow_device_calls = false), the kernel annotates +// itself for DEM generation of that stabilizer type's decoder model via +// cudaq::dem_from_kernel - one dem_from_kernel pass per decoder. The lock-in +// round declares the deterministic prepared-basis stabilizer singles, every +// later round declares cross-round detectors against the previous round, +// and, on the pass whose type matches the prepared basis, the logical +// observable is declared over the final data measurements at obs_indices +// (the data qubits are already rotated into that basis before the final +// measurement; the other basis' observable is not deterministic). +__qpu__ std::int64_t demo_circuit_qpu( + bool allow_device_calls, bool declare_detectors_z, bool declare_detectors_x, + const cudaq::qec::code::one_qubit_encoding &statePrep, bool is_on_Z_basis, + std::size_t numData, std::size_t numAncx, std::size_t numAncz, + std::size_t numRounds, std::size_t numLogical, + const std::vector &cnot_schedX_flat, + const std::vector &cnot_schedZ_flat, double p_spam, + bool apply_corrections, const std::vector &obs_indices) { #if PER_SHOT_DEBUG debug_start_shot(); #endif @@ -307,16 +367,24 @@ demo_circuit_qpu(bool allow_device_calls, statePrep(logical); } + // Do 1 stabilizer round to lock in the stabilizers. Its syndrome is the + // reference round for the first cross-round detectors when + // declare_detectors is set (DEM generation always uses numLogical = 1). + std::vector lockin_syndrome(numAncx + numAncz); { for (int i = 0; i < numLogical; i++) { auto subData = data.slice(i * numData, numData); auto subXstab_anc = xstab_anc.slice(i * numAncx, numAncx); auto subZstab_anc = zstab_anc.slice(i * numAncz, numAncz); - custom_memory_circuit_stabs( + std::vector no_prev(0); + auto syndrome = custom_memory_circuit_stabs( subData, subXstab_anc, subZstab_anc, /*numRounds=*/1, cnot_schedX_flat, cnot_schedZ_flat, /*enqueue_syndromes=*/allow_device_calls, - /*do_errors_after_non_last_rounds=*/false, p_spam, i, is_on_Z_basis); + /*do_errors_after_non_last_rounds=*/false, p_spam, i, is_on_Z_basis, + declare_detectors_z, declare_detectors_x, no_prev); + if (i == 0) + lockin_syndrome = syndrome; } } @@ -336,7 +404,8 @@ demo_circuit_qpu(bool allow_device_calls, custom_memory_circuit_stabs( subData, subXstab_anc, subZstab_anc, numRounds, cnot_schedX_flat, cnot_schedZ_flat, /*enqueue_syndromes=*/allow_device_calls, - /*do_errors_after_non_last_rounds=*/true, p_spam, i, is_on_Z_basis); + /*do_errors_after_non_last_rounds=*/true, p_spam, i, is_on_Z_basis, + declare_detectors_z, declare_detectors_x, lockin_syndrome); } std::uint16_t num_x_corrections = 0; @@ -387,12 +456,25 @@ demo_circuit_qpu(bool allow_device_calls, auto subMeas = mz(subData); + // The observable belongs to the pass declaring the prepared-basis + // detectors; the other basis' observable is not deterministic. + bool declare_observable = + is_on_Z_basis ? declare_detectors_z : declare_detectors_x; + if (declare_observable && i == 0) { + std::vector obs_log(obs_indices.size()); + for (std::size_t k = 0; k < obs_indices.size(); ++k) + obs_log[k] = subMeas[obs_indices[k]]; + cudaq::logical_observable(obs_log, /*observable_index=*/0); + } + int bit_offset = (numLogical - 1 - i) * numData; + // Pack the measured bits branch-free. Routing the measurement results + // through a branch (or a call like to_integer) would tag this kernel + // with qubitMeasurementFeedback, which cudaq::dem_from_kernel rejects. for (std::size_t j = 0; j < subMeas.size(); j++) { - if (subMeas[j]) { - ret |= (1ul << (bit_offset + j)); - } + std::uint64_t bitval = subMeas[j]; + ret |= bitval << (bit_offset + j); } } ret |= packed_counts << (numData * numLogical); @@ -449,6 +531,17 @@ void demo_circuit_host(const cudaq::qec::code &code, int distance, std::vector dem_z_vec, dem_x_vec; cudaq::qec::detector_error_model dem_z, dem_x; + // The prepared-basis logical observable's data-qubit support (row 0 of the + // matching observables matrix); demo_circuit_qpu declares the corresponding + // logical_observable over the final data measurements when generating the + // DEM. Only the prepared basis is deterministic, so only it is declared. + auto obs_matrix_on_basis = + is_on_Z_basis ? code.get_observables_z() : code.get_observables_x(); + std::vector obs_indices; + for (std::size_t col = 0; col < obs_matrix_on_basis.shape()[1]; ++col) + if (obs_matrix_on_basis.at({0, col})) + obs_indices.push_back(col); + if (load_dem) { // *** Pass numLogical to load_dem_from_file load_dem_from_file(dem_filename, dem_z_vec, dem_x_vec, numLogical); @@ -456,73 +549,9 @@ void demo_circuit_host(const cudaq::qec::code &code, int distance, dem_x = dem_x_vec.at(0); } else { if (p_spam == 0.0) { - printf("p_spam is 0.0, cannot get the MSM\n"); + printf("p_spam is 0.0, cannot generate the DEM\n"); exit(0); } - // ---------------DRY RUN TYPE-------------------------- - cudaq::ExecutionContext ctx_msm_size("msm_size"); - ctx_msm_size.noiseModel = &noise; - auto &platform = cudaq::get_platform(); - platform.with_execution_context(ctx_msm_size, [&] { - // Always use numLogical = 1 for the MSM - cudaq::qec::qpu::demo_circuit_qpu( - /*allow_device_calls=*/false, prep, is_on_Z_basis, numData, numAncx, - numAncz, numRounds, - /*numLogical=*/1, cnot_schedX_flat, cnot_schedZ_flat, p_spam, - /*apply_corrections=*/false); - }); - if (!ctx_msm_size.msm_dimensions.has_value()) { - throw std::runtime_error("No MSM dimensions found"); - } - if (ctx_msm_size.msm_dimensions.value().second == 0) { - throw std::runtime_error("No MSM dimensions found"); - } - - // -------------DATA COLLECT---------------------------- - - cudaq::ExecutionContext ctx_msm("msm"); - ctx_msm.noiseModel = &noise; - // Line that help us in preallocation - ctx_msm.msm_dimensions = ctx_msm_size.msm_dimensions; - platform.with_execution_context(ctx_msm, [&] { - // Always use numLogical = 1 for the MSM - cudaq::qec::qpu::demo_circuit_qpu( - /*allow_device_calls=*/false, prep, is_on_Z_basis, numData, numAncx, - numAncz, numRounds, - /*numLogical=*/1, cnot_schedX_flat, cnot_schedZ_flat, p_spam, - /*apply_corrections=*/false); - }); - - auto msm_as_strings = ctx_msm.result.sequential_data(); - printf("MSM Dimensions: %ld measurements x %ld error mechanisms\n", - ctx_msm.msm_dimensions.value().first, - ctx_msm.msm_dimensions.value().second); - - for (std::size_t i = 0; i < ctx_msm.msm_dimensions.value().first; i++) { - for (std::size_t j = 0; j < ctx_msm.msm_dimensions.value().second; j++) { - printf("%c", msm_as_strings[j][i] == '1' ? '1' : '.'); - } - printf("\n"); - } - - // ------------------------------------------------------------------------ - // + Split the DEM generation for X and Z errors - - // Populate error rates and error IDs - dem_z.error_rates = ctx_msm.msm_probabilities.value(); - dem_z.error_ids = ctx_msm.msm_prob_err_id.value(); - - dem_x.error_rates = ctx_msm.msm_probabilities.value(); - dem_x.error_ids = ctx_msm.msm_prob_err_id.value(); - - cudaqx::tensor mzTable(msm_as_strings); - mzTable = mzTable.transpose(); - printf("mzTable:\n"); - mzTable.dump_bits(); - // Subtract the number of data qubits to get the number of syndrome - // measurements. - std::size_t totalNumSyndromes = mzTable.shape()[0] - distance * distance; - std::size_t numNoiseMechs = mzTable.shape()[1]; std::size_t numSyndromesPerRound = distance * distance - 1; // in the case of the surface code, where each should be @@ -541,127 +570,37 @@ void demo_circuit_host(const cudaq::qec::code &code, int distance, } // ------------------------------------------------------------------------ - - std::size_t numRoundsOfSyndromData = - totalNumSyndromes / numSyndromesPerRound; - - if (numRoundsOfSyndromData != numRounds + 1) { - throw std::runtime_error("Num rounds of syndrome data [" + - std::to_string(numRoundsOfSyndromData) + - "] is not equal to the number of rounds + 1[" + - std::to_string(numRounds + 1) + "]"); - } - // ------------------------------------------------------------------------ - // NOTE: - // There should be (numRounds + 1) rounds of data in MSM. - // This corresponds to numRounds + measurements during state prep - // Not every measurement during stateprep is a detector, but some - // may be. - // In this Z-basis surface code case, the Z stabs during state prep - // are detectors. - // Skip the X stabs during the first round. - // ------------------------------------------------------------------------ - std::size_t numDetectors_z = - numSyndromesPerRound_z * - (is_on_Z_basis ? numRounds + 1 - : numRounds); // (numRounds +1 ) when on basis - std::size_t numDetectors_x = - numSyndromesPerRound_x * (is_on_Z_basis ? numRounds : numRounds + 1); - - // Build detector matrix for X errors (from Z stabs) - dem_z.detector_error_matrix = - tensor({numDetectors_z, numNoiseMechs}); - dem_x.detector_error_matrix = - tensor({numDetectors_x, numNoiseMechs}); - - // Build Z stabilizer detector matrix - std::size_t z_detector_idx = 0; - auto start_index_z = is_on_Z_basis ? 0 : 1; - for (std::size_t round = start_index_z; round < numRounds; round++) { - if (round == start_index_z) { - // Round 0: directly copy Z syndrome values without XORing - for (std::size_t z_stab = 0; z_stab < numSyndromesPerRound_z; - z_stab++) { - for (std::size_t mech = 0; mech < numNoiseMechs; mech++) { - dem_z.detector_error_matrix.at({z_detector_idx, mech}) = - mzTable.at({round * numSyndromesPerRound + z_stab, mech}); - } - z_detector_idx++; - } - } else { - // Rounds 1+: XOR current round with previous round for Z stabs - for (std::size_t z_stab = 0; z_stab < numSyndromesPerRound_z; - z_stab++) { - for (std::size_t mech = 0; mech < numNoiseMechs; mech++) { - dem_z.detector_error_matrix.at({z_detector_idx, mech}) = - mzTable.at({round * numSyndromesPerRound + z_stab, mech}) ^ - mzTable.at({(round - 1) * numSyndromesPerRound + z_stab, mech}); - } - z_detector_idx++; - } - } - } - - // Build X stabilizer detector matrix (skip round 0) - std::size_t x_detector_idx = 0; - std::size_t offset = numSyndromesPerRound_z; // X stabs start after Z stabs - auto start_index_x = is_on_Z_basis ? 1 : 0; - - for (std::size_t round = start_index_x; round < numRounds; - round++) { // start from round 1 - if (round == start_index_x) { - // Round 1: directly copy X syndrome values (first meaningful round) - for (std::size_t x_stab = 0; x_stab < numSyndromesPerRound_x; - x_stab++) { - for (std::size_t mech = 0; mech < numNoiseMechs; mech++) { - dem_x.detector_error_matrix.at({x_detector_idx, mech}) = mzTable.at( - {round * numSyndromesPerRound + offset + x_stab, mech}); - } - x_detector_idx++; - } - } else { - // Rounds 2+: XOR current round with previous round for X stabs - for (std::size_t x_stab = 0; x_stab < numSyndromesPerRound_x; - x_stab++) { - for (std::size_t mech = 0; mech < numNoiseMechs; mech++) { - dem_x.detector_error_matrix.at({x_detector_idx, mech}) = - mzTable.at( - {round * numSyndromesPerRound + offset + x_stab, mech}) ^ - mzTable.at( - {(round - 1) * numSyndromesPerRound + offset + x_stab, - mech}); - } - x_detector_idx++; - } - } - } - - // ------------------------------------------------------------------------ - // Data qubits Measurement Syndrome Matrix - auto first_data_row = (numRounds + 1) * numSyndromesPerRound; - // Here the (numRounds + 1) is to skip directly to the data measurements - cudaqx::tensor msm_obs( - {mzTable.shape()[0] - first_data_row, numNoiseMechs}); - for (std::size_t row = first_data_row; row < mzTable.shape()[0]; row++) - for (std::size_t col = 0; col < numNoiseMechs; col++) - msm_obs.at({row - first_data_row, col}) = mzTable.at({row, col}); - - // ------------------------------------------------------------------------ - - // Populate dem.observables_flips_matrix by converting the physical data - // qubit measurements to logical observables. - // We populate observables_flips_matrix for both Z and X DEMs - auto obs_matrix_z = code.get_observables_z(); - auto obs_matrix_x = code.get_observables_x(); - - printf("obs_matrix_z:\n"); - obs_matrix_z.dump_bits(); - printf("obs_matrix_x:\n"); - obs_matrix_x.dump_bits(); - - // Calculate observables flips for both DEMs - dem_z.observables_flips_matrix = obs_matrix_z.dot(msm_obs) % 2; - dem_x.observables_flips_matrix = obs_matrix_x.dot(msm_obs) % 2; + // One dem_from_kernel pass per decoder: the Z pass declares only the + // Z-stabilizer detectors (and, on a Z-basis prep, the observable), the X + // pass only the X-stabilizer ones. Each pass directly yields that + // decoder's model. + cudaq::M2DSparseMatrix m2d_z, m2d_x; + cudaq::M2OSparseMatrix m2o_z, m2o_x; + std::string dem_text_z = cudaq::dem_from_kernel( + cudaq::qec::qpu::demo_circuit_qpu, &noise, m2d_z, m2o_z, + /*allow_device_calls=*/false, + /*declare_detectors_z=*/true, /*declare_detectors_x=*/false, prep, + is_on_Z_basis, numData, numAncx, numAncz, numRounds, + /*numLogical=*/1, cnot_schedX_flat, cnot_schedZ_flat, p_spam, + /*apply_corrections=*/false, obs_indices); + dem_z = cudaq::qec::dem_from_stim_text(dem_text_z); + + std::string dem_text_x = cudaq::dem_from_kernel( + cudaq::qec::qpu::demo_circuit_qpu, &noise, m2d_x, m2o_x, + /*allow_device_calls=*/false, + /*declare_detectors_z=*/false, /*declare_detectors_x=*/true, prep, + is_on_Z_basis, numData, numAncx, numAncz, numRounds, + /*numLogical=*/1, cnot_schedX_flat, cnot_schedZ_flat, p_spam, + /*apply_corrections=*/false, obs_indices); + dem_x = cudaq::qec::dem_from_stim_text(dem_text_x); + + // Only the prepared-basis pass declares an observable (the other basis' + // observable is not deterministic), so the off-basis DEM comes back with + // zero observables. Give its decoder an explicit all-zero observables + // row; its corrections are never queried in this experiment. + auto &off_basis_dem = is_on_Z_basis ? dem_x : dem_z; + off_basis_dem.observables_flips_matrix = cudaqx::tensor( + {std::size_t{1}, off_basis_dem.detector_error_matrix.shape()[1]}); printf("numSyndromesPerRound_z: %ld\n", numSyndromesPerRound_z); printf("numSyndromesPerRound_x: %ld\n", numSyndromesPerRound_x); @@ -682,30 +621,30 @@ void demo_circuit_host(const cudaq::qec::code &code, int distance, printf("dem_x.observables_flips_matrix:\n"); dem_x.observables_flips_matrix.dump_bits(); - // In the X version we want the X mat to be 16X16 and the Z version 12X16. - - // Generate detector matrices, with a simplified API, in order to make the - // dynamic choice of which one should get the first one picked - std::vector det_mat_z = - cudaq::qec::generate_timelike_sparse_detector_matrix( - numSyndromesPerRound_z, numRounds + 1, is_on_Z_basis); - // numRounds + 1 we want the same num ber of columns - - std::vector det_mat_x = - cudaq::qec::generate_timelike_sparse_detector_matrix( - numSyndromesPerRound_x, numRounds + 1, !is_on_Z_basis); - - printf("detector_matrix_z :\n"); - for (int i = 0; i < det_mat_z.size(); i++) { - printf("%ld ", det_mat_z[i]); - } - printf("\n"); + // ------------------------------------------------------------------------ + // Per-decoder runtime detector matrices, straight from each pass' + // measurements-to-detectors map. m2d indexes the kernel's full + // chronological measurement stream ([Z..., X...] per round), while each + // decoder only receives its own stabilizer type, so remap every index + // into the decoder's own (Z-only or X-only) enqueue stream. + auto build_det_mat = [&](const cudaq::M2DSparseMatrix &m2d, bool z_type) { + std::vector out; + for (const auto &row : m2d.rows) { + for (auto g : row) { + const std::size_t round = g / numSyndromesPerRound; + const std::size_t pos = g % numSyndromesPerRound; + out.push_back(static_cast( + z_type ? round * numSyndromesPerRound_z + pos + : round * numSyndromesPerRound_x + + (pos - numSyndromesPerRound_z))); + } + out.push_back(-1); + } + return out; + }; + std::vector det_mat_z = build_det_mat(m2d_z, /*z_type=*/true); + std::vector det_mat_x = build_det_mat(m2d_x, /*z_type=*/false); - printf("detector_matrix_x :\n"); - for (int i = 0; i < det_mat_x.size(); i++) { - printf("%ld ", det_mat_x[i]); - } - printf("\n"); // ------------------------------------------------------------------------ if (save_dem) { @@ -739,15 +678,19 @@ void demo_circuit_host(const cudaq::qec::code &code, int distance, auto run_result = cudaq::get_platform().is_remote() ? cudaq::run(numShots, cudaq::qec::qpu::demo_circuit_qpu, - /*allow_device_calls=*/true, prep, is_on_Z_basis, + /*allow_device_calls=*/true, + /*declare_detectors_z=*/false, + /*declare_detectors_x=*/false, prep, is_on_Z_basis, numData, numAncx, numAncz, numRounds, numLogical, cnot_schedX_flat, cnot_schedZ_flat, p_spam, - /*apply_corrections=*/true) + /*apply_corrections=*/true, obs_indices) : cudaq::run(numShots, noise, cudaq::qec::qpu::demo_circuit_qpu, - /*allow_device_calls=*/true, prep, is_on_Z_basis, + /*allow_device_calls=*/true, + /*declare_detectors_z=*/false, + /*declare_detectors_x=*/false, prep, is_on_Z_basis, numData, numAncx, numAncz, numRounds, numLogical, cnot_schedX_flat, cnot_schedZ_flat, p_spam, - /*apply_corrections=*/true); + /*apply_corrections=*/true, obs_indices); // ------------------------------------------------------------------------ // Collecting data on experiment diff --git a/libs/qec/unittests/realtime/app_examples/surface_code_1.py b/libs/qec/unittests/realtime/app_examples/surface_code_1.py index 401699f50..70994a568 100644 --- a/libs/qec/unittests/realtime/app_examples/surface_code_1.py +++ b/libs/qec/unittests/realtime/app_examples/surface_code_1.py @@ -60,12 +60,11 @@ def sorted_stabilizer_ops_inplace_numpy(ops: List[cudaq.Operator]) -> None: ops[:] = [ops[i] for i in order] -def save_dem_to_file(dem, dem_filename, numSyndromesPerRound, num_logical): +def save_dem_to_file(dem, d_sparse, dem_filename, numSyndromesPerRound, + num_logical): multi_config = qec.multi_decoder_config() decoders = [] for i in range(num_logical): - # We actually send 1 additional round in this example, so add 1. - numRounds = dem.num_detectors() // numSyndromesPerRound + 1 config = qec.decoder_config() config.id = i config.type = "multi_error_lut" @@ -73,8 +72,7 @@ def save_dem_to_file(dem, dem_filename, numSyndromesPerRound, num_logical): config.syndrome_size = dem.num_detectors() config.H_sparse = qec.pcm_to_sparse_vec(dem.detector_error_matrix) config.O_sparse = qec.pcm_to_sparse_vec(dem.observables_flips_matrix) - config.D_sparse = qec.generate_timelike_sparse_detector_matrix( - numSyndromesPerRound, numRounds, False) + config.D_sparse = d_sparse multi_error_lut_config = qec.multi_error_lut_config() multi_error_lut_config.lut_error_depth = 2 config.set_decoder_custom_args(multi_error_lut_config) @@ -210,6 +208,14 @@ def se_x_ft(logical_qubit: patch, return results +# Runs the stabilizer measurement rounds. When declare_detectors is set (DEM +# generation), every round after the first declares one cross-round detector +# per syndrome bit against the previous round, in the same [Z..., X...] order +# the syndromes are enqueued; the first round is the reference round and +# declares none of its own. All detector bookkeeping (the measure-handle +# lists) is guarded by declare_detectors - a synthesized constant on hardware +# targets - so it folds away entirely for live runs: the adaptive QIR profile +# rejects the memory traffic that live measure-handle lists leave behind. @cudaq.kernel def custom_memory_circuit_stabs( data: cudaq.qview, @@ -224,27 +230,35 @@ def custom_memory_circuit_stabs( logical_qubit_idx: int, decoder_window: int, manually_inject_errors: bool, + declare_detectors: bool, ) -> None: # Create the logical patch logical = patch(data, xstab_anc, zstab_anc) - # Mirror the C++ idiom `std::vector(N)`: pre-allocate a - # list of unbound `measure_handle`s, each slot is overwritten with a real - # handle from `se_z_ft`/`se_x_ft` before any discrimination occurs. - # Can't use `measure_handle`: https://github.com/NVIDIA/cuda-quantum/issues/4527 - combined_syndrome = [False for i in range(len(xstab_anc) + len(zstab_anc))] - # Handle the stabilizer lock-in round (numRounds == 1) + # Bool mirror of the current round's syndrome for the enqueue calls, + # rebuilt in place each round with constant-bound loops (size queries on + # measure-handle lists, e.g. to_bools, do not survive the full loop + # unrolling that hardware targets require). + combined_bools = [False for i in range(len(xstab_anc) + len(zstab_anc))] + # The previous round's measurement handles (cross-round detectors only). + have_prev = False + prev = [ + cudaq.measure_handle() for _ in range(len(xstab_anc) + len(zstab_anc)) + ] + + # Handle the stabilizer lock-in round (numRounds == 1). The syndrome + # lists are read with constant-bound indexed loops (len of the qviews): + # iterating the returned lists directly gives loops bounded by the list + # size, which does not survive the full loop unrolling that hardware + # targets require. if num_rounds == 1: syndrome_z = se_z_ft(logical, cnot_schedZ_flat) syndrome_x = se_x_ft(logical, cnot_schedX_flat) - i = 0 - for s in syndrome_z: - combined_syndrome[i] = bool(s) - i += 1 - for s in syndrome_x: - combined_syndrome[i] = bool(s) - i += 1 + for k in range(len(zstab_anc)): + combined_bools[k] = bool(syndrome_z[k]) + for k in range(len(xstab_anc)): + combined_bools[len(zstab_anc) + k] = bool(syndrome_x[k]) if enqueue_synd: - qec.enqueue_syndromes_test(logical_qubit_idx, combined_syndrome, 0) + qec.enqueue_syndromes_test(logical_qubit_idx, combined_bools, 0) return # Process rounds window by window for the main measurement rounds @@ -253,24 +267,34 @@ def custom_memory_circuit_stabs( for window_idx in range(num_rounds // decoder_window): # For window_idx > 0, enqueue the last syndrome from previous window first if window_idx > 0 and enqueue_synd: - qec.enqueue_syndromes_test(logical_qubit_idx, combined_syndrome, 0) + qec.enqueue_syndromes_test(logical_qubit_idx, combined_bools, 0) # Process the current window rounds for round_idx in range(window_idx * decoder_window, (window_idx + 1) * decoder_window): syndrome_z = se_z_ft(logical, cnot_schedZ_flat) syndrome_x = se_x_ft(logical, cnot_schedX_flat) - i = 0 - for s in syndrome_z: - combined_syndrome[i] = bool(s) - i += 1 - for s in syndrome_x: - combined_syndrome[i] = bool(s) - i += 1 + for k in range(len(zstab_anc)): + combined_bools[k] = bool(syndrome_z[k]) + for k in range(len(xstab_anc)): + combined_bools[len(zstab_anc) + k] = bool(syndrome_x[k]) if enqueue_synd: - qec.enqueue_syndromes_test(logical_qubit_idx, combined_syndrome, - 0) + qec.enqueue_syndromes_test(logical_qubit_idx, combined_bools, 0) + if declare_detectors: + combined_syndrome = [ + cudaq.measure_handle() + for _ in range(len(xstab_anc) + len(zstab_anc)) + ] + for k in range(len(zstab_anc)): + combined_syndrome[k] = syndrome_z[k] + for k in range(len(xstab_anc)): + combined_syndrome[len(zstab_anc) + k] = syndrome_x[k] + if have_prev: + for k in range(len(xstab_anc) + len(zstab_anc)): + cudaq.detector(prev[k], combined_syndrome[k]) + prev = combined_syndrome + have_prev = True if do_errors_after_non_last_rounds and round_idx < ( window_idx + 1) * decoder_window - 1: @@ -281,9 +305,16 @@ def custom_memory_circuit_stabs( x(logical.data[3]) +# When declare_detectors is set (only meaningful with num_logical = 1 and +# allow_device_calls = False), the kernel annotates itself for DEM generation +# via cudaq.dem_from_kernel: every stabilizer round declares cross-round +# detectors against the previous round (the lock-in round is the first +# reference), and the Z logical observable is declared over the final data +# measurements at z_obs_indices. @cudaq.kernel def demo_circuit_qpu( allow_device_calls: bool, + declare_detectors: bool, #state_prep: Callable[[patch], None], num_data: int, num_ancx: int, @@ -296,6 +327,7 @@ def demo_circuit_qpu( apply_corrections: bool, decoder_window: int, manually_inject_errors: bool, + z_obs_indices: List[int], ) -> int: # if PER_SHOT_DEBUG: # debug_start_shot() @@ -321,56 +353,90 @@ def demo_circuit_qpu( logical = patch(sub_data, sub_x, sub_z) prep_0(logical) # FIXME: replace with state_prep(logical) - # One stabilizer round to lock in - for i in range(num_logical): - sub_data = data[i * num_data:(i + 1) * - num_data] # FIXME: all sub_data are incorrect - sub_x = xstab_anc[i * num_ancx:(i + 1) * num_ancx] # same other vectors - sub_z = zstab_anc[i * num_ancz:(i + 1) * num_ancz] + if declare_detectors: + # DEM-generation slice (num_logical == 1, no device calls): a single + # call covering the stabilizer lock-in round plus one decoder window + # (num_rounds + 1 rounds, with spam after every non-last round). This + # reproduces the exact gate and noise sequence of the live path in + # the else branch - the injected error round there equals the + # per-round spam - while keeping all measurement handles inside one + # kernel call, so the lock-in round can serve as the reference round + # for the first cross-round detectors. Keep the two paths in + # lockstep. + sub_data = data[0:num_data] + sub_x = xstab_anc[0:num_ancx] + sub_z = zstab_anc[0:num_ancz] custom_memory_circuit_stabs( sub_data, sub_x, sub_z, - 1, + num_rounds + 1, cnot_schedX_flat, cnot_schedZ_flat, - allow_device_calls, - False, + False, # enqueue_synd + True, # do_errors_after_non_last_rounds p_spam, - i, - decoder_window, - manually_inject_errors, + 0, + num_rounds + 1, # decoder_window: a single window + False, # manually_inject_errors + True, # declare_detectors ) + else: + # One stabilizer round to lock in + for i in range(num_logical): + sub_data = data[i * num_data:(i + 1) * + num_data] # FIXME: all sub_data are incorrect + sub_x = xstab_anc[i * num_ancx:(i + 1) * + num_ancx] # same other vectors + sub_z = zstab_anc[i * num_ancz:(i + 1) * num_ancz] + custom_memory_circuit_stabs( + sub_data, + sub_x, + sub_z, + 1, + cnot_schedX_flat, + cnot_schedZ_flat, + allow_device_calls, + False, + p_spam, + i, + decoder_window, + manually_inject_errors, + False, # declare_detectors + ) - # Inject errors - for i in range(num_logical): - sub_data = data[i * num_data:(i + 1) * - num_data] # FIXME: all sub_data are incorrect - sub_x = xstab_anc[i * num_ancx:(i + 1) * num_ancx] # same other vectors - sub_z = zstab_anc[i * num_ancz:(i + 1) * num_ancz] - logical = patch(sub_data, sub_x, sub_z) - spam_error(logical, p_spam, 0.0, 0.0) + # Inject errors + for i in range(num_logical): + sub_data = data[i * num_data:(i + 1) * + num_data] # FIXME: all sub_data are incorrect + sub_x = xstab_anc[i * num_ancx:(i + 1) * + num_ancx] # same other vectors + sub_z = zstab_anc[i * num_ancz:(i + 1) * num_ancz] + logical = patch(sub_data, sub_x, sub_z) + spam_error(logical, p_spam, 0.0, 0.0) - # Do stabilizer rounds - for i in range(num_logical): - sub_data = data[i * num_data:(i + 1) * - num_data] # FIXME: all sub_data are incorrect - sub_x = xstab_anc[i * num_ancx:(i + 1) * num_ancx] # same other vectors - sub_z = zstab_anc[i * num_ancz:(i + 1) * num_ancz] - custom_memory_circuit_stabs( - sub_data, - sub_x, - sub_z, - num_rounds, - cnot_schedX_flat, - cnot_schedZ_flat, - allow_device_calls, - True, - p_spam, - i, - decoder_window, - manually_inject_errors, - ) + # Do stabilizer rounds + for i in range(num_logical): + sub_data = data[i * num_data:(i + 1) * + num_data] # FIXME: all sub_data are incorrect + sub_x = xstab_anc[i * num_ancx:(i + 1) * + num_ancx] # same other vectors + sub_z = zstab_anc[i * num_ancz:(i + 1) * num_ancz] + custom_memory_circuit_stabs( + sub_data, + sub_x, + sub_z, + num_rounds, + cnot_schedX_flat, + cnot_schedZ_flat, + allow_device_calls, + True, + p_spam, + i, + decoder_window, + manually_inject_errors, + False, # declare_detectors + ) # Only apply corrections after processing all windows if allow_device_calls and apply_corrections: @@ -395,7 +461,19 @@ def demo_circuit_qpu( ret = ret << num_data sub_data = data[i * num_data:(i + 1) * num_data] sub_meas = mz(sub_data) - ret |= cudaq.to_integer(cudaq.to_bools(sub_meas)) + if declare_detectors and i == 0: + zlog = [cudaq.measure_handle() for _ in range(len(z_obs_indices))] + for k in range(len(z_obs_indices)): + zlog[k] = sub_meas[z_obs_indices[k]] + cudaq.logical_observable(zlog, observable_index=0) + # Pack the measured bits branch-free (bit j = data qubit j, the same + # LSB-first order as cudaq.to_integer). Routing the measurement + # results through a call (to_bools/to_integer) or a branch would tag + # this kernel with qubitMeasurementFeedback, which + # cudaq.dem_from_kernel rejects. + for j in range(num_data): + bitval = sub_meas[j] + ret = ret | (bitval << j) # The remaining bits are allocated to the number of corrections. ret = ret | (num_corrections << (num_data * num_logical)) @@ -453,7 +531,16 @@ def demo_circuit_host(code_obj: qec.code, noise = cudaq.NoiseModel() - # Build or load DEM (MSM path) + # The Z logical observable's data-qubit support (row 0 of the Z + # observables matrix); demo_circuit_qpu declares the matching + # logical_observable over the final data measurements when generating the + # DEM. + obs_matrix = code_obj.get_observables_z() + z_obs_indices = [ + int(col) for col in range(obs_matrix.shape[1]) if obs_matrix[0, col] + ] + + # Build or load DEM dem = qec.DetectorErrorModel() if load_dem: @@ -463,101 +550,59 @@ def demo_circuit_host(code_obj: qec.code, print(f"Preparing DEM to save to {dem_filename}") # Always use stim to build the DEM cudaq.set_target("stim") - cudaq.set_noise(noise) if p_spam == 0.0: - raise RuntimeError( - "Cannot build a DEM with p_spam = 0.0 (cannot get the MSM).") - # Always use numLogical = 1 for the MSM - (msm_as_strings, msm_dimensions, msm_probabilities, msm_prob_err_id - ) = qec.compute_msm( - lambda: demo_circuit_qpu( - False, - num_data, - num_ancx, - num_ancz, - decoder_window, # Use decoder_window instead of numRounds for DEM generation - 1, # numLogical - cnot_schedX_flat, - cnot_schedZ_flat, - p_spam, - False, # applyCorrections - decoder_window, - False, # manuallyInjectErrors - ), - True) - - print("MSM result obtained.") - # print(f"MSM dimensions: {msm_dimensions}") - # print(f"MSM probabilities: {msm_probabilities}") - # print(f"MSM probability error ID: {msm_prob_err_id}") - - # Populate error rates and error IDs - dem.error_rates = msm_probabilities - dem.error_ids = msm_prob_err_id - mzTable = qec.construct_mz_table(msm_as_strings) - print("mzTable:", mzTable) - # Subtract the number of data qubits to get the number of syndrome measurements. - totalNumSyndromes = mzTable.shape[0] - distance * distance - numNoiseMechs = mzTable.shape[1] - numSyndromesPerRound = distance * distance - 1 - if (totalNumSyndromes % numSyndromesPerRound != 0): - raise RuntimeError("Num syndromes per round is not a divisor of " - "the number of syndrome measurements") - - numRoundsOfSyndromData = totalNumSyndromes // numSyndromesPerRound - if (numRoundsOfSyndromData != decoder_window + - 1): # Use decoder_window instead of numRounds - raise RuntimeError("Num rounds of syndrome data [" + - str(numRoundsOfSyndromData) + - "] is not equal to the decoder_window + 1[" + - str(decoder_window + 1) + "]") - detector_error_matrix = np.zeros( - (decoder_window * numSyndromesPerRound, numNoiseMechs), - dtype=np.uint8) - # There should be (decoder_window + 1) rounds of data in MSM. - # TODO: [feature] Good candidate. Auto-generating the detector error - # matrix. Currently, we need to manually construct the detector error - # matrix by copying the measurements from the MSM. - for round in range( - decoder_window): # Use decoder_window instead of numRounds - for syndrome in range(numSyndromesPerRound): - for noise_mech in range(numNoiseMechs): - detector_error_matrix[ - round * numSyndromesPerRound + syndrome, - noise_mech] = mzTable[ - (round + 0) * numSyndromesPerRound + syndrome, - noise_mech] ^ mzTable[ - (round + 1) * numSyndromesPerRound + syndrome, - noise_mech] - dem.detector_error_matrix = detector_error_matrix - print("detector_error_matrix:", dem.detector_error_matrix) - - first_data_row = ( - decoder_window + - 1) * numSyndromesPerRound # Use decoder_window instead of numRounds + raise RuntimeError("Cannot build a DEM with p_spam = 0.0.") + # Analyze the same demo_circuit_qpu kernel the shots run, with + # declare_detectors so it annotates its detectors and observable. + # Always use numLogical = 1, and decoder_window rounds instead of + # numRounds: the decoder consumes one window at a time. + dem_text, m2d, m2o = cudaq.dem_from_kernel( + demo_circuit_qpu, + False, # allow_device_calls + True, # declare_detectors + num_data, + num_ancx, + num_ancz, + decoder_window, # Use decoder_window instead of numRounds for DEM generation + 1, # numLogical + cnot_schedX_flat, + cnot_schedZ_flat, + p_spam, + False, # applyCorrections + decoder_window, + False, # manuallyInjectErrors + z_obs_indices, + noise_model=noise, + return_measurement_matrices=True) + dem = qec.dem_from_stim_text(dem_text) + numSyndromesPerRound = distance * distance - 1 - msm_obs = np.zeros((mzTable.shape[0] - first_data_row, numNoiseMechs), - dtype=np.uint8) - for row in range(first_data_row, mzTable.shape[0]): - for col in range(numNoiseMechs): - msm_obs[row - first_data_row, col] = mzTable[row, col] - - print("msm_obs:", msm_obs) - - # Populate dem.observables_flips_matrix by converting the physical data - # qubit measurements to logical observables. - obs_matrix = code_obj.get_observables_z() - print("obs_matrix:", obs_matrix) - dem.observables_flips_matrix = (obs_matrix @ msm_obs) % 2 + if dem.num_detectors() != decoder_window * numSyndromesPerRound: + raise RuntimeError( + "Number of detectors [" + str(dem.num_detectors()) + + "] is not equal to decoder_window * numSyndromesPerRound [" + + str(decoder_window * numSyndromesPerRound) + "]") print("numSyndromesPerRound:", numSyndromesPerRound) dem.canonicalize_for_rounds(numSyndromesPerRound, remove_zero_syndrome_errors=True) + # The runtime detector matrix comes straight from the analysis' + # measurements-to-detectors map: row d lists the (chronological, and + # thus enqueue-ordered) measurement indices whose XOR forms detector + # d. + m2d = m2d.tocsr() + d_sparse = [] + for r in range(m2d.shape[0]): + row = m2d.indices[m2d.indptr[r]:m2d.indptr[r + 1]] + d_sparse.extend(sorted(int(c) for c in row)) + d_sparse.append(-1) + print("dem.detector_error_matrix:") print(dem.detector_error_matrix) print("dem.observables_flips_matrix:") print(dem.observables_flips_matrix) - save_dem_to_file(dem, dem_filename, numSyndromesPerRound, num_logical) + save_dem_to_file(dem, d_sparse, dem_filename, numSyndromesPerRound, + num_logical) return # Actual run @@ -614,7 +659,8 @@ def demo_circuit_host(code_obj: qec.code, # Run shots run_result = cudaq.run( demo_circuit_qpu, - True, + True, # allow_device_calls + False, # declare_detectors # prep_0, num_data, num_ancx, @@ -627,13 +673,13 @@ def demo_circuit_host(code_obj: qec.code, True, decoder_window, manually_inject_errors, + z_obs_indices, shots_count=num_shots, noise_model=cudaq.NoiseModel() if not is_remote_qpu else None) print("Done with cudaq.run!") # print(f"Result: {len(run_result)}") - obs_matrix = code_obj.get_observables_z() num_non_zero = 0 num_corrections = 0 print("Result size: " + str(len(run_result)))