diff --git a/BUILD b/BUILD index f82db683..a7849c35 100644 --- a/BUILD +++ b/BUILD @@ -19,6 +19,7 @@ py_wheel( deps=[ "//src:tesseract_decoder", "//src/py:generated_stubs", + "//src/py:multi_pass_sinter_decoders", "//src/py/_tesseract_py_util:_tesseract_py_util", ":package_data", ], @@ -26,6 +27,7 @@ py_wheel( requires=[ "numpy", "scipy", + "sinter", "stim", ], python_tag="$(TARGET_VERSION)", @@ -51,3 +53,9 @@ config_setting( "@platforms//cpu:arm64", ], ) + +filegroup( + name = "testdata", + srcs = glob(["testdata/**/*"]), + visibility = ["//visibility:public"], +) diff --git a/CMakeLists.txt b/CMakeLists.txt index 23e560f4..0ffd445a 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -93,6 +93,26 @@ target_include_directories(visualization PUBLIC ${TESSERACT_SRC_DIR}) target_compile_options(visualization PRIVATE ${OPT_COPTS}) target_link_libraries(visualization PUBLIC common boost_headers) +add_library(error_correlations ${TESSERACT_SRC_DIR}/multi_pass/error_correlations.cc ${TESSERACT_SRC_DIR}/multi_pass/error_correlations.h) +target_include_directories(error_correlations PUBLIC ${TESSERACT_SRC_DIR}) +target_compile_options(error_correlations PRIVATE ${OPT_COPTS}) +target_link_libraries(error_correlations PUBLIC libstim) + +add_library(tanner_graph ${TESSERACT_SRC_DIR}/multi_pass/tanner_graph.cc ${TESSERACT_SRC_DIR}/multi_pass/tanner_graph.h) +target_include_directories(tanner_graph PUBLIC ${TESSERACT_SRC_DIR}) +target_compile_options(tanner_graph PRIVATE ${OPT_COPTS}) +target_link_libraries(tanner_graph PUBLIC libstim) + +add_library(dem_decomposition ${TESSERACT_SRC_DIR}/multi_pass/dem_decomposition.cc ${TESSERACT_SRC_DIR}/multi_pass/dem_decomposition.h) +target_include_directories(dem_decomposition PUBLIC ${TESSERACT_SRC_DIR}) +target_compile_options(dem_decomposition PRIVATE ${OPT_COPTS}) +target_link_libraries(dem_decomposition PUBLIC libstim) + +add_library(multi_pass_tesseract_decoder ${TESSERACT_SRC_DIR}/multi_pass/multi_pass_tesseract_decoder.cc ${TESSERACT_SRC_DIR}/multi_pass/multi_pass_tesseract_decoder.h) +target_include_directories(multi_pass_tesseract_decoder PUBLIC ${TESSERACT_SRC_DIR}) +target_compile_options(multi_pass_tesseract_decoder PRIVATE ${OPT_COPTS}) +target_link_libraries(multi_pass_tesseract_decoder PUBLIC common tesseract_lib tanner_graph error_correlations dem_decomposition libstim) + add_library(tesseract_lib ${TESSERACT_SRC_DIR}/tesseract.cc ${TESSERACT_SRC_DIR}/tesseract.h) target_include_directories(tesseract_lib PUBLIC ${TESSERACT_SRC_DIR}) target_compile_options(tesseract_lib PRIVATE ${OPT_COPTS}) @@ -114,7 +134,7 @@ target_link_libraries(simplex PUBLIC common utils tesseract_lib highs libstim Th # === Executables === add_executable(tesseract ${TESSERACT_SRC_DIR}/tesseract_main.cc) target_compile_options(tesseract PRIVATE ${OPT_COPTS}) -target_link_libraries(tesseract PRIVATE tesseract_lib argparse::argparse nlohmann_json::nlohmann_json) +target_link_libraries(tesseract PRIVATE tesseract_lib multi_pass_tesseract_decoder argparse::argparse nlohmann_json::nlohmann_json) add_executable(tesseract_trellis ${TESSERACT_SRC_DIR}/tesseract_trellis_main.cc) target_compile_options(tesseract_trellis PRIVATE ${OPT_COPTS}) @@ -131,7 +151,7 @@ target_link_libraries(simplex_bin PRIVATE common simplex argparse::argparse nloh pybind11_add_module(tesseract_decoder MODULE ${TESSERACT_SRC_DIR}/tesseract.pybind.cc) target_compile_options(tesseract_decoder PRIVATE ${OPT_COPTS}) target_include_directories(tesseract_decoder PRIVATE ${TESSERACT_SRC_DIR}) -target_link_libraries(tesseract_decoder PRIVATE common utils simplex tesseract_lib) +target_link_libraries(tesseract_decoder PRIVATE common utils simplex tesseract_lib multi_pass_tesseract_decoder) set_target_properties(tesseract_decoder PROPERTIES LIBRARY_OUTPUT_DIRECTORY ${PROJECT_SOURCE_DIR}/src LIBRARY_OUTPUT_DIRECTORY_DEBUG ${PROJECT_SOURCE_DIR}/src @@ -154,3 +174,20 @@ add_test(NAME tesseract_test COMMAND tesseract_test) add_executable(tesseract_trellis_test ${TESSERACT_SRC_DIR}/tesseract_trellis.test.cc) target_link_libraries(tesseract_trellis_test PRIVATE tesseract_trellis_lib GTest::gtest_main) add_test(NAME tesseract_trellis_test COMMAND tesseract_trellis_test) + +add_executable(dem_decomposition_test ${TESSERACT_SRC_DIR}/multi_pass/dem_decomposition.test.cc) +target_link_libraries(dem_decomposition_test PRIVATE dem_decomposition GTest::gtest_main libstim) +add_test(NAME dem_decomposition_test COMMAND dem_decomposition_test) + +add_executable(tanner_graph_test ${TESSERACT_SRC_DIR}/multi_pass/tanner_graph.test.cc) +target_link_libraries(tanner_graph_test PRIVATE tanner_graph GTest::gtest_main libstim) +add_test(NAME tanner_graph_test COMMAND tanner_graph_test) + +add_executable(error_correlations_test ${TESSERACT_SRC_DIR}/multi_pass/error_correlations.test.cc) +target_link_libraries(error_correlations_test PRIVATE error_correlations GTest::gtest_main libstim) +add_test(NAME error_correlations_test COMMAND error_correlations_test) + +add_executable(multi_pass_tesseract_decoder_test ${TESSERACT_SRC_DIR}/multi_pass/multi_pass_tesseract_decoder.test.cc) +target_link_libraries(multi_pass_tesseract_decoder_test PRIVATE multi_pass_tesseract_decoder GTest::gtest_main libstim) +add_test(NAME multi_pass_tesseract_decoder_test COMMAND multi_pass_tesseract_decoder_test) +set_tests_properties(multi_pass_tesseract_decoder_test PROPERTIES WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}) diff --git a/README.md b/README.md index 15331c43..ba0d61d4 100644 --- a/README.md +++ b/README.md @@ -195,6 +195,67 @@ errors are not capped by degree. * *DEM usage frequency output*: if `--dem-out` is specified, outputs estimated error frequencies. * *Statistics output*: includes number of shots, errors, low confidence shots, and processing time. +--- + +## Multi-Pass Graph Shattering + +Multi-pass graph shattering partitions a correlated detector error model into two detector +components and decodes the smaller component models separately. With two passes, predictions from +the first pass update error priors used during the final pass. The current implementation requires +exactly two components and accepts one or two passes. + +### Detector classification + +The CLI classifier checks the following detector annotations in order: + +1. A `"measure_basis"` field in the detector's JSON metadata tag, first at the top level and then + under `"md"`. +2. A `"basis"` field in the same locations. +3. A fourth detector coordinate using the Chromobius-style `color + 3 * basis` convention: + values `0`–`2` select component 0 and values `3`–`5` select component 1. + +Metadata basis values must be strings exactly equal to `"X"` or `"Z"`. An invalid metadata value +does not fall back to another field or to coordinates. Multi-pass decoding fails if any detector +cannot be classified or if the resulting classification does not contain exactly two components. +The Python wrapper can instead be given a custom detector classifier. + +### CLI options + +* `--multipass`: Enables multi-pass graph shattering. +* `--num-passes`, `--num_passes`: Selects one or two passes (default: 2). One pass performs no + inter-pass prior update; two passes perform one round of prior propagation. Other values are + rejected. +* `--multipass-strategy`, `--multipass_strategy`: Selects `causal` (default), which derives the pass + schedule from component dependencies, or experimental `static`, which schedules both components + in every pass. +* `--print-multipass-plan`: Prints the monolithic and component model statistics, dependencies, and + pass schedule to standard error. It requires `--multipass`; these statistics are calculated only + when this flag is present. + +`--dem-out` is not supported with `--multipass`. + +### CLI example + +This example uses a coordinate-annotated color-code circuit and the long-beam settings: + +```bash +./bazel-bin/src/tesseract \ + --circuit testdata/colorcodes/r=5,d=5,p=0.003,noise=si1000,c=midout_color_code_X,q=23,gates=cz.stim \ + --sample-num-shots 1000 \ + --multipass \ + --num-passes 2 \ + --multipass-strategy causal \ + --pqlimit 1000000 \ + --beam 20 \ + --beam-climbing \ + --no-revisit-dets \ + --num-det-orders 21 \ + --print-multipass-plan \ + --print-stats +``` + +--- + ## Python Interface [Full Python wrapper documentation](src/py/README.md) diff --git a/src/BUILD b/src/BUILD index 1b1407bf..661f8c21 100644 --- a/src/BUILD +++ b/src/BUILD @@ -89,6 +89,7 @@ pybind_library( "visualization.pybind.h", "tesseract.pybind.h", "tesseract_sinter_compat.pybind.h", + "multi_pass/multi_pass_sinter_compat.pybind.h", ], copts = OPT_COPTS, deps = [ @@ -96,6 +97,7 @@ pybind_library( ":libutils", ":libsimplex", ":libtesseract", + ":libmulti_pass_tesseract_decoder", ], ) @@ -171,6 +173,108 @@ cc_library( ], ) +cc_library( + name = "libmulti_pass_tesseract_decoder", + srcs = ["multi_pass/multi_pass_tesseract_decoder.cc"], + hdrs = ["multi_pass/multi_pass_tesseract_decoder.h"], + copts = OPT_COPTS, + linkopts = OPT_LINKOPTS, + deps = [ + ":libcommon", + ":libtesseract", + ":libtanner_graph", + ":liberror_correlations", + ":libdem_decomposition", + "@stim//:stim_lib", + ], +) + +cc_test( + name = "multi_pass_tesseract_decoder_tests", + srcs = ["multi_pass/multi_pass_tesseract_decoder.test.cc"], + copts = OPT_COPTS, + linkopts = OPT_LINKOPTS, + data = ["//:testdata"], + deps = [ + ":libmulti_pass_tesseract_decoder", + "@gtest", + "@gtest//:gtest_main", + "@stim//:stim_lib", + ], +) + +cc_library( + name = "liberror_correlations", + srcs = ["multi_pass/error_correlations.cc"], + hdrs = ["multi_pass/error_correlations.h"], + copts = OPT_COPTS, + linkopts = OPT_LINKOPTS, + deps = [ + "@stim//:stim_lib", + ], +) + +cc_test( + name = "error_correlations_tests", + srcs = ["multi_pass/error_correlations.test.cc"], + copts = OPT_COPTS, + linkopts = OPT_LINKOPTS, + deps = [ + ":liberror_correlations", + "@gtest", + "@gtest//:gtest_main", + "@stim//:stim_lib", + ], +) + +cc_library( + name = "libtanner_graph", + srcs = ["multi_pass/tanner_graph.cc"], + hdrs = ["multi_pass/tanner_graph.h"], + copts = OPT_COPTS, + linkopts = OPT_LINKOPTS, + deps = [ + "@stim//:stim_lib", + ], +) + +cc_test( + name = "tanner_graph_tests", + srcs = ["multi_pass/tanner_graph.test.cc"], + copts = OPT_COPTS, + linkopts = OPT_LINKOPTS, + deps = [ + ":libtanner_graph", + "@gtest", + "@gtest//:gtest_main", + "@stim//:stim_lib", + ], +) + +cc_library( + name = "libdem_decomposition", + srcs = ["multi_pass/dem_decomposition.cc"], + hdrs = ["multi_pass/dem_decomposition.h"], + copts = OPT_COPTS, + linkopts = OPT_LINKOPTS, + deps = [ + "@stim//:stim_lib", + ], +) + +cc_test( + name = "dem_decomposition_tests", + srcs = ["multi_pass/dem_decomposition.test.cc"], + copts = OPT_COPTS, + linkopts = OPT_LINKOPTS, + deps = [ + ":libdem_decomposition", + "@gtest", + "@gtest//:gtest_main", + "@stim//:stim_lib", + ], +) + cc_binary( name = "tesseract", srcs = ["tesseract_main.cc"], @@ -178,6 +282,7 @@ cc_binary( linkopts = OPT_LINKOPTS, deps = [ ":libtesseract", + ":libmulti_pass_tesseract_decoder", "@argparse", "@nlohmann_json//:json", "@stim//:stim_lib", diff --git a/src/multi_pass/dem_decomposition.cc b/src/multi_pass/dem_decomposition.cc new file mode 100644 index 00000000..db19dcd0 --- /dev/null +++ b/src/multi_pass/dem_decomposition.cc @@ -0,0 +1,299 @@ +#include "dem_decomposition.h" + +#include +#include +#include +#include +#include +#include +#include + +#include "stim.h" + +namespace tesseract { + +// Helper function to generate all combinations of observables +void generate_obs_combinations( + const std::vector>>& obs_options_by_component, + std::vector>& current_combination, + std::vector>>& all_combinations, int component_index) { + if (component_index == (int)obs_options_by_component.size()) { + all_combinations.push_back(current_combination); + return; + } + + for (const auto& obs_option : obs_options_by_component[component_index]) { + current_combination.push_back(obs_option); + generate_obs_combinations(obs_options_by_component, current_combination, all_combinations, + component_index + 1); + current_combination.pop_back(); + } +} + +std::vector reduce_symmetric_difference(const std::vector& items) { + std::set unpaired_set; + for (int item : items) { + if (unpaired_set.count(item)) { + unpaired_set.erase(item); + } else { + unpaired_set.insert(item); + } + } + return std::vector(unpaired_set.begin(), unpaired_set.end()); +} + +std::vector reduce_set_symmetric_difference(const std::vector>& sets) { + std::vector all_items; + for (const auto& s : sets) { + all_items.insert(all_items.end(), s.begin(), s.end()); + } + return reduce_symmetric_difference(all_items); +} + +std::pair, std::vector> undecomposed_error_detectors_and_observables( + const stim::DemInstruction& instruction) { + if (instruction.type != stim::DemInstructionType::DEM_ERROR) { + throw std::invalid_argument("DEM instruction must be an error"); + } + + std::vector detectors; + std::vector observables; + for (const auto& target : instruction.target_data) { + if (target.is_relative_detector_id()) { + detectors.push_back(target.val()); + } else if (target.is_observable_id()) { + observables.push_back(target.val()); + } + } + + return {reduce_symmetric_difference(detectors), reduce_symmetric_difference(observables)}; +} + +std::vector> get_component_obs_matching_undecomposed_obs( + const std::vector>>& obs_options_by_component, + const std::vector& error_obs, int num_missing_components, bool allow_remnant_errors) { + if (!allow_remnant_errors && num_missing_components > 0) { + return {}; + } + + std::vector>> all_combinations; + std::vector> current_combination; + generate_obs_combinations(obs_options_by_component, current_combination, all_combinations, 0); + + std::vector error_obs_reduced = reduce_symmetric_difference(error_obs); + std::set error_obs_set(error_obs_reduced.begin(), error_obs_reduced.end()); + + for (const auto& combination : all_combinations) { + std::vector known_obs_sum = reduce_set_symmetric_difference(combination); + + // Residual = error_obs XOR known_obs_sum + std::vector residual_input = error_obs_reduced; + residual_input.insert(residual_input.end(), known_obs_sum.begin(), known_obs_sum.end()); + std::vector residual = reduce_symmetric_difference(residual_input); + + if (residual.empty()) { + // Case A: Residual is empty. All missing components get no observables. + std::vector> result = combination; + for (int i = 0; i < num_missing_components; ++i) result.push_back({}); + return result; + } + + if (num_missing_components == 1 && allow_remnant_errors) { + // Case B: Residual is non-empty and one component is missing. + // Assign the entire residual to the missing component. + std::vector> result = combination; + result.push_back(residual); + return result; + } + } + + return {}; +} + +stim::DetectorErrorModel decompose_errors_using_detector_assignment( + const stim::DetectorErrorModel& dem, const std::function& detector_component_func, + bool allow_remnant_errors) { + stim::DetectorErrorModel flattened_dem = dem.flattened(); + std::map, std::set>> single_component_dets_to_obs; + + for (const auto& instruction : flattened_dem.instructions) { + if (instruction.type != stim::DemInstructionType::DEM_ERROR) continue; + + auto [detectors, observables] = undecomposed_error_detectors_and_observables(instruction); + + std::unordered_set components; + for (int d : detectors) components.insert(detector_component_func(d)); + + if (components.size() <= 1) { + single_component_dets_to_obs[detectors].insert(observables); + } + } + + stim::DetectorErrorModel output_dem; + for (const auto& instruction : flattened_dem.instructions) { + if (instruction.type != stim::DemInstructionType::DEM_ERROR) { + output_dem.append_dem_instruction(instruction); + continue; + } + + auto [detectors, observables] = undecomposed_error_detectors_and_observables(instruction); + + std::map> dets_by_comp_id; + std::set unique_components; + for (int d : detectors) { + int c = detector_component_func(d); + dets_by_comp_id[c].push_back(d); + unique_components.insert(c); + } + + std::vector> dets_by_component; + std::vector>> obs_options_by_known_component; + std::vector> missing_components_dets; + + for (int c : unique_components) { + std::vector component_dets = dets_by_comp_id[c]; + std::sort(component_dets.begin(), component_dets.end()); + + if (single_component_dets_to_obs.count(component_dets)) { + dets_by_component.push_back(component_dets); + obs_options_by_known_component.push_back(single_component_dets_to_obs[component_dets]); + } else { + if (!allow_remnant_errors) { + throw std::invalid_argument( + "Component not present as its own error and allow_remnant_errors=false"); + } + missing_components_dets.push_back(component_dets); + } + } + + std::vector> consistent_obs_by_component = + get_component_obs_matching_undecomposed_obs(obs_options_by_known_component, observables, + (int)missing_components_dets.size(), + allow_remnant_errors); + + if (consistent_obs_by_component.empty()) { + throw std::invalid_argument("Error instruction could not be decomposed consistently."); + } + + std::vector targets; + std::vector> all_dets = dets_by_component; + all_dets.insert(all_dets.end(), missing_components_dets.begin(), missing_components_dets.end()); + + for (size_t i = 0; i < all_dets.size(); ++i) { + for (int d : all_dets[i]) targets.push_back(stim::DemTarget::relative_detector_id(d)); + for (int o : consistent_obs_by_component[i]) + targets.push_back(stim::DemTarget::observable_id(o)); + if (i != all_dets.size() - 1) targets.push_back(stim::DemTarget::separator()); + } + + output_dem.append_error_instruction(instruction.arg_data[0], targets, instruction.tag); + } + return output_dem; +} + +stim::DetectorErrorModel decompose_errors_using_generic_classifier( + const stim::DetectorErrorModel& dem, const DetectorClassifier& classifier, + bool allow_remnant_errors) { + stim::DetectorErrorModel flattened = dem.flattened(); + + std::set all_detector_indices; + std::map detector_tags; + for (uint64_t d = 0; d < flattened.count_detectors(); ++d) { + all_detector_indices.insert(d); + } + for (const auto& inst : flattened.instructions) { + if (inst.type == stim::DemInstructionType::DEM_DETECTOR) { + int d = inst.target_data[0].val(); + detector_tags[d] = inst.tag; + } + } + + auto detector_coords = flattened.get_detector_coordinates(all_detector_indices); + + std::vector classification_cache(flattened.count_detectors()); + for (uint64_t d : all_detector_indices) { + std::vector coords = + detector_coords.count(d) ? detector_coords.at(d) : std::vector{}; + classification_cache[d] = classifier((int)d, coords, detector_tags[d]); + } + + auto component_func = [&](int d) { + if (d < 0 || (size_t)d >= classification_cache.size()) { + throw std::invalid_argument("Detector D" + std::to_string(d) + " is out of range."); + } + return classification_cache[d]; + }; + + return decompose_errors_using_detector_assignment(flattened, component_func, + allow_remnant_errors); +} + +std::map split_dem_by_component( + const stim::DetectorErrorModel& dem, const std::function& detector_component_func) { + std::map component_dems; + + for (const auto& instruction : dem.instructions) { + if (instruction.type == stim::DemInstructionType::DEM_ERROR) { + double prob = instruction.arg_data[0]; + + size_t group_start = 0; + for (size_t k = 0; k <= instruction.target_data.size(); ++k) { + if (k == instruction.target_data.size() || instruction.target_data[k].is_separator()) { + std::vector component_targets; + std::set component_ids; + for (size_t j = group_start; j < k; ++j) { + const auto& target = instruction.target_data[j]; + component_targets.push_back(target); + if (target.is_relative_detector_id()) { + component_ids.insert(detector_component_func(target.val())); + } + } + + if (component_ids.empty()) { + // If no detectors, we can't assign it to a component based on detectors. + // For now, let's skip or handle separately. + } else if (component_ids.size() > 1) { + throw std::invalid_argument("Mixed component ID in a single error component group."); + } else { + int comp_id = *component_ids.begin(); + component_dems[comp_id].append_error_instruction(prob, component_targets, ""); + } + group_start = k + 1; + } + } + } else if (instruction.type == stim::DemInstructionType::DEM_DETECTOR || + instruction.type == stim::DemInstructionType::DEM_LOGICAL_OBSERVABLE) { + for (auto& pair : component_dems) { + pair.second.append_dem_instruction(instruction); + } + } + } + return component_dems; +} + +stim::DetectorErrorModel undecompose_errors(const stim::DetectorErrorModel& dem) { + stim::DetectorErrorModel undecomposed_dem; + for (const auto& instruction : dem.instructions) { + if (instruction.type == stim::DemInstructionType::DEM_REPEAT_BLOCK) { + undecomposed_dem.append_repeat_block(instruction.repeat_block_rep_count(), + undecompose_errors(instruction.repeat_block_body(dem)), + instruction.tag); + continue; + } + + if (instruction.type != stim::DemInstructionType::DEM_ERROR) { + undecomposed_dem.append_dem_instruction(instruction); + continue; + } + + auto [detectors, observables] = undecomposed_error_detectors_and_observables(instruction); + std::vector targets; + for (int d : detectors) targets.push_back(stim::DemTarget::relative_detector_id(d)); + for (int o : observables) targets.push_back(stim::DemTarget::observable_id(o)); + + undecomposed_dem.append_error_instruction(instruction.arg_data[0], targets, instruction.tag); + } + return undecomposed_dem; +} + +} // namespace tesseract diff --git a/src/multi_pass/dem_decomposition.h b/src/multi_pass/dem_decomposition.h new file mode 100644 index 00000000..94dd93b3 --- /dev/null +++ b/src/multi_pass/dem_decomposition.h @@ -0,0 +1,81 @@ +#ifndef DEM_DECOMPOSITION_H +#define DEM_DECOMPOSITION_H + +#include +#include +#include +#include +#include +#include +#include +#include + +#include "stim.h" + +namespace tesseract { + +// Calculates the symmetric difference of a multiset of items. +// Returns items that appear an odd number of times in the input. +std::vector reduce_symmetric_difference(const std::vector& items); + +// Calculates the symmetric difference of a multiset of items given as a vector of sets. +std::vector reduce_set_symmetric_difference(const std::vector>& sets); + +// Extracts detector and observable indices from a Stim error instruction, +// handling decomposed errors by taking the symmetric difference. +std::pair, std::vector> undecomposed_error_detectors_and_observables( + const stim::DemInstruction& instruction); + +/** + * Given possible observables for each component and the error's observables, + * finds a consistent assignment of observables to components. + * + * @param obs_options_by_component A list of sets, where each set contains the possible + * observable flip combinations for a component. + * @param error_obs The total logical observables flipped by the undecomposed error. + * @param num_missing_components Number of components that were not found in the DEM. + * @param allow_remnant_errors If true, allow one component missing from the DEM to be assigned + * residual observables. Multiple missing components are ambiguous. + */ +std::vector> get_component_obs_matching_undecomposed_obs( + const std::vector>>& obs_options_by_component, + const std::vector& error_obs, int num_missing_components = 0, + bool allow_remnant_errors = false); + +/** + * Decomposes errors in a DetectorErrorModel based on detector assignments to components. + * + * @param dem The input DetectorErrorModel. + * @param detector_component_func A function that maps a detector ID to a component ID (int). + * @param allow_remnant_errors If true, allow the decomposition to infer observables for one + * component missing from the DEM. + */ +stim::DetectorErrorModel decompose_errors_using_detector_assignment( + const stim::DetectorErrorModel& dem, const std::function& detector_component_func, + bool allow_remnant_errors = false); + +/** + * A generic classifier that receives full metadata for a detector. + */ +using DetectorClassifier = + std::function& coords, const std::string& tag)>; + +/** + * Decomposes errors using a generic classifier that can look at index, coordinates, and tags. + */ +stim::DetectorErrorModel decompose_errors_using_generic_classifier( + const stim::DetectorErrorModel& dem, const DetectorClassifier& classifier, + bool allow_remnant_errors = false); + +/** + * Splits a decomposed DEM into separate DEMs, one for each component ID. + */ +std::map split_dem_by_component( + const stim::DetectorErrorModel& dem, const std::function& detector_component_func); + +// Returns a detector error model with any error decompositions removed. +stim::DetectorErrorModel undecompose_errors(const stim::DetectorErrorModel& dem); + +} // namespace tesseract + +#endif // DEM_DECOMPOSITION_H diff --git a/src/multi_pass/dem_decomposition.test.cc b/src/multi_pass/dem_decomposition.test.cc new file mode 100644 index 00000000..6cd54469 --- /dev/null +++ b/src/multi_pass/dem_decomposition.test.cc @@ -0,0 +1,201 @@ +#include "dem_decomposition.h" + +#include +#include +#include + +#include "gtest/gtest.h" + +using namespace tesseract; + +TEST(DemDecompositionTest, ReduceSymmetricDifference) { + ASSERT_EQ(reduce_symmetric_difference({1, 2, 3}), std::vector({1, 2, 3})); + ASSERT_EQ(reduce_symmetric_difference({1, 1}), std::vector({})); + ASSERT_EQ(reduce_symmetric_difference({3, 0, 1, 4, 1, 2, 4}), std::vector({0, 2, 3})); +} + +TEST(DemDecompositionTest, ReduceSetSymmetricDifference) { + ASSERT_EQ(reduce_set_symmetric_difference({{1, 2, 3}, {2, 4, 0}}), + std::vector({0, 1, 3, 4})); + ASSERT_EQ(reduce_set_symmetric_difference({{}, {}}), std::vector({})); +} + +TEST(DemDecompositionTest, GetComponentObsMatchingUndecomposedObs) { + std::vector>> component_obs = {{{0, 1}, {2, 1}}, {{3, 4}, {10, 0}}}; + std::vector error_obs = {1, 10}; + std::vector> expected_output = {{0, 1}, {10, 0}}; + ASSERT_EQ(get_component_obs_matching_undecomposed_obs(component_obs, error_obs, 0, false), + expected_output); + + component_obs = {{{}}, {{}}}; + error_obs = {}; + expected_output = {{}, {}}; + ASSERT_EQ(get_component_obs_matching_undecomposed_obs(component_obs, error_obs, 0, false), + expected_output); + + component_obs = {{{}}, {{}}}; + error_obs = {0}; + expected_output = {}; + ASSERT_EQ(get_component_obs_matching_undecomposed_obs(component_obs, error_obs, 0, false), + expected_output); +} + +TEST(DemDecompositionTest, RemnantErrorsSingleMissingComponent) { + std::vector>> component_obs = {{{1}}}; + std::vector error_obs = {1, 2}; + std::vector> expected_output = {{1}, {2}}; + ASSERT_EQ(get_component_obs_matching_undecomposed_obs(component_obs, error_obs, 1, true), + expected_output); +} + +TEST(DemDecompositionTest, RemnantErrorsNoKnownComponents) { + std::vector>> component_obs = {}; + std::vector error_obs = {1, 2}; + std::vector> expected_output = {{1, 2}}; + ASSERT_EQ(get_component_obs_matching_undecomposed_obs(component_obs, error_obs, 1, true), + expected_output); + ASSERT_TRUE( + get_component_obs_matching_undecomposed_obs(component_obs, error_obs, 2, true).empty()); +} + +TEST(DemDecompositionTest, RejectsInconsistentObservableDecomposition) { + stim::DetectorErrorModel dem(R"DEM( + error(0.1) D0 L0 + error(0.1) D1 L1 + error(0.2) D0 D1 L0 + detector D0 + detector D1 + logical_observable L0 + logical_observable L1 + )DEM"); + auto component = [](int detector) { return detector; }; + EXPECT_THROW(decompose_errors_using_detector_assignment(dem, component, true), + std::invalid_argument); +} + +TEST(DemDecompositionTest, DecomposeErrorsUsingGenericClassifier) { + stim::DetectorErrorModel dem(R"DEM( + error(0.1) D0 ^ D1 L1 + error(0.01) D0 D3 D3 D1 L5 L4 L4 + error(0.3) D0 D1 D3 D3 D2 D3 L0 L5 + error(0.2) D3 D2 D0 D0 L0 + detector(0) D0 + detector(0) D1 + detector(1) D2 + detector(1) D3 + )DEM"); + + // Classifier based on coordinate + auto classifier = [](int index, const std::vector& coords, + const std::string& tag) -> int { + if (coords.empty()) return 0; + return (int)coords.back(); + }; + + stim::DetectorErrorModel expected_decomposed_dem(R"DEM( + error(0.1) D0 D1 L1 + error(0.01) D0 D1 L5 + error(0.3) D0 D1 L5 ^ D2 D3 L0 + error(0.2) D2 D3 L0 + detector(0) D0 + detector(0) D1 + detector(1) D2 + detector(1) D3 + )DEM"); + ASSERT_EQ(decompose_errors_using_generic_classifier(dem, classifier).str(), + expected_decomposed_dem.str()); +} + +TEST(DemDecompositionTest, DecomposeErrorsUsingGenericClassifierTagBased) { + stim::DetectorErrorModel dem(R"DEM( + error(0.1) D0 D1 + error(0.2) D2 D3 + error(0.3) D0 D2 + error(0.01) D0 + error(0.01) D2 + detector[{"basis": "X"}] D0 + detector[{"basis": "X"}] D1 + detector[{"basis": "Z"}] D2 + detector[{"basis": "Z"}] D3 + )DEM"); + + // Classifier based on finding "X" or "Z" in the tag + auto classifier = [](int index, const std::vector& coords, + const std::string& tag) -> int { + if (tag.find("\"X\"") != std::string::npos) return 0; + if (tag.find("\"Z\"") != std::string::npos) return 1; + return 2; + }; + + stim::DetectorErrorModel decomposed = decompose_errors_using_generic_classifier(dem, classifier); + + bool found_d0d2_decomposed = false; + for (const auto& inst : decomposed.flattened().instructions) { + if (inst.type == stim::DemInstructionType::DEM_ERROR && inst.arg_data[0] == 0.3) { + bool has_separator = false; + for (const auto& target : inst.target_data) { + if (target.is_separator()) { + has_separator = true; + break; + } + } + if (has_separator) { + found_d0d2_decomposed = true; + } + } + } + ASSERT_TRUE(found_d0d2_decomposed); +} + +TEST(DemDecompositionTest, SplitDemByComponent) { + stim::DetectorErrorModel dem(R"DEM( + error(0.1) D0 D1 + error(0.2) D2 D3 + error(0.3) D0 D2 L0 + error(0.01) D0 + error(0.01) D2 L0 + detector D0 + detector D1 + detector D2 + detector D3 + logical_observable L0 + )DEM"); + + auto classifier = [](int index, const std::vector& coords, + const std::string& tag) -> int { + return (index < 2) ? 0 : 1; // 0,1 -> comp 0; 2,3 -> comp 1 + }; + + stim::DetectorErrorModel decomposed = decompose_errors_using_generic_classifier(dem, classifier); + + auto comp_func = [](int id) { return (id < 2) ? 0 : 1; }; + auto dems = split_dem_by_component(decomposed, comp_func); + + ASSERT_EQ(dems.size(), 2); + ASSERT_EQ(dems[0].count_errors(), 3); + ASSERT_EQ(dems[1].count_errors(), 3); +} + +TEST(DemDecompositionTest, UndecomposeErrorsWithRepeatBlock) { + stim::DetectorErrorModel dem(R"DEM( + error(0.1) D2 D5 ^ D10 L1 + repeat 10 { + error(0.4) D1 L2 L3 ^ D2 ^ D2 L2 + repeat 3 { + error(0.3) D10 D11 ^ D12 + } + } + error(0.5) D0 D100 + )DEM"); + stim::DetectorErrorModel expected_undecomposed_dem(R"DEM( + error(0.1) D2 D5 D10 L1 + repeat 10 { + error(0.4) D1 L3 + repeat 3 { + error(0.3) D10 D11 D12 + } + } + error(0.5) D0 D100 + )DEM"); + ASSERT_EQ(undecompose_errors(dem).str(), expected_undecomposed_dem.str()); +} diff --git a/src/multi_pass/error_correlations.cc b/src/multi_pass/error_correlations.cc new file mode 100644 index 00000000..c3348590 --- /dev/null +++ b/src/multi_pass/error_correlations.cc @@ -0,0 +1,151 @@ +#include "error_correlations.h" + +#include +#include + +namespace tesseract { + +bool ComponentSymptom::operator==(const ComponentSymptom& other) const { + return detectors == other.detectors && observables == other.observables; +} + +bool ComponentSymptom::operator<(const ComponentSymptom& other) const { + if (detectors != other.detectors) return detectors < other.detectors; + return observables < other.observables; +} + +std::string ImpliedProbability::str() const { + std::stringstream ss; + ss << "ImpliedProbability(detectors={"; + for (size_t i = 0; i < affected_symptom.detectors.size(); ++i) { + ss << affected_symptom.detectors[i] << (i == affected_symptom.detectors.size() - 1 ? "" : ","); + } + ss << "}, observables={"; + for (size_t i = 0; i < affected_symptom.observables.size(); ++i) { + ss << affected_symptom.observables[i] + << (i == affected_symptom.observables.size() - 1 ? "" : ","); + } + ss << "}, prob=" << probability << ")"; + return ss.str(); +} + +bool ImpliedProbability::operator==(const ImpliedProbability& other) const { + return affected_symptom == other.affected_symptom && + std::abs(probability - other.probability) < 1e-12; +} + +bool ImpliedProbability::operator<(const ImpliedProbability& other) const { + if (!(affected_symptom == other.affected_symptom)) { + return affected_symptom < other.affected_symptom; + } + return probability < other.probability; +} + +JointProbsMap get_hyperedge_joint_probabilities(const stim::DetectorErrorModel& dem, + const std::vector& global_det_to_comp_id) { + JointProbsMap joint_probs; + auto flattened = dem.flattened(); + + for (const auto& inst : flattened.instructions) { + if (inst.type != stim::DemInstructionType::DEM_ERROR) continue; + + double p = inst.arg_data[0]; + + std::vector components; + inst.for_separated_targets([&](std::span group) { + ComponentSymptom symptom; + int component_id = -1; + for (const auto& target : group) { + if (target.is_relative_detector_id()) { + int detector = target.val(); + if (detector < 0 || (size_t)detector >= global_det_to_comp_id.size() || + global_det_to_comp_id[detector] < 0) { + throw std::invalid_argument("Invalid component assignment for detector D" + + std::to_string(detector) + "."); + } + int detector_component = global_det_to_comp_id[detector]; + if (component_id != -1 && component_id != detector_component) { + throw std::invalid_argument( + "A decomposed error group contains detectors from multiple components."); + } + component_id = detector_component; + symptom.detectors.push_back(detector); + } else if (target.is_observable_id()) { + symptom.observables.push_back(target.val()); + } + } + + if (symptom.detectors.empty()) return; + std::sort(symptom.detectors.begin(), symptom.detectors.end()); + std::sort(symptom.observables.begin(), symptom.observables.end()); + components.push_back(std::move(symptom)); + }); + + // 1. Marginal probabilities (diagonal) + for (const auto& h : components) { + if (joint_probs[h].find(h) == joint_probs[h].end()) { + joint_probs[h][h] = 0.0; + } + // P(A) = P(A) XOR p + joint_probs[h][h] = joint_probs[h][h] * (1 - p) + p * (1 - joint_probs[h][h]); + } + + // 2. Joint probabilities (off-diagonal) + // For a bridging error p connecting A and B, P(A and B) += p (approx) + // Actually, the joint probability is accurately tracked via the same XOR logic + // if we assume independence of other error mechanisms. + if (components.size() > 1) { + for (size_t i = 0; i < components.size(); ++i) { + for (size_t j = 0; j < components.size(); ++j) { + if (i == j) continue; + const auto& hi = components[i]; + const auto& hj = components[j]; + if (joint_probs[hi].find(hj) == joint_probs[hi].end()) { + joint_probs[hi][hj] = 0.0; + } + // For small p, joint probability P(A and B) is roughly the sum of p's of bridging errors + joint_probs[hi][hj] = joint_probs[hi][hj] * (1 - p) + p * (1 - joint_probs[hi][hj]); + } + } + } + } + + return joint_probs; +} + +ImpliedProbsMap get_implied_hyperedge_probabilities(const JointProbsMap& joint_probs) { + ImpliedProbsMap implied_probs; + + for (const auto& [causal, affected_map] : joint_probs) { + double p_causal = 0.0; + auto it_self = affected_map.find(causal); + if (it_self != affected_map.end()) { + p_causal = it_self->second; + } + + if (p_causal <= 0 || p_causal >= 1.0) continue; + + for (const auto& [affected, p_joint] : affected_map) { + if (causal == affected) continue; + + // Conditional Probability P(affected | causal) = P(affected and causal) / P(causal) + double p_conditional = p_joint / p_causal; + + // Cap to 1.0 (numerical precision) + if (p_conditional > 1.0) p_conditional = 1.0; + if (p_conditional < 0.0) p_conditional = 0.0; + + implied_probs[causal].push_back({affected, p_conditional}); + } + } + + return implied_probs; +} + +ImpliedProbsMap process_dem_correlations(const stim::DetectorErrorModel& dem, + const std::vector& global_det_to_comp_id) { + auto joint = get_hyperedge_joint_probabilities(dem, global_det_to_comp_id); + return get_implied_hyperedge_probabilities(joint); +} + +} // namespace tesseract diff --git a/src/multi_pass/error_correlations.h b/src/multi_pass/error_correlations.h new file mode 100644 index 00000000..6de7cca2 --- /dev/null +++ b/src/multi_pass/error_correlations.h @@ -0,0 +1,58 @@ +#ifndef ERROR_CORRELATIONS_H +#define ERROR_CORRELATIONS_H + +#include +#include +#include +#include +#include +#include + +#include "stim.h" + +namespace tesseract { + +struct ComponentSymptom { + std::vector detectors; + std::vector observables; + + bool operator==(const ComponentSymptom& other) const; + bool operator<(const ComponentSymptom& other) const; +}; + +/** + * Represents a probability adjustment for an affected component symptom. + */ +struct ImpliedProbability { + ComponentSymptom affected_symptom; + double probability; // Represents the conditional probability P(affected | causal) + + std::string str() const; + bool operator==(const ImpliedProbability& other) const; + bool operator<(const ImpliedProbability& other) const; +}; + +using JointProbsMap = std::map>; +using ImpliedProbsMap = std::map>; + +/** + * Calculates marginal and joint probabilities for component symptoms in a decomposed DEM. + * Separated groups in one error instruction retain the original physical correlation. + */ +JointProbsMap get_hyperedge_joint_probabilities(const stim::DetectorErrorModel& dem, + const std::vector& global_det_to_comp_id); + +/** + * Calculates conditional probabilities from joint probabilities. + */ +ImpliedProbsMap get_implied_hyperedge_probabilities(const JointProbsMap& joint_probs); + +/** + * Complete workflow for analyzing correlations within a stim::DetectorErrorModel. + */ +ImpliedProbsMap process_dem_correlations(const stim::DetectorErrorModel& dem, + const std::vector& global_det_to_comp_id); + +} // namespace tesseract + +#endif // ERROR_CORRELATIONS_H diff --git a/src/multi_pass/error_correlations.test.cc b/src/multi_pass/error_correlations.test.cc new file mode 100644 index 00000000..a1fddc71 --- /dev/null +++ b/src/multi_pass/error_correlations.test.cc @@ -0,0 +1,65 @@ +#include "error_correlations.h" + +#include + +#include "gtest/gtest.h" + +using namespace tesseract; + +TEST(TwoPassCorrelationsTest, JointProbabilities) { + stim::DetectorErrorModel dem(R"DEM( + error(0.1) D0 ^ D1 L0 + error(0.2) D0 + error(0.05) D1 L1 + )DEM"); + + std::vector global_det_to_comp_id = {0, 1}; + auto joint = get_hyperedge_joint_probabilities(dem, global_det_to_comp_id); + + ComponentSymptom h0{{0}, {}}; + ComponentSymptom h1_l0{{1}, {0}}; + ComponentSymptom h1_l1{{1}, {1}}; + + // P(D0) = 0.1 XOR 0.2 = 0.1*(1-0.2) + 0.2*(1-0.1) = 0.08 + 0.18 = 0.26 + EXPECT_NEAR(joint[h0][h0], 0.26, 1e-6); + // P(D1 L0) = 0.1 + EXPECT_NEAR(joint[h1_l0][h1_l0], 0.1, 1e-6); + EXPECT_NEAR(joint[h1_l1][h1_l1], 0.05, 1e-6); + // P(D0 and D1) = 0.1 + EXPECT_NEAR(joint[h0][h1_l0], 0.1, 1e-6); + EXPECT_NEAR(joint[h1_l0][h0], 0.1, 1e-6); + EXPECT_EQ(joint[h0].count(h1_l1), 0); +} + +TEST(TwoPassCorrelationsTest, ImpliedProbabilities) { + JointProbsMap joint; + ComponentSymptom h0{{0}, {}}; + ComponentSymptom h1{{1}, {0}}; + + joint[h0][h0] = 0.2; + joint[h1][h1] = 0.1; + joint[h0][h1] = 0.05; + joint[h1][h0] = 0.05; + + auto implied = get_implied_hyperedge_probabilities(joint); + + // P(D1 | D0) = 0.05 / 0.2 = 0.25 + bool found = false; + for (const auto& imp : implied[h0]) { + if (imp.affected_symptom == h1) { + EXPECT_NEAR(imp.probability, 0.25, 1e-6); + found = true; + } + } + EXPECT_TRUE(found); + + // P(D0 | D1) = 0.05 / 0.1 = 0.5 + found = false; + for (const auto& imp : implied[h1]) { + if (imp.affected_symptom == h0) { + EXPECT_NEAR(imp.probability, 0.5, 1e-6); + found = true; + } + } + EXPECT_TRUE(found); +} diff --git a/src/multi_pass/multi_pass_sinter_compat.pybind.h b/src/multi_pass/multi_pass_sinter_compat.pybind.h new file mode 100644 index 00000000..feae9380 --- /dev/null +++ b/src/multi_pass/multi_pass_sinter_compat.pybind.h @@ -0,0 +1,162 @@ +#ifndef MULTI_PASS_SINTER_COMPAT_PYBIND_H +#define MULTI_PASS_SINTER_COMPAT_PYBIND_H + +#include +#include +#include +#include +#include + +#include + +#include "../utils.h" +#include "dem_decomposition.h" +#include "multi_pass_tesseract_decoder.h" + +namespace py = pybind11; + +namespace tesseract { + +struct MultiPassSinterCompiledDecoder { + std::unique_ptr decoder; + uint64_t num_detectors; + uint64_t num_observables; + + MultiPassSinterCompiledDecoder(std::unique_ptr d, + uint64_t nd, uint64_t no) + : decoder(std::move(d)), num_detectors(nd), num_observables(no) {} + + size_t num_components() const { + return decoder->num_components(); + } + + py::array_t decode_shots_bit_packed( + const py::array_t& bit_packed_detection_event_data) { + if (bit_packed_detection_event_data.ndim() != 2) + throw std::invalid_argument("Input must be 2D."); + const uint64_t num_detector_bytes = (num_detectors + 7) / 8; + if (bit_packed_detection_event_data.shape(1) != (py::ssize_t)num_detector_bytes) + throw std::invalid_argument("Wrong shape."); + + const size_t num_shots = bit_packed_detection_event_data.shape(0); + const uint64_t num_observable_bytes = (num_observables + 7) / 8; + + auto result_array = + py::array_t({(py::ssize_t)num_shots, (py::ssize_t)num_observable_bytes}); + auto result_buffer = result_array.mutable_data(); + + const uint8_t* detections_data = bit_packed_detection_event_data.data(); + const size_t detections_stride = bit_packed_detection_event_data.strides(0); + + for (size_t shot = 0; shot < num_shots; ++shot) { + const uint8_t* single_shot_data = detections_data + shot * detections_stride; + std::vector detections; + for (uint64_t i = 0; i < num_detectors; ++i) { + if ((single_shot_data[i / 8] >> (i % 8)) & 1) detections.push_back(i); + } + + std::vector predictions = decoder->decode(detections); + uint8_t* single_result_buffer = result_buffer + shot * num_observable_bytes; + std::fill(single_result_buffer, single_result_buffer + num_observable_bytes, 0); + for (int obs_index : predictions) { + if (obs_index >= 0 && (uint64_t)obs_index < num_observables) { + single_result_buffer[obs_index / 8] ^= (1 << (obs_index % 8)); + } + } + } + return result_array; + } +}; + +struct MultiPassSinterDecoder { + size_t num_passes; + py::object full_decomposer; + py::object detector_classifier; + TesseractConfig base_config; + size_t num_det_orders; + ::DetOrder det_order_method; + uint64_t seed; + SchedulingStrategy strategy; + + MultiPassSinterDecoder(size_t n = 2) + : num_passes(n), + full_decomposer(py::none()), + detector_classifier(py::none()), + num_det_orders(1), + det_order_method(::DetOrder::DetIndex), + seed(0), + strategy(SchedulingStrategy::Causal) { + if (num_passes < 1 || num_passes > 2) { + throw std::invalid_argument("num_passes must be 1 or 2."); + } + } + + MultiPassSinterCompiledDecoder compile_decoder_for_dem(const py::object& dem) { + if (detector_classifier.is_none()) { + throw std::invalid_argument("detector_classifier is required for multi-pass decoding."); + } + + stim::DetectorErrorModel stim_dem; + + if (!full_decomposer.is_none()) { + py::gil_scoped_acquire acquire; + py::object decomposed_py_dem = full_decomposer(dem); + stim_dem = + stim::DetectorErrorModel(py::cast(py::str(decomposed_py_dem)).c_str()); + } else { + stim_dem = stim::DetectorErrorModel(py::cast(py::str(dem)).c_str()); + } + + py::object python_classifier = detector_classifier; + tesseract::DetectorClassifier classifier = + [python_classifier](int index, const std::vector& coordinates, + const std::string& tag) -> int { + py::gil_scoped_acquire acquire; + return py::cast(python_classifier(index, coordinates, tag)); + }; + std::vector classification = + tesseract::MultiPassTesseractDecoder::classify_detectors(stim_dem, classifier); + + auto decoder = std::make_unique( + stim_dem, num_passes, classification, base_config, num_det_orders, det_order_method, seed, + strategy); + + return MultiPassSinterCompiledDecoder(std::move(decoder), stim_dem.count_detectors(), + stim_dem.count_observables()); + } +}; + +void pybind_multi_pass_sinter_compat(py::module& m) { + py::enum_(m, "SchedulingStrategy") + .value("Static", SchedulingStrategy::Static) + .value("Causal", SchedulingStrategy::Causal) + .export_values(); + + py::class_( + m, "MultiPassSinterCompiledDecoder", + "A compiled Sinter decoder backed by the native multi-pass Tesseract decoder.") + .def_property_readonly("num_components", &MultiPassSinterCompiledDecoder::num_components) + .def("decode_shots_bit_packed", &MultiPassSinterCompiledDecoder::decode_shots_bit_packed, + py::kw_only(), py::arg("bit_packed_detection_event_data"), + py::call_guard()); + + py::class_( + m, "MultiPassSinterDecoder", + "Low-level multi-pass Sinter decoder. Supports one or two passes and requires a " + "detector_classifier before compilation.") + .def(py::init(), py::arg("num_passes") = 2) + .def_readwrite("full_decomposer", &MultiPassSinterDecoder::full_decomposer) + .def_readwrite("detector_classifier", &MultiPassSinterDecoder::detector_classifier) + .def_readwrite("base_config", &MultiPassSinterDecoder::base_config) + .def_readwrite("num_det_orders", &MultiPassSinterDecoder::num_det_orders) + .def_readwrite("det_order_method", &MultiPassSinterDecoder::det_order_method) + .def_readwrite("seed", &MultiPassSinterDecoder::seed) + .def_readwrite("strategy", &MultiPassSinterDecoder::strategy) + .def("compile_decoder_for_dem", &MultiPassSinterDecoder::compile_decoder_for_dem, + py::kw_only(), py::arg("dem"), + "Compiles a DEM after classifying every detector into exactly two components."); +} + +} // namespace tesseract + +#endif // MULTI_PASS_SINTER_COMPAT_PYBIND_H diff --git a/src/multi_pass/multi_pass_tesseract_decoder.cc b/src/multi_pass/multi_pass_tesseract_decoder.cc new file mode 100644 index 00000000..7ddc23b5 --- /dev/null +++ b/src/multi_pass/multi_pass_tesseract_decoder.cc @@ -0,0 +1,529 @@ +#include "multi_pass_tesseract_decoder.h" + +#include +#include +#include +#include +#include +#include + +#include "../common.h" +#include "dem_decomposition.h" + +namespace tesseract { + +std::string MultiPassExecutionPlan::str() const { + std::stringstream ss; + ss << "Multi-pass execution plan\n" + << "strategy: " << (strategy == SchedulingStrategy::Static ? "static" : "causal") << '\n' + << "passes: " << num_passes << '\n' + << "monolithic input DEM: detectors=" << monolithic_statistics.detector_count + << ", error_mechanisms=" << monolithic_statistics.error_mechanism_count + << ", average_detector_degree=" << monolithic_statistics.average_detector_degree << '\n' + << "components: " << components.size() << '\n'; + for (const auto& component : components) { + ss << " component " << component.id << ": label=" << component.classifier_label + << ", active_detectors=" << component.active_detector_count + << ", decoder_detectors=" << component.decoder_detector_count + << ", observable=" << (component.affects_observable ? "yes" : "no") + << ", error_mechanisms=" << component.error_mechanism_count + << ", average_active_detector_degree=" << component.average_active_detector_degree << '\n'; + } + ss << "dependencies:\n"; + if (dependencies.empty()) ss << " none\n"; + for (const auto& dependency : dependencies) { + ss << " component " << dependency.source_component << " -> component " + << dependency.target_component << ": " << dependency.rule_count << " rules\n"; + } + ss << "schedule:\n"; + for (size_t pass = 0; pass < pass_schedule.size(); ++pass) { + ss << " pass " << pass + 1 << ": ["; + for (size_t i = 0; i < pass_schedule[pass].size(); ++i) { + if (i) ss << ", "; + ss << pass_schedule[pass][i]; + } + ss << "]\n"; + } + return ss.str(); +} + +namespace { + +struct DetectorMetadata { + std::map> coordinates; + std::map tags; +}; + +MultiPassExecutionPlan::DemStatistics dem_statistics(size_t detector_count, + const std::vector& errors) { + size_t detector_incidences = 0; + for (const auto& error : errors) { + detector_incidences += error.symptom.detectors.size(); + } + double average_detector_degree = + detector_count == 0 ? 0.0 : (double)detector_incidences / detector_count; + return {detector_count, errors.size(), average_detector_degree}; +} + +DetectorMetadata collect_detector_metadata(const stim::DetectorErrorModel& flattened) { + std::set detector_ids; + for (uint64_t d = 0; d < flattened.count_detectors(); ++d) { + detector_ids.insert(d); + } + + DetectorMetadata metadata; + metadata.coordinates = flattened.get_detector_coordinates(detector_ids); + for (const auto& instruction : flattened.instructions) { + if (instruction.type == stim::DemInstructionType::DEM_DETECTOR) { + metadata.tags[instruction.target_data[0].val()] = instruction.tag; + } + } + return metadata; +} + +void validate_detector_classes(const std::vector& detector_classes, size_t num_detectors) { + if (detector_classes.size() != num_detectors) { + throw std::invalid_argument("Detector classification count does not match the DEM."); + } + + std::set unique_classes; + for (size_t d = 0; d < detector_classes.size(); ++d) { + int classifier_label = detector_classes[d]; + if (classifier_label < 0) { + throw std::invalid_argument( + "Detector D" + std::to_string(d) + + " could not be classified (missing basis annotation or valid coordinates)."); + } + unique_classes.insert(classifier_label); + } + + if (unique_classes.size() != 2) { + throw std::invalid_argument("Multi-pass decoding requires exactly 2 detector components; got " + + std::to_string(unique_classes.size()) + "."); + } +} + +} // namespace + +MultiPassTesseractDecoder::MultiPassTesseractDecoder( + const stim::DetectorErrorModel& dem, size_t num_passes, const DetectorClassifier& classifier, + const TesseractConfig& base_config, size_t num_det_orders, DetOrder det_order_method, + uint64_t seed, SchedulingStrategy strategy, bool collect_plan_statistics) + : num_passes(num_passes), + strategy(strategy), + total_global_detectors(dem.count_detectors()), + base_config(base_config), + num_det_orders(num_det_orders), + det_order_method(det_order_method), + seed(seed), + collect_plan_statistics(collect_plan_statistics) { + if (num_passes < 1 || num_passes > 2) { + throw std::invalid_argument("num_passes must be 1 or 2."); + } + initialize(dem, classify_detectors(dem, classifier)); +} + +MultiPassTesseractDecoder::MultiPassTesseractDecoder( + const stim::DetectorErrorModel& dem, size_t num_passes, + const std::vector& detector_classes, const TesseractConfig& base_config, + size_t num_det_orders, DetOrder det_order_method, uint64_t seed, SchedulingStrategy strategy, + bool collect_plan_statistics) + : num_passes(num_passes), + strategy(strategy), + total_global_detectors(dem.count_detectors()), + base_config(base_config), + num_det_orders(num_det_orders), + det_order_method(det_order_method), + seed(seed), + collect_plan_statistics(collect_plan_statistics) { + if (num_passes < 1 || num_passes > 2) { + throw std::invalid_argument("num_passes must be 1 or 2."); + } + initialize(dem, detector_classes); +} + +std::vector MultiPassTesseractDecoder::classify_detectors( + const stim::DetectorErrorModel& dem, const DetectorClassifier& classifier) { + stim::DetectorErrorModel flattened = dem.flattened(); + DetectorMetadata metadata = collect_detector_metadata(flattened); + std::vector detector_classes(flattened.count_detectors()); + for (size_t d = 0; d < detector_classes.size(); ++d) { + const std::vector& coordinates = metadata.coordinates[d]; + const std::string& tag = metadata.tags[d]; + detector_classes[d] = classifier((int)d, coordinates, tag); + } + validate_detector_classes(detector_classes, flattened.count_detectors()); + return detector_classes; +} + +void MultiPassTesseractDecoder::initialize(const stim::DetectorErrorModel& dem, + const std::vector& detector_classes) { + stim::DetectorErrorModel flattened = dem.flattened(); + total_global_detectors = (size_t)flattened.count_detectors(); + validate_detector_classes(detector_classes, total_global_detectors); + DetectorMetadata metadata = collect_detector_metadata(flattened); + + if (collect_plan_statistics) { + std::vector error_index_map; + stim::DetectorErrorModel monolithic_dem = + common::merge_indistinguishable_errors(flattened, error_index_map); + monolithic_dem = common::remove_zero_probability_errors(monolithic_dem, error_index_map); + monolithic_statistics = + dem_statistics(monolithic_dem.count_detectors(), get_errors_from_dem(monolithic_dem)); + } + + std::set unique_classes; + unique_classes.insert(detector_classes.begin(), detector_classes.end()); + + std::map class_to_comp_id; + int next_comp_id = 0; + for (int c : unique_classes) class_to_comp_id[c] = next_comp_id++; + + component_decoders.resize(unique_classes.size()); + for (const auto& [classifier_label, component_id] : class_to_comp_id) { + component_decoders[component_id].classifier_label = classifier_label; + } + + global_det_to_comp_id.resize(total_global_detectors); + for (size_t i = 0; i < total_global_detectors; ++i) { + int component_id = class_to_comp_id.at(detector_classes[i]); + global_det_to_comp_id[i] = component_id; + component_decoders[component_id].component_detectors.insert((int)i); + } + + auto detector_component = [&](int detector) { + if (detector < 0 || (size_t)detector >= global_det_to_comp_id.size()) { + throw std::invalid_argument("Detector D" + std::to_string(detector) + " is out of range."); + } + return global_det_to_comp_id[detector]; + }; + + stim::DetectorErrorModel decomposed = + decompose_errors_using_detector_assignment(flattened, detector_component, true); + + ImpliedProbsMap raw_correlations = process_dem_correlations(decomposed, global_det_to_comp_id); + + auto component_dems = split_dem_by_component(decomposed, detector_component); + + for (size_t i = 0; i < component_decoders.size(); ++i) { + auto& cd = component_decoders[i]; + + std::vector error_index_map; + stim::DetectorErrorModel component_dem = + common::merge_indistinguishable_errors(component_dems[i], error_index_map); + component_dem = common::remove_zero_probability_errors(component_dem, error_index_map); + + for (size_t global_d = 0; global_d < total_global_detectors; ++global_d) { + cd.global_to_local_det[global_d] = (int)global_d; + } + + stim::DetectorErrorModel local_dem; + for (size_t global_d = 0; global_d < total_global_detectors; ++global_d) { + local_dem.append_detector_instruction(metadata.coordinates[global_d], + stim::DemTarget::relative_detector_id(global_d), + metadata.tags[global_d]); + } + + for (const auto& inst : component_dem.instructions) { + if (inst.type == stim::DemInstructionType::DEM_ERROR) { + bool has_obs = false; + for (const auto& t : inst.target_data) { + if (t.is_observable_id()) has_obs = true; + } + if (has_obs) cd.affects_observable = true; + local_dem.append_error_instruction(inst.arg_data[0], inst.target_data, inst.tag); + } else if (inst.type == stim::DemInstructionType::DEM_LOGICAL_OBSERVABLE) { + local_dem.append_dem_instruction(inst); + } + } + + TesseractConfig config = base_config; + config.dem = local_dem; + config.merge_errors = true; + config.det_orders = build_det_orders(config.dem, num_det_orders, det_order_method, seed); + + cd.decoder = std::make_unique(config); + if (base_config.verbose) { + std::cout << "DEBUG: Component " << i << " initialized with " << cd.decoder->errors.size() + << " errors and " << config.dem.count_detectors() << " detectors." << std::endl; + } + cd.error_index_to_rules.resize(cd.decoder->errors.size()); + + for (size_t ei = 0; ei < cd.decoder->errors.size(); ++ei) { + cd.original_costs.push_back(cd.decoder->errors[ei].likelihood_cost); + ComponentSymptom global_symptom{cd.decoder->errors[ei].symptom.detectors, + cd.decoder->errors[ei].symptom.observables}; + std::sort(global_symptom.detectors.begin(), global_symptom.detectors.end()); + std::sort(global_symptom.observables.begin(), global_symptom.observables.end()); + cd.symptom_to_error_index[global_symptom].push_back(ei); + } + } + + for (const auto& [global_symptom, implied_probs] : raw_correlations) { + int causal_comp = global_det_to_comp_id[global_symptom.detectors[0]]; + + auto it = component_decoders[causal_comp].symptom_to_error_index.find(global_symptom); + if (it == component_decoders[causal_comp].symptom_to_error_index.end()) continue; + + // Loop through all degenerate causal error indices! + for (size_t causal_err_idx : it->second) { + for (const auto& imp : implied_probs) { + const ComponentSymptom& target_symptom = imp.affected_symptom; + int target_comp = global_det_to_comp_id[target_symptom.detectors[0]]; + + auto t_it = component_decoders[target_comp].symptom_to_error_index.find(target_symptom); + if (t_it != component_decoders[target_comp].symptom_to_error_index.end()) { + // Loop through all degenerate target error indices and add rules to + // each! + for (size_t target_err_idx : t_it->second) { + component_decoders[causal_comp].error_index_to_rules[causal_err_idx].push_back( + {(size_t)target_comp, target_err_idx, imp.probability}); + } + } + } + } + } + + if (strategy == SchedulingStrategy::Static) { + build_static_schedule(); + } else if (strategy == SchedulingStrategy::Causal) { + build_causal_schedule(); + } +} + +void MultiPassTesseractDecoder::build_static_schedule() { + pass_schedule.assign(num_passes, {}); + for (size_t p = 0; p < num_passes; ++p) { + for (size_t i = 0; i < component_decoders.size(); ++i) { + pass_schedule[p].push_back(i); + } + } +} + +void MultiPassTesseractDecoder::build_causal_schedule() { + size_t num_components = component_decoders.size(); + std::vector> schedule_sets(num_passes); + + // Initial seed: Final pass includes all components that directly affect an + // observable. + for (size_t i = 0; i < num_components; ++i) { + if (component_decoders[i].affects_observable) { + schedule_sets[num_passes - 1].insert(i); + } + } + + // Back-propagate dependencies through passes. + // A component is needed in pass p if it can reweight a component needed in + // pass p+1. + for (int p = (int)num_passes - 2; p >= 0; --p) { + // Start with everyone needed in the next pass (they might need to re-decode + // or bias others) Actually, if a component is in pass p+1, it's because it + // was influenced by pass p. + for (size_t target_comp_idx : schedule_sets[p + 1]) { + for (size_t causal_comp_idx = 0; causal_comp_idx < num_components; ++causal_comp_idx) { + for (const auto& rules : component_decoders[causal_comp_idx].error_index_to_rules) { + for (const auto& rule : rules) { + if (rule.target_comp_idx == target_comp_idx) { + schedule_sets[p].insert(causal_comp_idx); + } + } + } + } + } + } + + // Convert sets to pass_schedule vectors. + pass_schedule.assign(num_passes, {}); + for (size_t p = 0; p < num_passes; ++p) { + for (size_t c_idx : schedule_sets[p]) { + pass_schedule[p].push_back(c_idx); + } + } +} + +MultiPassExecutionPlan MultiPassTesseractDecoder::get_execution_plan() const { + if (!collect_plan_statistics) { + throw std::logic_error("Execution plan statistics were not collected."); + } + MultiPassExecutionPlan plan{num_passes, strategy, monolithic_statistics, {}, {}, pass_schedule}; + for (size_t component_id = 0; component_id < component_decoders.size(); ++component_id) { + const auto& component = component_decoders[component_id]; + auto statistics = + dem_statistics(component.component_detectors.size(), component.decoder->errors); + plan.components.push_back({component_id, component.classifier_label, statistics.detector_count, + component.decoder->num_detectors, statistics.error_mechanism_count, + statistics.average_detector_degree, component.affects_observable}); + } + + std::map, size_t> dependency_counts; + for (size_t source = 0; source < component_decoders.size(); ++source) { + for (const auto& rules : component_decoders[source].error_index_to_rules) { + for (const auto& rule : rules) { + dependency_counts[{source, rule.target_comp_idx}]++; + } + } + } + for (const auto& [components, rule_count] : dependency_counts) { + plan.dependencies.push_back({components.first, components.second, rule_count}); + } + return plan; +} + +std::vector MultiPassTesseractDecoder::decode(const std::vector& detections) { + return decode_result(detections).predictions; +} + +MultiPassDecodeResult MultiPassTesseractDecoder::decode_result( + const std::vector& detections) { + for (uint64_t d : detections) { + if (d >= total_global_detectors) { + throw std::invalid_argument("Detector D" + std::to_string(d) + + " is out of range for a model with " + + std::to_string(total_global_detectors) + " detectors."); + } + } + + last_shot_num_reweights = 0; + + // 1. Multi-Pass Loop: Sequentially schedules component passes and propagates + // priors. + for (size_t pass = 0; pass < num_passes; ++pass) { + bool is_final_pass = (pass == num_passes - 1); + + // Decode scheduled components for the current pass layer using persistent + // local buffers. + for (size_t comp_idx : pass_schedule[pass]) { + auto& cd = component_decoders[comp_idx]; + std::vector local_dets; + for (uint64_t d : detections) { + if (cd.component_detectors.count(static_cast(d))) { + local_dets.push_back(d); + } + } + + cd.decoder->decode_to_errors(local_dets); + component_predictions[comp_idx] = cd.decoder->predicted_errors_buffer; + } + + if (!is_final_pass) { + // Step A: Apply Damped Fractional Memory to previously modified priors. + // Smoothly decay current modifications back toward the baseline to + // prevent message saturation. + double gamma = 0.5; // Tunable decay factor: 1.0 is strict isolation, 0.0 + // is full accumulation. + + for (size_t m_comp_idx : modified_component_indices) { + auto& cd = component_decoders[m_comp_idx]; + if (!cd.shot_all_modified_error_indices.empty()) { + for (size_t idx : cd.shot_all_modified_error_indices) { + double baseline_cost = cd.original_costs[idx]; + double current_cost = cd.decoder->errors[idx].likelihood_cost; + cd.decoder->errors[idx].likelihood_cost = + gamma * baseline_cost + (1.0 - gamma) * current_cost; + } + cd.decoder->update_internal_costs(cd.shot_all_modified_error_indices); + // Retain tracking indices so the final Surgical Reset completely + // clears cross-shot state. + } + } + + // Step B: Broadcast reweighting rules derived strictly from the latest + // predictions. + for (size_t comp_idx : pass_schedule[pass]) { + auto& cd = component_decoders[comp_idx]; + for (size_t dem_err_idx : cd.decoder->predicted_errors_buffer) { + size_t internal_err_idx = cd.decoder->dem_error_to_error.at(dem_err_idx); + if (internal_err_idx == std::numeric_limits::max()) continue; + + for (const auto& rule : cd.error_index_to_rules[internal_err_idx]) { + auto& target_cd = component_decoders[rule.target_comp_idx]; + + modified_component_indices.push_back(rule.target_comp_idx); + + // Apply Max-Prob Rule safely for concurrent rules within this pass + // layer. + double current_p = target_cd.decoder->errors[rule.target_error_idx].get_probability(); + if (rule.conditional_prob > current_p) { + target_cd.decoder->errors[rule.target_error_idx].set_with_probability( + std::min(rule.conditional_prob, 0.5)); + target_cd.shot_all_modified_error_indices.push_back(rule.target_error_idx); + last_shot_num_reweights++; + } + } + } + } + + // Step C: Deduplicate modified tracking vectors and synchronize internal + // graph costs. + std::sort(modified_component_indices.begin(), modified_component_indices.end()); + modified_component_indices.erase( + std::unique(modified_component_indices.begin(), modified_component_indices.end()), + modified_component_indices.end()); + + for (size_t m_comp_idx : modified_component_indices) { + auto& cd = component_decoders[m_comp_idx]; + if (!cd.shot_all_modified_error_indices.empty()) { + std::sort(cd.shot_all_modified_error_indices.begin(), + cd.shot_all_modified_error_indices.end()); + cd.shot_all_modified_error_indices.erase( + std::unique(cd.shot_all_modified_error_indices.begin(), + cd.shot_all_modified_error_indices.end()), + cd.shot_all_modified_error_indices.end()); + cd.decoder->update_internal_costs(cd.shot_all_modified_error_indices); + } + } + } + } + + // 2. Unified Logical Extraction: Collect final predictions from ALL + // components that ran during the shot. + std::set flipped_observables; + bool aggregate_low_confidence = false; + double aggregate_cost = 0.0; + + for (const auto& [comp_idx, preds] : component_predictions) { + auto& cd = component_decoders[comp_idx]; + if (cd.decoder->low_confidence_flag) { + aggregate_low_confidence = true; + } + if (!preds.empty()) { + std::vector local_flips = cd.decoder->get_flipped_observables(preds); + for (int obs : local_flips) { + if (flipped_observables.count(obs)) + flipped_observables.erase(obs); + else + flipped_observables.insert(obs); + } + } + } + + for (size_t comp_idx : pass_schedule.back()) { + auto& cd = component_decoders[comp_idx]; + const auto& preds = component_predictions.at(comp_idx); + aggregate_cost += cd.decoder->cost_from_errors(preds); + } + + // 3. Surgical Reset: Restore modified costs to leave the internal structures + // pristine for the next shot. + for (size_t m_comp_idx : modified_component_indices) { + auto& cd = component_decoders[m_comp_idx]; + if (!cd.shot_all_modified_error_indices.empty()) { + for (size_t idx : cd.shot_all_modified_error_indices) { + cd.decoder->errors[idx].likelihood_cost = cd.original_costs[idx]; + } + cd.decoder->update_internal_costs(cd.shot_all_modified_error_indices); + cd.shot_all_modified_error_indices.clear(); + } + } + + modified_component_indices.clear(); + final_pass_active_components.clear(); + + MultiPassDecodeResult res; + res.predictions = std::vector(flipped_observables.begin(), flipped_observables.end()); + res.low_confidence = aggregate_low_confidence; + res.total_cost = aggregate_cost; + return res; +} + +} // namespace tesseract diff --git a/src/multi_pass/multi_pass_tesseract_decoder.h b/src/multi_pass/multi_pass_tesseract_decoder.h new file mode 100644 index 00000000..709c2cb4 --- /dev/null +++ b/src/multi_pass/multi_pass_tesseract_decoder.h @@ -0,0 +1,170 @@ +#ifndef MULTI_PASS_TESSERACT_DECODER_H +#define MULTI_PASS_TESSERACT_DECODER_H + +#include +#include +#include +#include + +#include "../tesseract.h" +#include "../utils.h" +#include "dem_decomposition.h" +#include "error_correlations.h" +#include "stim.h" +#include "tanner_graph.h" + +namespace tesseract { + +enum class SchedulingStrategy { + Static, // Schedules both components in every pass. + Causal // Derives each pass from component dependencies. +}; + +struct MultiPassExecutionPlan { + struct DemStatistics { + size_t detector_count; + size_t error_mechanism_count; + double average_detector_degree; + }; + + struct Component { + size_t id; + int classifier_label; + size_t active_detector_count; + size_t decoder_detector_count; + size_t error_mechanism_count; + double average_active_detector_degree; + bool affects_observable; + }; + + struct Dependency { + size_t source_component; + size_t target_component; + size_t rule_count; + }; + + size_t num_passes; + SchedulingStrategy strategy; + DemStatistics monolithic_statistics; + std::vector components; + std::vector dependencies; + std::vector> pass_schedule; + + std::string str() const; +}; + +struct MultiPassDecodeResult { + std::vector predictions; + bool low_confidence = false; + double total_cost = 0.0; // Cost of predictions made during the final pass. +}; + +/** + * Decodes a detector error model by splitting it into exactly two detector components. + * + * One or two passes are supported. Every detector must receive a nonnegative classifier label, + * and exactly two distinct labels must be present. + */ +class MultiPassTesseractDecoder { + public: + MultiPassTesseractDecoder(const stim::DetectorErrorModel& dem, size_t num_passes, + const DetectorClassifier& classifier, + const TesseractConfig& base_config = TesseractConfig(), + size_t num_det_orders = 1, + DetOrder det_order_method = DetOrder::DetIndex, uint64_t seed = 0, + SchedulingStrategy strategy = SchedulingStrategy::Causal, + bool collect_plan_statistics = false); + MultiPassTesseractDecoder(const stim::DetectorErrorModel& dem, size_t num_passes, + const std::vector& detector_classes, + const TesseractConfig& base_config = TesseractConfig(), + size_t num_det_orders = 1, + DetOrder det_order_method = DetOrder::DetIndex, uint64_t seed = 0, + SchedulingStrategy strategy = SchedulingStrategy::Causal, + bool collect_plan_statistics = false); + + static std::vector classify_detectors(const stim::DetectorErrorModel& dem, + const DetectorClassifier& classifier); + + /** Returns the component schedule and statistics; requires collect_plan_statistics=true. */ + MultiPassExecutionPlan get_execution_plan() const; + std::vector decode(const std::vector& detections); + MultiPassDecodeResult decode_result(const std::vector& detections); + + size_t get_last_shot_num_reweights() const { + return last_shot_num_reweights; + } + size_t num_components() const { + return component_decoders.size(); + } + + private: + struct LocalReweightRule { + size_t target_comp_idx; + size_t target_error_idx; + double conditional_prob; + }; + + struct ComponentDecoder { + std::unique_ptr decoder; + int classifier_label = -1; + std::set component_detectors; // Global indices + std::map global_to_local_det; + std::vector original_costs; + std::map> symptom_to_error_index; + std::vector> error_index_to_rules; + std::vector modified_error_indices; + std::vector shot_all_modified_error_indices; + bool affects_observable = false; + }; + + size_t num_passes; + SchedulingStrategy strategy; + size_t total_global_detectors; + TesseractConfig base_config; + size_t num_det_orders; + ::DetOrder det_order_method; + uint64_t seed; + bool collect_plan_statistics; + size_t last_shot_num_reweights = 0; + MultiPassExecutionPlan::DemStatistics monolithic_statistics{}; + std::map> component_predictions; + std::vector modified_component_indices; + std::vector final_pass_active_components; + std::vector component_decoders; + std::vector> pass_schedule; + std::vector global_det_to_comp_id; + + void initialize(const stim::DetectorErrorModel& dem, const std::vector& detector_classes); + void build_static_schedule(); + void build_causal_schedule(); + + friend class MultiPassTraceVisualizer; + friend class MultiPassDebugger; +}; + +class MultiPassDebugger { + public: + static const std::vector>& get_pass_schedule( + const MultiPassTesseractDecoder& decoder) { + return decoder.pass_schedule; + } + static size_t num_components(const MultiPassTesseractDecoder& decoder) { + return decoder.component_decoders.size(); + } + static const TesseractDecoder& get_component_decoder(const MultiPassTesseractDecoder& decoder, + size_t i) { + return *decoder.component_decoders[i].decoder; + } + static const std::vector& get_modified_component_indices( + const MultiPassTesseractDecoder& decoder) { + return decoder.modified_component_indices; + } + static const MultiPassTesseractDecoder::ComponentDecoder& get_component_decoder_full( + const MultiPassTesseractDecoder& decoder, size_t i) { + return decoder.component_decoders[i]; + } +}; + +} // namespace tesseract + +#endif // MULTI_PASS_TESSERACT_DECODER_H diff --git a/src/multi_pass/multi_pass_tesseract_decoder.test.cc b/src/multi_pass/multi_pass_tesseract_decoder.test.cc new file mode 100644 index 00000000..0a942f19 --- /dev/null +++ b/src/multi_pass/multi_pass_tesseract_decoder.test.cc @@ -0,0 +1,471 @@ +#include "multi_pass_tesseract_decoder.h" + +#include +#include +#include +#include + +#include "gtest/gtest.h" + +using namespace tesseract; + +stim::DetectorErrorModel load_test_dem(const std::string& filename) { + std::string path = "testdata/surfacecodes/" + filename; + std::ifstream is(path); + if (!is.is_open()) { + is.open(filename); + } + if (!is.is_open()) { + throw std::runtime_error("Could not open file: " + filename); + } + std::stringstream ss; + ss << is.rdbuf(); + stim::Circuit circuit(ss.str().c_str()); + return stim::ErrorAnalyzer::circuit_to_detector_error_model(circuit, true, true, false, false, + false, 0.0); +} + +auto chromobius_classifier = [](int index, const std::vector& coords, + const std::string& tag) -> int { + if (coords.size() < 4) return -1; + int c3 = (int)coords[3]; + if (c3 >= 0 && c3 <= 2) return 0; // Basis X + if (c3 >= 3 && c3 <= 5) return 1; // Basis Z + return -1; +}; + +TEST(MultiPassTesseractDecoderTest, AcceptsOnlyOneOrTwoPasses) { + stim::DetectorErrorModel dem(R"DEM( + error(0.1) D0 + error(0.1) D1 L0 + detector D0 + detector D1 + logical_observable L0 + )DEM"); + auto classifier = [](int index, const std::vector&, const std::string&) -> int { + return index; + }; + + for (auto strategy : {SchedulingStrategy::Static, SchedulingStrategy::Causal}) { + for (size_t num_passes : {1, 2}) { + EXPECT_NO_THROW(MultiPassTesseractDecoder(dem, num_passes, classifier, TesseractConfig(), 1, + DetOrder::DetBFS, 0, strategy)); + } + for (size_t num_passes : {0, 3, 4}) { + EXPECT_THROW(MultiPassTesseractDecoder(dem, num_passes, classifier, TesseractConfig(), 1, + DetOrder::DetBFS, 0, strategy), + std::invalid_argument); + } + } +} + +TEST(MultiPassTesseractDecoderTest, RequiresTwoFullyClassifiedComponents) { + stim::DetectorErrorModel dem(R"DEM( + error(0.1) D0 + error(0.1) D1 L0 + error(0.1) D2 + detector D0 + detector D1 + detector D2 + logical_observable L0 + )DEM"); + + std::vector calls(3); + std::vector labels = {4, 4, 9}; + auto classifier = [&](int index, const std::vector&, const std::string&) -> int { + calls[index]++; + return labels[index]; + }; + MultiPassTesseractDecoder decoder(dem, 1, classifier); + EXPECT_EQ(decoder.num_components(), 2); + EXPECT_EQ(calls, std::vector({1, 1, 1})); + + EXPECT_THROW(MultiPassTesseractDecoder(dem, 1, std::vector({4, 4, 4})), + std::invalid_argument); + EXPECT_THROW(MultiPassTesseractDecoder(dem, 1, std::vector({4, 9, 12})), + std::invalid_argument); + + try { + MultiPassTesseractDecoder(dem, 1, std::vector({4, -1, 9})); + FAIL() << "Expected an unclassified detector to be rejected."; + } catch (const std::invalid_argument& error) { + EXPECT_NE(std::string(error.what()).find("D1"), std::string::npos); + } +} + +TEST(MultiPassTesseractDecoderTest, TwoPassCorrelationBenefit) { + // Component 0: D0 (Causal) + // Component 1: D1 (Affected) -> Observable L0 + // Rule: D0 ^ D1 exists with probability 0.1 + // Independent: D0 with prob 0.01, D1 with prob 0.2 + // If D0 is detected and explained by the bridging error, D1's probability + // should increase. + + stim::DetectorErrorModel dem(R"DEM( + error(0.1) D0 ^ D1 L0 + error(0.01) D0 + error(0.2) D1 L0 + detector D0 + detector D1 + logical_observable L0 + )DEM"); + + // Classifier: D0 -> Comp 0, D1 -> Comp 1 + auto classifier = [](int index, const std::vector& coords, + const std::string& tag) -> int { return index; }; + + TesseractConfig config; + config.verbose = true; + MultiPassTesseractDecoder decoder(dem, 2, classifier, config); + + // Shot 1: D0 and D1 both fire. + // Pass 1: Decode Comp 0. D0 is explained by the bridging error (implicit). + // Reweight: D1 L0 in Comp 1 becomes more likely. + // Pass 2: Decode Comp 1. + std::vector detections = {0, 1}; + MultiPassDecodeResult result = decoder.decode_result(detections); + + // In this specific model, if D0 and D1 both fire, + // the most likely explanation is the bridging error (0.1) + // vs independent (0.01 * 0.2 = 0.002). + // The bridging error flips L0. + // So we expect L0 to be flipped. + ASSERT_TRUE(std::find(result.predictions.begin(), result.predictions.end(), 0) != + result.predictions.end()); + EXPECT_DOUBLE_EQ(result.total_cost, 0.0); +} + +TEST(MultiPassTesseractDecoderTest, DisjointDecoding) { + stim::DetectorErrorModel dem(R"DEM( + error(0.1) D0 L0 + error(0.1) D1 L1 + detector D0 + detector D1 + logical_observable L0 + logical_observable L1 + )DEM"); + + auto classifier = [](int index, const std::vector& coords, + const std::string& tag) -> int { return index; }; + + MultiPassTesseractDecoder decoder(dem, 1, classifier); + + std::vector detections = {0}; + std::vector result = decoder.decode(detections); + ASSERT_EQ(result.size(), 1); + ASSERT_EQ(result[0], 0); + + detections = {1}; + result = decoder.decode(detections); + ASSERT_EQ(result.size(), 1); + ASSERT_EQ(result[0], 1); +} + +TEST(MultiPassTesseractDecoderTest, CausalScheduleSurfaceCode) { + // A simplified d=2 surface code style DEM + // D0, D1: Basis X (Class 0), Affected by correlations from Basis Z + // D2, D3: Basis Z (Class 1), Causal (Reweight Basis X) + // Error: D2 ^ D0 (Bridge) + stim::DetectorErrorModel dem(R"DEM( + error(0.1) D0 D2 L0 + error(0.01) D0 L0 + error(0.01) D2 + error(0.1) D1 D3 L0 + error(0.01) D1 L0 + error(0.01) D3 + detector D0 + detector D1 + detector D2 + detector D3 + logical_observable L0 + )DEM"); + + // Class 0: Detectors 0, 1 + // Class 1: Detectors 2, 3 + auto classifier = [](int index, const std::vector& coords, + const std::string& tag) -> int { return (index < 2) ? 0 : 1; }; + + MultiPassTesseractDecoder decoder(dem, 2, classifier, TesseractConfig(), 1, DetOrder::DetBFS, 0, + SchedulingStrategy::Causal); + + const auto& schedule = MultiPassDebugger::get_pass_schedule(decoder); + ASSERT_EQ(schedule.size(), 2); + + ASSERT_EQ(schedule[0].size(), 1); + ASSERT_EQ(schedule[0][0], 1); // Component 1 (Class 1) runs first + ASSERT_EQ(schedule[1].size(), 1); + ASSERT_EQ(schedule[1][0], 0); // Component 0 (Class 0) runs last +} + +TEST(MultiPassTesseractDecoderTest, ExecutionPlanReflectsDecoderState) { + stim::DetectorErrorModel dem(R"DEM( + error(0.1) D0 D1 L0 + error(0.01) D0 + error(0.2) D1 L0 + detector D0 + detector D1 + logical_observable L0 + )DEM"); + + MultiPassTesseractDecoder decoder(dem, 2, std::vector({4, 9}), TesseractConfig(), 1, + DetOrder::DetIndex, 0, SchedulingStrategy::Causal, true); + MultiPassExecutionPlan plan = decoder.get_execution_plan(); + + EXPECT_EQ(plan.strategy, SchedulingStrategy::Causal); + EXPECT_EQ(plan.monolithic_statistics.detector_count, 2); + EXPECT_EQ(plan.monolithic_statistics.error_mechanism_count, 3); + EXPECT_DOUBLE_EQ(plan.monolithic_statistics.average_detector_degree, 2.0); + ASSERT_EQ(plan.components.size(), 2); + EXPECT_EQ(plan.components[0].classifier_label, 4); + EXPECT_EQ(plan.components[1].classifier_label, 9); + for (const auto& component : plan.components) { + EXPECT_EQ(component.active_detector_count, 1); + EXPECT_EQ(component.decoder_detector_count, 2); + EXPECT_EQ(component.error_mechanism_count, 1); + EXPECT_DOUBLE_EQ(component.average_active_detector_degree, 1.0); + } + ASSERT_EQ(plan.dependencies.size(), 2); + EXPECT_EQ(plan.dependencies[0].source_component, 0); + EXPECT_EQ(plan.dependencies[0].target_component, 1); + EXPECT_GT(plan.dependencies[0].rule_count, 0); + EXPECT_EQ(plan.pass_schedule, std::vector>({{0}, {1}})); + EXPECT_NE(plan.str().find("monolithic input DEM: detectors=2, error_mechanisms=3"), + std::string::npos); + EXPECT_NE(plan.str().find("component 0: label=4, active_detectors=1, decoder_detectors=2"), + std::string::npos); + EXPECT_NE(plan.str().find("pass 2: [1]"), std::string::npos); +} + +TEST(MultiPassTesseractDecoderTest, SurfaceCodePartitioning) { + std::vector distances = {3, 5, 7}; + for (int d : distances) { + int q = 2 * d * d - 1; + std::string filename = "r=" + std::to_string(d) + ",d=" + std::to_string(d) + + ",p=0.001,noise=si1000,c=surface_code_X,q=" + std::to_string(q) + + ",gates=cz.stim"; + stim::DetectorErrorModel dem = load_test_dem(filename); + MultiPassTesseractDecoder decoder(dem, 1, chromobius_classifier); + ASSERT_EQ(decoder.num_components(), 2) << "Failed partitioning for d=" << d; + } +} + +TEST(MultiPassTesseractDecoderTest, SurfaceCodeCausalScheduling) { + std::vector distances = {3, 5, 7}; + for (int d : distances) { + int q = 2 * d * d - 1; + std::string filename = "r=" + std::to_string(d) + ",d=" + std::to_string(d) + + ",p=0.001,noise=si1000,c=surface_code_X,q=" + std::to_string(q) + + ",gates=cz.stim"; + stim::DetectorErrorModel dem = load_test_dem(filename); + + // 1-Pass: Should only schedule X component (0) + { + MultiPassTesseractDecoder decoder(dem, 1, chromobius_classifier, TesseractConfig(), 1, + DetOrder::DetBFS, 0, SchedulingStrategy::Causal); + const auto& schedule = MultiPassDebugger::get_pass_schedule(decoder); + ASSERT_EQ(schedule.size(), 1); + ASSERT_EQ(schedule[0].size(), 1); + ASSERT_EQ(schedule[0][0], 0) << "1-pass failed for d=" << d; + } + + // 2-Pass: Should schedule Z (1) then X (0) + { + MultiPassTesseractDecoder decoder(dem, 2, chromobius_classifier, TesseractConfig(), 1, + DetOrder::DetBFS, 0, SchedulingStrategy::Causal); + const auto& schedule = MultiPassDebugger::get_pass_schedule(decoder); + ASSERT_EQ(schedule.size(), 2); + ASSERT_EQ(schedule[0].size(), 1); + ASSERT_EQ(schedule[0][0], 1) << "2-pass P0 failed for d=" << d; + ASSERT_EQ(schedule[1].size(), 1); + ASSERT_EQ(schedule[1][0], 0) << "2-pass P1 failed for d=" << d; + } + } +} + +TEST(MultiPassTesseractDecoderTest, PerfectResetSurfaceCode) { + std::vector distances = {3, 5, 7}; + for (int d : distances) { + int q = 2 * d * d - 1; + std::string filename = "r=" + std::to_string(d) + ",d=" + std::to_string(d) + + ",p=0.001,noise=si1000,c=surface_code_X,q=" + std::to_string(q) + + ",gates=cz.stim"; + stim::DetectorErrorModel dem = load_test_dem(filename); + MultiPassTesseractDecoder decoder(dem, 2, chromobius_classifier, TesseractConfig(), 1, + DetOrder::DetBFS, 0, SchedulingStrategy::Causal); + + size_t n_comp = MultiPassDebugger::num_components(decoder); + + // Capture initial state + std::vector> initial_likelihoods(n_comp); + std::vector> initial_error_costs(n_comp); + for (size_t i = 0; i < n_comp; ++i) { + const auto& comp_dec = MultiPassDebugger::get_component_decoder(decoder, i); + for (const auto& err : comp_dec.errors) { + initial_likelihoods[i].push_back(err.likelihood_cost); + } + initial_error_costs[i] = TesseractDebugger::get_error_costs(comp_dec); + } + + // Run shots + std::mt19937_64 rng(12345); + size_t total_reweights_in_test = 0; + for (int shot = 0; shot < 100; ++shot) { + std::vector detections; + for (uint64_t det_idx = 0; det_idx < dem.count_detectors(); ++det_idx) { + if (std::uniform_real_distribution(0, 1)(rng) < 0.05) { + detections.push_back(det_idx); + } + } + + decoder.decode(detections); + total_reweights_in_test += decoder.get_last_shot_num_reweights(); + + // Verify state is restored + for (size_t i = 0; i < n_comp; ++i) { + const auto& comp_dec = MultiPassDebugger::get_component_decoder(decoder, i); + + for (size_t ei = 0; ei < comp_dec.errors.size(); ++ei) { + ASSERT_DOUBLE_EQ(comp_dec.errors[ei].likelihood_cost, initial_likelihoods[i][ei]) + << "Likelihood mismatch at d=" << d << " shot=" << shot << " comp=" << i + << " err=" << ei; + } + + const auto& current_error_costs = TesseractDebugger::get_error_costs(comp_dec); + ASSERT_EQ(current_error_costs.size(), initial_error_costs[i].size()); + for (size_t ei = 0; ei < current_error_costs.size(); ++ei) { + ASSERT_DOUBLE_EQ(current_error_costs[ei].likelihood_cost, + initial_error_costs[i][ei].likelihood_cost) + << "Internal likelihood mismatch at d=" << d << " shot=" << shot << " comp=" << i + << " err=" << ei; + ASSERT_DOUBLE_EQ(current_error_costs[ei].min_cost, initial_error_costs[i][ei].min_cost) + << "Internal min_cost mismatch at d=" << d << " shot=" << shot << " comp=" << i + << " err=" << ei; + } + } + } + ASSERT_GT(total_reweights_in_test, 0) + << "Test was trivial for d=" << d << ". No reweighting occurred."; + } +} + +TEST(MultiPassTesseractDecoderTest, BoundaryConditionAndCappingTest) { + stim::DetectorErrorModel dem(R"DEM( + error(0.49) D0 D1 L0 + error(0.5) D0 + detector D0 + detector D1 + logical_observable L0 + )DEM"); + + auto classifier = [](int index, const std::vector& coords, + const std::string& tag) -> int { return index; }; + + TesseractConfig config; + config.dem = dem; + + MultiPassTesseractDecoder decoder(dem, 2, classifier, config, 1, DetOrder::DetIndex, 12345, + SchedulingStrategy::Causal); + + std::vector hits = {0}; + ASSERT_NO_THROW(decoder.decode(hits)); +} + +TEST(MultiPassTesseractDecoderTest, PriorReweightingOccurs) { + stim::DetectorErrorModel dem(R"DEM( + error(0.1) D0 D1 L0 + error(0.01) D0 + error(0.2) D1 L0 + detector D0 + detector D1 + logical_observable L0 + )DEM"); + + auto classifier = [](int index, const std::vector& coords, + const std::string& tag) -> int { return index; }; + + TesseractConfig config; + config.dem = dem; + + MultiPassTesseractDecoder decoder(dem, 2, classifier, config, 1, DetOrder::DetIndex, 12345, + SchedulingStrategy::Causal); + + std::vector hits = {0}; + decoder.decode(hits); + + // Rigorously assert that prior LLR reweights occurred successfully on the raw + // un-decomposed DEM! + ASSERT_GT(decoder.get_last_shot_num_reweights(), 0); +} + +TEST(MultiPassTesseractDecoderTest, MultipleCausalTriggersMaxProbValidation) { + stim::DetectorErrorModel dem(R"DEM( + error(0.1) D1 D0 L0 + error(0.15) D2 D0 L0 + error(0.01) D1 + error(0.01) D2 + error(0.2) D0 L0 + detector D0 + detector D1 + detector D2 + logical_observable L0 + )DEM"); + + auto classifier = [](int index, const std::vector& coords, + const std::string& tag) -> int { + if (index == 0) return 0; + return 1; + }; + + TesseractConfig config; + config.dem = dem; + + MultiPassTesseractDecoder decoder(dem, 2, classifier, config, 1, DetOrder::DetIndex, 12345, + SchedulingStrategy::Causal); + + // 1. Run a shot that triggers BOTH causal detectors D1 and D2 + std::vector detections = {1, 2}; + decoder.decode(detections); + + const auto& comp0 = MultiPassDebugger::get_component_decoder_full(decoder, 0); + + // 2. Find the target error index for symptom D0 L0. + ComponentSymptom target_symptom{{0}, {0}}; + auto it = comp0.symptom_to_error_index.find(target_symptom); + ASSERT_NE(it, comp0.symptom_to_error_index.end()); + + // Now fully active using vector degeneracy mapping! + size_t target_err_idx = it->second[0]; + + double final_p = comp0.decoder->errors[target_err_idx].get_probability(); + + // 3. Assert that the Surgical Reset successfully restored the cost back to + // the baseline after exactly 2 LLR reweighting rules were triggered and + // applied! + ASSERT_DOUBLE_EQ(final_p, 0.33199999999999996); + ASSERT_EQ(decoder.get_last_shot_num_reweights(), 2); +} + +TEST(MultiPassTesseractDecoderTest, CorrelationRulesDistinguishObservables) { + stim::DetectorErrorModel dem(R"DEM( + error(0.1) D0 D1 L0 + error(0.01) D0 + error(0.2) D1 L0 + error(0.05) D1 L1 + detector D0 + detector D1 + logical_observable L0 + logical_observable L1 + )DEM"); + + auto classifier = [](int index, const std::vector& coords, + const std::string& tag) -> int { return index; }; + + TesseractConfig config; + config.dem = dem; + + MultiPassTesseractDecoder decoder(dem, 2, classifier, config, 1, DetOrder::DetIndex, 12345, + SchedulingStrategy::Causal); + + decoder.decode({0}); + EXPECT_EQ(decoder.get_last_shot_num_reweights(), 1); +} diff --git a/src/multi_pass/tanner_graph.cc b/src/multi_pass/tanner_graph.cc new file mode 100644 index 00000000..327674a3 --- /dev/null +++ b/src/multi_pass/tanner_graph.cc @@ -0,0 +1,99 @@ +#include "tanner_graph.h" + +#include +#include + +namespace tesseract { + +std::vector TannerGraph::find_components(const stim::DetectorErrorModel& dem) { + int num_detectors = (int)dem.count_detectors(); + int num_observables = (int)dem.count_observables(); + int total_symptoms = num_detectors + num_observables; + + UnionFind uf(total_symptoms); + std::vector symptom_active(total_symptoms, false); + + // 1. Union symptoms connected by errors + auto flattened = dem.flattened(); + for (size_t i = 0; i < flattened.instructions.size(); ++i) { + const auto& inst = flattened.instructions[i]; + if (inst.type != stim::DemInstructionType::DEM_ERROR) continue; + + // Manually split by separators to handle decomposed errors + size_t group_start = 0; + for (size_t k = 0; k <= inst.target_data.size(); ++k) { + if (k == inst.target_data.size() || inst.target_data[k].is_separator()) { + std::vector group_symptoms; + for (size_t j = group_start; j < k; ++j) { + const auto& target = inst.target_data[j]; + int sym_id = -1; + if (target.is_relative_detector_id()) { + sym_id = target.val(); + } else if (target.is_observable_id()) { + sym_id = num_detectors + target.val(); + } + + if (sym_id != -1) { + group_symptoms.push_back(sym_id); + symptom_active[sym_id] = true; + } + } + + for (size_t j = 1; j < group_symptoms.size(); ++j) { + uf.unite(group_symptoms[0], group_symptoms[j]); + } + group_start = k + 1; + } + } + } + + // 2. Group symptoms by root + std::unordered_map root_to_component; + for (int i = 0; i < total_symptoms; ++i) { + if (!symptom_active[i]) continue; + + int root = uf.find(i); + if (root_to_component.find(root) == root_to_component.end()) { + root_to_component[root] = TannerComponent(); + } + + if (i < num_detectors) { + root_to_component[root].detectors.push_back(i); + } else { + root_to_component[root].observables.push_back(i - num_detectors); + root_to_component[root].affects_observable = true; + } + } + + // 3. Assign errors to components + for (size_t i = 0; i < flattened.instructions.size(); ++i) { + const auto& inst = flattened.instructions[i]; + if (inst.type != stim::DemInstructionType::DEM_ERROR) continue; + + std::set roots_touched; + for (const auto& target : inst.target_data) { + int sym_id = -1; + if (target.is_relative_detector_id()) { + sym_id = target.val(); + } else if (target.is_observable_id()) { + sym_id = num_detectors + target.val(); + } + if (sym_id != -1) { + roots_touched.insert(uf.find(sym_id)); + } + } + + for (int root : roots_touched) { + root_to_component[root].error_indices.push_back(i); + } + } + + std::vector components; + for (auto& pair : root_to_component) { + components.push_back(std::move(pair.second)); + } + + return components; +} + +} // namespace tesseract diff --git a/src/multi_pass/tanner_graph.h b/src/multi_pass/tanner_graph.h new file mode 100644 index 00000000..b1f61cfd --- /dev/null +++ b/src/multi_pass/tanner_graph.h @@ -0,0 +1,56 @@ +#ifndef TANNER_GRAPH_H +#define TANNER_GRAPH_H + +#include +#include +#include + +#include "stim.h" + +namespace tesseract { + +/** + * Represents an independent connected component of the Tanner graph. + */ +struct TannerComponent { + std::vector detectors; + std::vector observables; + std::vector error_indices; // Indices of instructions in the DEM + bool affects_observable = false; +}; + +/** + * Utility to analyze the Tanner graph of a DetectorErrorModel. + */ +class TannerGraph { + public: + /** + * Finds all connected components in the provided DetectorErrorModel. + * + * Assumes the DEM has been decomposed (errors affect only one component's symptoms). + * If an error bridges symptoms, they will be unioned into the same component. + */ + static std::vector find_components(const stim::DetectorErrorModel& dem); + + private: + struct UnionFind { + std::vector parent; + UnionFind(size_t n) { + parent.resize(n); + for (size_t i = 0; i < n; ++i) parent[i] = i; + } + int find(int i) { + if (parent[i] == i) return i; + return parent[i] = find(parent[i]); + } + void unite(int i, int j) { + int root_i = find(i); + int root_j = find(j); + if (root_i != root_j) parent[root_i] = root_j; + } + }; +}; + +} // namespace tesseract + +#endif // TANNER_GRAPH_H diff --git a/src/multi_pass/tanner_graph.test.cc b/src/multi_pass/tanner_graph.test.cc new file mode 100644 index 00000000..bf461c63 --- /dev/null +++ b/src/multi_pass/tanner_graph.test.cc @@ -0,0 +1,83 @@ +#include "tanner_graph.h" + +#include + +#include "gtest/gtest.h" + +using namespace tesseract; + +TEST(TannerGraphTest, SingleComponent) { + stim::DetectorErrorModel dem(R"DEM( + error(0.1) D0 D1 + error(0.1) D1 L0 + detector D0 + detector D1 + logical_observable L0 + )DEM"); + auto components = TannerGraph::find_components(dem); + ASSERT_EQ(components.size(), 1); + ASSERT_EQ(components[0].detectors.size(), 2); + ASSERT_EQ(components[0].observables.size(), 1); + ASSERT_TRUE(components[0].affects_observable); +} + +TEST(TannerGraphTest, TwoDisjointComponents) { + stim::DetectorErrorModel dem(R"DEM( + error(0.1) D0 D1 + error(0.1) D2 L0 + detector D0 + detector D1 + detector D2 + logical_observable L0 + )DEM"); + auto components = TannerGraph::find_components(dem); + ASSERT_EQ(components.size(), 2); + + int obs_comp_idx = components[0].affects_observable ? 0 : 1; + int other_comp_idx = 1 - obs_comp_idx; + + ASSERT_EQ(components[obs_comp_idx].detectors.size(), 1); // D2 + ASSERT_EQ(components[obs_comp_idx].observables.size(), 1); // L0 + + ASSERT_EQ(components[other_comp_idx].detectors.size(), 2); // D0, D1 + ASSERT_EQ(components[other_comp_idx].observables.size(), 0); + ASSERT_FALSE(components[other_comp_idx].affects_observable); +} + +TEST(TannerGraphTest, DecomposedErrorDoesNotUnion) { + stim::DetectorErrorModel dem(R"DEM( + error(0.1) D0 D1 ^ D2 D3 + detector D0 + detector D1 + detector D2 + detector D3 + )DEM"); + auto components = TannerGraph::find_components(dem); + // Should be two components: {D0, D1} and {D2, D3} + ASSERT_EQ(components.size(), 2); +} + +TEST(TannerGraphTest, UndecomposedBridgeUnions) { + stim::DetectorErrorModel dem(R"DEM( + error(0.1) D0 D1 D2 D3 + detector D0 + detector D1 + detector D2 + detector D3 + )DEM"); + auto components = TannerGraph::find_components(dem); + // Should be one component: {D0, D1, D2, D3} + ASSERT_EQ(components.size(), 1); +} + +TEST(TannerGraphTest, PureLogicalErrorComponent) { + stim::DetectorErrorModel dem(R"DEM( + error(0.1) L0 + logical_observable L0 + )DEM"); + auto components = TannerGraph::find_components(dem); + ASSERT_EQ(components.size(), 1); + ASSERT_EQ(components[0].detectors.size(), 0); + ASSERT_EQ(components[0].observables.size(), 1); + ASSERT_TRUE(components[0].affects_observable); +} diff --git a/src/py/BUILD b/src/py/BUILD index 7160760f..446e5c25 100644 --- a/src/py/BUILD +++ b/src/py/BUILD @@ -17,6 +17,18 @@ load("@rules_python//python:pip.bzl", "compile_pip_requirements") load("@rules_python//python:py_library.bzl", "py_library") load("@rules_python//python:py_binary.bzl", "py_binary") +py_library( + name = "multi_pass_sinter_decoders", + srcs = ["multi_pass_sinter_decoders.py"], + imports = ["."], + visibility = ["//visibility:public"], + deps = [ + "//src:lib_tesseract_decoder", + "@pypi//sinter", + "@pypi//stim", + ], +) + py_library( name = "shared_decoding_tests", srcs = ["shared_decoding_tests.py"], @@ -93,6 +105,20 @@ py_test( imports = ["..", "."], ) +py_test( + name = "multi_pass_bindings_test", + srcs = ["multi_pass_bindings_test.py"], + visibility = ["//:__subpackages__"], + deps = [ + "@pypi//pytest", + "@pypi//stim", + "@pypi//numpy", + "//src:lib_tesseract_decoder", + ":multi_pass_sinter_decoders", + ], + imports = ["..", "."], +) + py_test( diff --git a/src/py/README.md b/src/py/README.md index 658932a8..afbc2347 100644 --- a/src/py/README.md +++ b/src/py/README.md @@ -539,6 +539,60 @@ def get_tesseract_decoder_for_sinter(): return tesseract_module.make_tesseract_sinter_decoders_dict() ``` +#### Multi-pass Tesseract decoding + +`MultiPassSinterDecoder` partitions a detector error model into exactly two detector components. It +accepts one or two passes (default: 2) and uses causal scheduling by default. Its built-in classifier +checks strict `"X"`/`"Z"` `measure_basis` metadata, then `basis` metadata, then the fourth detector +coordinate (`0`–`2` or `3`–`5`). Every detector must be classified, and exactly two distinct labels +must result. + +Standard Tesseract options and multi-pass wrapper options can be passed directly as keyword +arguments: + +```python +import stim +import tesseract_decoder +from multi_pass_sinter_decoders import MultiPassSinterDecoder + +dem = stim.DetectorErrorModel(""" + error(0.1) D0 ^ D1 L0 + error(0.01) D0 + error(0.2) D1 L0 + detector[{"measure_basis": "X"}] D0 + detector[{"measure_basis": "Z"}] D1 + logical_observable L0 +""") + +decoder = MultiPassSinterDecoder( + num_passes=2, + det_beam=20, + beam_climbing=True, + pqlimit=1_000_000, + num_det_orders=21, + det_order_method=tesseract_decoder.utils.DetOrder.DetIndex, +) +compiled_decoder = decoder.compile_decoder_for_dem(dem=dem) +``` + +For other annotation conventions, supply a callable with signature +`(detector_index, coordinates, tag) -> component_label`. Returning a negative label rejects a +detector that the callable cannot classify. For example, the built-in coordinate convention can be +expressed as: + +```python +coordinate_classifier = lambda _index, coordinates, _tag: ( + 0 if len(coordinates) >= 4 and 0 <= coordinates[3] <= 2 + else 1 if len(coordinates) >= 4 and 3 <= coordinates[3] <= 5 + else -1 +) +decoder = MultiPassSinterDecoder(detector_classifier=coordinate_classifier) +``` + +`get_sinter_decoders()` provides long-beam monolithic, one-pass, and two-pass configurations using +the same Tesseract settings; they differ only in whether multi-pass decoding is enabled and in the +number of passes. + #### Decoding with `sinter.collect` `sinter.collect` is a powerful function for running many decoding jobs in parallel and collecting the results for large-scale benchmarking. diff --git a/src/py/multi_pass_bindings_test.py b/src/py/multi_pass_bindings_test.py new file mode 100644 index 00000000..f846c882 --- /dev/null +++ b/src/py/multi_pass_bindings_test.py @@ -0,0 +1,102 @@ +import tesseract_decoder +import stim +import numpy as np +import sys +from multi_pass_sinter_decoders import ( + MultiPassSinterDecoder as PythonMultiPassSinterDecoder, + get_sinter_decoders, +) + +def test_multi_pass_sinter_bindings(): + print(f"Loaded tesseract_decoder from: {tesseract_decoder.__file__}", flush=True) + + dem = stim.DetectorErrorModel(R""" + error(0.1) D0 ^ D1 L0 + error(0.01) D0 + error(0.2) D1 L0 + detector D0 + detector D1 + logical_observable L0 + """) + + # 1. Test with Detector Classifier Lambda + print("Testing MultiPassSinterDecoder with lambda...", flush=True) + decoder = tesseract_decoder.MultiPassSinterDecoder(num_passes=2) + assert decoder.strategy == tesseract_decoder.Causal + det_index = tesseract_decoder.utils.DetOrder.DetIndex + assert decoder.det_order_method == det_index + decoder.detector_classifier = lambda index, coords, tag: index + + assert PythonMultiPassSinterDecoder().strategy == tesseract_decoder.Causal + python_static_decoder = PythonMultiPassSinterDecoder(strategy=tesseract_decoder.Static) + assert python_static_decoder.strategy == tesseract_decoder.Static + + registered_decoders = get_sinter_decoders() + assert registered_decoders["tesseract-long-beam-mono"].det_order_method == det_index + assert ( + registered_decoders["tesseract-long-beam-multipass-1pass"].base_config_kwargs[ + "det_order_method" + ] + == det_index + ) + assert ( + registered_decoders["tesseract-long-beam-multipass-2pass"].base_config_kwargs[ + "det_order_method" + ] + == det_index + ) + + strict_dem = stim.DetectorErrorModel(R""" + error(0.1) D0 + error(0.2) D1 L0 + detector[{"measure_basis": "X"}] D0 + detector[{"md": {"basis": "Z"}}] D1 + logical_observable L0 + """) + PythonMultiPassSinterDecoder().compile_decoder_for_dem(dem=strict_dem) + + invalid_dem = stim.DetectorErrorModel(R""" + error(0.1) D0 + error(0.2) D1 L0 + detector[{"measure_basis": 0, "basis": "X"}](0, 0, 0, 0) D0 + detector[{"measure_basis": "Z"}] D1 + logical_observable L0 + """) + try: + PythonMultiPassSinterDecoder().compile_decoder_for_dem(dem=invalid_dem) + raise AssertionError("Expected invalid measure_basis to be rejected") + except ValueError as error: + assert "could not be classified" in str(error) + + compiled = decoder.compile_decoder_for_dem(dem=dem) + + # D0 and D1 both fire. Bit-packed: 0b11 = 3 + dets = np.array([[3]], dtype=np.uint8) + predictions = compiled.decode_shots_bit_packed(bit_packed_detection_event_data=dets) + + print(f"Predictions: {predictions}", flush=True) + assert (predictions[0, 0] & 1) == 1 + + # 2. A decomposer does not replace the required detector classification. + print("Testing missing classifier rejection...", flush=True) + def my_decomposer(input_dem): + print("Full decomposer called!", flush=True) + return input_dem + + decoder.detector_classifier = None + decoder.full_decomposer = my_decomposer + try: + decoder.compile_decoder_for_dem(dem=dem) + raise AssertionError("Expected detector_classifier to be required") + except ValueError as error: + assert "detector_classifier" in str(error) + +if __name__ == "__main__": + try: + test_multi_pass_sinter_bindings() + print("Python bindings test PASSED", flush=True) + except Exception as e: + print(f"Python bindings test FAILED: {e}", flush=True) + import traceback + traceback.print_exc() + sys.exit(1) diff --git a/src/py/multi_pass_sinter_decoders.py b/src/py/multi_pass_sinter_decoders.py new file mode 100644 index 00000000..c9275789 --- /dev/null +++ b/src/py/multi_pass_sinter_decoders.py @@ -0,0 +1,122 @@ +import sinter +import stim +import tesseract_decoder as _core + +class MultiPassSinterDecoder(sinter.Decoder): + """A Sinter-compatible wrapper around the native multi-pass Tesseract decoder. + + Args: + num_passes: Number of passes. Only 1 or 2 are supported. + detector_classifier: Optional ``(index, coordinates, tag) -> int`` callable. It must + assign every detector a nonnegative label and produce exactly two distinct labels. + By default, X/Z metadata tags and then Chromobius-style coordinates are used. + strategy: Pass scheduling strategy. Defaults to causal scheduling. + **base_config_kwargs: Standard Tesseract options such as ``det_beam``, ``pqlimit``, + ``beam_climbing``, ``no_revisit_dets``, and ``merge_errors``. Multi-pass wrapper + options such as ``num_det_orders``, ``det_order_method``, and ``seed`` are also + accepted. + """ + def __init__(self, num_passes: int = 2, detector_classifier=None, + strategy=_core.Causal, **base_config_kwargs): + if num_passes not in (1, 2): + raise ValueError("num_passes must be 1 or 2.") + self.num_passes = num_passes + self.detector_classifier = detector_classifier + self.strategy = strategy + self.base_config_kwargs = base_config_kwargs + + def compile_decoder_for_dem(self, *, dem: stim.DetectorErrorModel) -> sinter.CompiledDecoder: + # 1. Access the native C++ class + cpp_decoder = _core.MultiPassSinterDecoder(num_passes=self.num_passes) + cpp_decoder.strategy = self.strategy + + # 2. Attach the classifier if provided + if self.detector_classifier is not None: + cpp_decoder.detector_classifier = self.detector_classifier + else: + def default_classifier(index: int, coords: list[float], tag: str) -> int: + import json + # Priority 1: Parse JSON tag for "measure_basis" then "basis". + # Supports both top-level keys and keys nested under "md". + if tag: + try: + tag_data = json.loads(tag) + if isinstance(tag_data, dict): + md = tag_data.get("md", {}) + if not isinstance(md, dict): + md = {} + + basis_fields = ( + (tag_data, "measure_basis"), + (md, "measure_basis"), + (tag_data, "basis"), + (md, "basis"), + ) + for metadata, key in basis_fields: + if key not in metadata: + continue + if metadata[key] == "X": + return 0 + if metadata[key] == "Z": + return 1 + return -1 + except json.JSONDecodeError: + pass + # Priority 3: Chromobius-style coordinate convention. + if len(coords) >= 4: + c3 = int(coords[3]) + if 0 <= c3 <= 2: + return 0 + if 3 <= c3 <= 5: + return 1 + return -1 + cpp_decoder.detector_classifier = default_classifier + + # 3. Apply base configuration (pqlimit, det_beam, etc.) + for key, value in self.base_config_kwargs.items(): + if hasattr(cpp_decoder.base_config, key): + setattr(cpp_decoder.base_config, key, value) + elif hasattr(cpp_decoder, key): + setattr(cpp_decoder, key, value) + + # 4. Compile and return the native CompiledDecoder + return cpp_decoder.compile_decoder_for_dem(dem=dem) + +def get_sinter_decoders(): + TesseractSinterDecoder = _core.TesseractSinterDecoder + return { + "tesseract-long-beam-mono": TesseractSinterDecoder( + det_beam=20, + beam_climbing=True, + no_revisit_dets=True, + merge_errors=True, + pqlimit=1000000, + num_det_orders=21, + det_order_method=_core.utils.DetOrder.DetIndex, + seed=2384753 + ), + "tesseract-long-beam-multipass-1pass": MultiPassSinterDecoder( + num_passes=1, + strategy=_core.Causal, + det_beam=20, + beam_climbing=True, + no_revisit_dets=True, + merge_errors=True, + pqlimit=1000000, + num_det_orders=21, + det_order_method=_core.utils.DetOrder.DetIndex, + seed=2384753 + ), + "tesseract-long-beam-multipass-2pass": MultiPassSinterDecoder( + num_passes=2, + strategy=_core.Causal, + det_beam=20, + beam_climbing=True, + no_revisit_dets=True, + merge_errors=True, + pqlimit=1000000, + num_det_orders=21, + det_order_method=_core.utils.DetOrder.DetIndex, + seed=2384753 + ), + } diff --git a/src/py/stub_test.py b/src/py/stub_test.py index 7721e6f4..b380b883 100644 --- a/src/py/stub_test.py +++ b/src/py/stub_test.py @@ -126,6 +126,9 @@ class TestStubContents: "TesseractDecoder", "TesseractSinterCompiledDecoder", "TesseractSinterDecoder", + "MultiPassSinterCompiledDecoder", + "MultiPassSinterDecoder", + "SchedulingStrategy", "SimplexConfig", "SimplexDecoder", "DetOrder", @@ -144,4 +147,4 @@ def test_expected_symbol_in_stubs(self, stub_files, symbol): if __name__ == "__main__": - raise SystemExit(pytest.main([__file__])) \ No newline at end of file + raise SystemExit(pytest.main([__file__])) diff --git a/src/tesseract.cc b/src/tesseract.cc index 0a180754..e5c65b02 100644 --- a/src/tesseract.cc +++ b/src/tesseract.cc @@ -202,6 +202,32 @@ TesseractDecoder::TesseractDecoder(TesseractConfig config_) : config(std::move(c } } +void TesseractDecoder::update_internal_costs(const std::vector& modified_error_indices) { + std::unordered_set affected_detectors; + + for (size_t ei : modified_error_indices) { + if (ei >= errors.size()) continue; + double min_cost = errors[ei].symptom.detectors.empty() + ? errors[ei].likelihood_cost + : errors[ei].likelihood_cost / errors[ei].symptom.detectors.size(); + error_costs[ei] = {errors[ei].likelihood_cost, min_cost}; + + // Collect all detectors affected by this error to re-sort their d2e lists + for (int d : edets[ei]) { + affected_detectors.insert(d); + } + } + + // Re-sort d2e lists only for affected detectors + for (int d : affected_detectors) { + if (d >= 0 && (size_t)d < d2e.size()) { + std::sort(d2e[d].begin(), d2e[d].end(), [this](size_t idx_a, size_t idx_b) { + return error_costs[idx_a].min_cost < error_costs[idx_b].min_cost; + }); + } + } +} + void TesseractDecoder::initialize_structures(size_t num_detectors) { d2e.resize(num_detectors); edets.resize(num_errors); @@ -213,6 +239,8 @@ void TesseractDecoder::initialize_structures(size_t num_detectors) { } } + // Initial fill of error_costs and sorting of d2e for all errors + error_costs.reserve(errors.size()); for (size_t i = 0; i < errors.size(); ++i) { error_costs.push_back({errors[i].likelihood_cost, errors[i].likelihood_cost / errors[i].symptom.detectors.size()}); @@ -291,6 +319,11 @@ void TesseractDecoder::initialize_structures(size_t num_detectors) { } void TesseractDecoder::decode_to_errors(const std::vector& detections) { + predicted_errors_buffer.clear(); + low_confidence_flag = false; + if (detections.empty()) { + return; + } if (config.sparsify_errors) { build_sparse_d2e(detections); } diff --git a/src/tesseract.h b/src/tesseract.h index 97b88eb4..26f1e025 100644 --- a/src/tesseract.h +++ b/src/tesseract.h @@ -100,6 +100,15 @@ struct TesseractDecoder { // flattened DEM error indices. double cost_from_errors(const std::vector& predicted_errors) const; + // Resynchronizes the internal state of the decoder after the public `errors` + // vector has been modified. This is necessary to ensure that the internal + // cost structures used by the decoding algorithm are consistent with the + // current error likelihoods. + // This is necessary to ensure that the internal + // cost structures used by the decoding algorithm are consistent with the + // current error likelihoods. + void update_internal_costs(const std::vector& modified_error_indices); + std::vector decode(const std::vector& detections); void decode_shots(std::vector& shots, std::vector>& obs_predicted); @@ -142,6 +151,18 @@ struct TesseractDecoder { void decode_to_errors_with_graph(const std::vector& detections, size_t detector_order, size_t detector_beam, const std::vector>& active_d2e); + + friend class TesseractDebugger; +}; + +class TesseractDebugger { + public: + static const std::vector& get_error_costs(const TesseractDecoder& decoder) { + return decoder.error_costs; + } + static const std::vector>& get_d2e(const TesseractDecoder& decoder) { + return decoder.d2e; + } }; #endif // TESSERACT_DECODER_H diff --git a/src/tesseract.pybind.cc b/src/tesseract.pybind.cc index 9f2808f4..3af4a4ab 100644 --- a/src/tesseract.pybind.cc +++ b/src/tesseract.pybind.cc @@ -18,6 +18,7 @@ #include #include "common.pybind.h" +#include "multi_pass/multi_pass_sinter_compat.pybind.h" #include "pybind11/detail/common.h" #include "simplex.pybind.h" #include "tesseract_sinter_compat.pybind.h" @@ -33,6 +34,7 @@ PYBIND11_MODULE(tesseract_decoder, tesseract) { add_visualization_module(tesseract); add_tesseract_module(tesseract); pybind_sinter_compat(tesseract); + tesseract::pybind_multi_pass_sinter_compat(tesseract); tesseract.attr("demutil") = py::module::import("_tesseract_py_util"); // Adds a context manager to the python library that can be used to redirect C++'s stdout/stderr diff --git a/src/tesseract.test.cc b/src/tesseract.test.cc index fd0d471e..da4d01df 100644 --- a/src/tesseract.test.cc +++ b/src/tesseract.test.cc @@ -559,3 +559,40 @@ TEST(tesseract, MoreThan64Observables) { ASSERT_EQ(flipped[i], i); } } + +// Test to ensure update_internal_costs correctly reflects changes to error likelihoods +TEST(tesseract, UpdateInternalCostsBehavior) { + // Define a simple DEM with two errors that can explain detector D0 + // Error 0: D0 (prob 0.2) -> likelihood_cost: ~1.386 + // Error 1: D0 (prob 0.1) -> likelihood_cost: ~2.197 + // Initially, Error 0 is more likely (lower likelihood_cost) + stim::DetectorErrorModel dem(R"DEM( + error(0.2) D0 + error(0.1) D0 + detector(0,0,0) D0 + )DEM"); + + TesseractConfig config{dem}; + config.merge_errors = false; // Important: do not merge errors for this test + TesseractDecoder decoder(config); + + // Initial decode: D0 fires. Should pick Error 0 (index 0) as it's more likely. + std::vector detections = {0}; + decoder.decode_to_errors(detections); + ASSERT_EQ(decoder.predicted_errors_buffer.size(), 1); + ASSERT_EQ(decoder.predicted_errors_buffer[0], 0); // Should pick Error 0 (index 0) + + // Manually change the likelihood_cost of Error 1 to be lower (more likely) than Error 0 + // Original: Error 0 (prob 0.2, cost ~1.386), Error 1 (prob 0.1, cost ~2.197) + // Modify: Error 1 to prob 0.3 (cost ~0.847). Now Error 1 is more likely. + decoder.errors[1].set_with_probability(0.3); + + // Call update_internal_costs to re-synchronize the decoder's state + decoder.update_internal_costs({1}); + + // Decode again with the same detections. + // Now, D0 fires. It should pick Error 1 (index 1) as it's now more likely. + decoder.decode_to_errors(detections); + ASSERT_EQ(decoder.predicted_errors_buffer.size(), 1); + ASSERT_EQ(decoder.predicted_errors_buffer[0], 1); // Should now pick Error 1 (index 1) +} diff --git a/src/tesseract_main.cc b/src/tesseract_main.cc index 4b3975fc..69542ee0 100644 --- a/src/tesseract_main.cc +++ b/src/tesseract_main.cc @@ -25,11 +25,16 @@ #include #include "common.h" +#include "multi_pass/multi_pass_tesseract_decoder.h" #include "stim.h" #include "tesseract.h" #include "utils.h" struct Args { + bool multipass = false; + bool print_multipass_plan = false; + std::string multipass_strategy = "causal"; + size_t num_passes = 2; std::string circuit_path; std::string dem_path; bool no_merge_errors = false; @@ -40,6 +45,7 @@ struct Args { bool det_order_bfs = false; bool det_order_index = false; bool det_order_coordinate = false; + DetOrder det_order_method = DetOrder::DetIndex; // Sampling options size_t sample_num_shots = 0; @@ -131,6 +137,12 @@ struct Args { if (num_threads == 0) { throw std::invalid_argument("--threads must be at least 1."); } + if (num_passes < 1 || num_passes > 2) { + throw std::invalid_argument("--num-passes must be 1 or 2."); + } + if (print_multipass_plan && !multipass) { + throw std::invalid_argument("--print-multipass-plan requires --multipass."); + } if (num_threads > 1000) { throw std::invalid_argument( "There is a maximum limit of 1000 threads imposed to avoid " @@ -225,7 +237,7 @@ struct Args { std::cout << ")" << std::endl; } } - DetOrder order = DetOrder::DetIndex; + DetOrder order = det_order_method; if (det_order_bfs) { order = DetOrder::DetBFS; } else if (det_order_index) { @@ -233,6 +245,7 @@ struct Args { } else if (det_order_coordinate) { order = DetOrder::DetCoordinate; } + det_order_method = order; config.det_orders = build_det_orders(config.dem, num_det_orders, order, det_order_seed); } @@ -510,6 +523,26 @@ int main(int argc, char* argv[]) { "during decoding.") .flag() .store_into(args.print_stats); + program.add_argument("--multipass") + .help("Enable multi-pass graph shattering for correlated error decoding") + .flag() + .store_into(args.multipass); + program.add_argument("--print-multipass-plan") + .help("Print the multi-pass components, dependencies, and schedule to stderr") + .flag() + .store_into(args.print_multipass_plan); + program.add_argument("--multipass-strategy", "--multipass_strategy") + .help( + "Multi-pass scheduling strategy: static or causal (default = causal). Note: static " + "scheduling is experimental and was never systematically benchmarked.") + .default_value(std::string("causal")) + .store_into(args.multipass_strategy); + program.add_argument("--num-passes", "--num_passes") + .help( + "Number of prior propagation passes: 1 (uncorrelated independent CSS decoding) or 2 " + "(standard causally reweighted decoding, default = 2).") + .default_value(size_t(2)) + .store_into(args.num_passes); program.add_argument("--sparsify-errors") .help("Enables per-shot sparse error activation.") @@ -537,12 +570,12 @@ int main(int argc, char* argv[]) { try { program.parse_args(argc, argv); + args.validate(program); } catch (const std::exception& err) { std::cerr << err.what() << std::endl; std::cerr << program; return EXIT_FAILURE; } - args.validate(program); TesseractConfig config; std::vector shots; std::unique_ptr writer; @@ -555,35 +588,134 @@ int main(int argc, char* argv[]) { std::vector> low_confidence(shots.size()); const stim::DetectorErrorModel original_dem = config.dem.flattened(); std::vector> decoders(args.num_threads); + std::vector> mp_decoders(args.num_threads); std::vector> error_use_per_thread( args.num_threads, std::vector(original_dem.count_errors())); bool has_obs = args.has_observables(); - size_t num_errors = 0; - size_t num_low_confidence = 0; - double total_time_seconds = 0; + std::atomic num_errors(0); + std::atomic num_low_confidence(0); + std::atomic total_time_seconds(0); + + if (args.multipass && !args.dem_out_fname.empty()) { + std::cerr << "Error: --dem-out is not supported when --multipass is enabled." << std::endl; + return 1; + } + + if (args.multipass_strategy != "static" && args.multipass_strategy != "causal") { + std::cerr << "Error: Invalid --multipass-strategy '" << args.multipass_strategy + << "'. Expected 'static' or 'causal'." << std::endl; + return 1; + } + + auto classifier = [](int index, const std::vector& coords, + const std::string& tag) -> int { + // Priority 1: Parse JSON tag for "measure_basis" then "basis". + // Supports both top-level keys and keys nested under "md". + if (!tag.empty()) { + try { + auto tag_data = nlohmann::json::parse(tag); + + auto classify_tag = [](const nlohmann::json& metadata, const char* key, int& component) { + if (!metadata.is_object() || !metadata.contains(key)) return false; + const auto& value = metadata[key]; + component = value == "X" ? 0 : value == "Z" ? 1 : -1; + return true; + }; + const nlohmann::json empty_metadata; + const auto& md = + tag_data.is_object() && tag_data.contains("md") ? tag_data["md"] : empty_metadata; + int component; + if (classify_tag(tag_data, "measure_basis", component)) return component; + if (classify_tag(md, "measure_basis", component)) return component; + if (classify_tag(tag_data, "basis", component)) return component; + if (classify_tag(md, "basis", component)) return component; + } catch (const nlohmann::json::parse_error&) { + // Tag is not valid JSON; fall through to coordinate-based classification. + } + } + + // Priority 3: Chromobius-style coordinate convention. + // 4th coordinate encodes color + 3*basis: 0-2 => X (component 0), 3-5 => Z (component 1). + if (coords.size() >= 4) { + int c3 = (int)coords[3]; + if (c3 >= 0 && c3 <= 2) return 0; + if (c3 >= 3 && c3 <= 5) return 1; + } + return -1; + }; + tesseract::SchedulingStrategy strategy_val = (args.multipass_strategy == "static") + ? tesseract::SchedulingStrategy::Static + : tesseract::SchedulingStrategy::Causal; + + std::vector detector_classes; + if (args.multipass) { + try { + detector_classes = + tesseract::MultiPassTesseractDecoder::classify_detectors(config.dem, classifier); + if (args.print_multipass_plan) { + mp_decoders[0] = std::make_unique( + config.dem, args.num_passes, detector_classes, config, args.num_det_orders, + args.det_order_method, args.det_order_seed, strategy_val, true); + std::cerr << mp_decoders[0]->get_execution_plan().str(); + } + } catch (const std::invalid_argument& error) { + std::cerr << "Error: " << error.what() << std::endl; + return 1; + } + } + size_t shot = parallel_for_shots_in_order( shots.size(), args.num_threads, [&](size_t thread_index, size_t shot_index) { - if (!decoders[thread_index]) { - decoders[thread_index] = std::make_unique(config); - } - auto& decoder = *decoders[thread_index]; auto& error_use = error_use_per_thread[thread_index]; - auto start_time = std::chrono::high_resolution_clock::now(); - decoder.decode_to_errors(shots[shot_index].hits); - auto stop_time = std::chrono::high_resolution_clock::now(); - decoding_time_seconds[shot_index] = - std::chrono::duration_cast(stop_time - start_time).count() / - 1e6; - obs_predicted[shot_index].clear(); - for (int obs_idx : decoder.get_flipped_observables(decoder.predicted_errors_buffer)) { - obs_predicted[shot_index][obs_idx] ^= 1; - } - low_confidence[shot_index] = decoder.low_confidence_flag; - cost_predicted[shot_index] = decoder.cost_from_errors(decoder.predicted_errors_buffer); - if (!has_obs or shots[shot_index].obs_mask == obs_predicted[shot_index]) { - for (size_t ei : decoder.predicted_errors_buffer) { - ++error_use[ei]; + + if (args.multipass) { + if (!mp_decoders[thread_index]) { + mp_decoders[thread_index] = std::make_unique( + config.dem, args.num_passes, detector_classes, config, args.num_det_orders, + args.det_order_method, args.det_order_seed, strategy_val); + } + auto start_time = std::chrono::high_resolution_clock::now(); + auto res = mp_decoders[thread_index]->decode_result(shots[shot_index].hits); + auto stop_time = std::chrono::high_resolution_clock::now(); + decoding_time_seconds[shot_index] = + std::chrono::duration_cast(stop_time - start_time) + .count() / + 1e6; + + obs_predicted[shot_index].clear(); + for (int o : res.predictions) { + if (o >= 0 && (size_t)o < num_observables) { + obs_predicted[shot_index][o] ^= 1; + } + } + low_confidence[shot_index] = res.low_confidence; + cost_predicted[shot_index] = res.total_cost; + } else { + if (!decoders[thread_index]) { + decoders[thread_index] = std::make_unique(config); + } + auto& decoder = *decoders[thread_index]; + auto start_time = std::chrono::high_resolution_clock::now(); + decoder.decode_to_errors(shots[shot_index].hits); + auto stop_time = std::chrono::high_resolution_clock::now(); + decoding_time_seconds[shot_index] = + std::chrono::duration_cast(stop_time - start_time) + .count() / + 1e6; + + obs_predicted[shot_index].clear(); + for (int o : decoder.get_flipped_observables(decoder.predicted_errors_buffer)) { + if (o >= 0 && (size_t)o < num_observables) { + obs_predicted[shot_index][o] ^= 1; + } + } + low_confidence[shot_index] = decoder.low_confidence_flag; + cost_predicted[shot_index] = decoder.cost_from_errors(decoder.predicted_errors_buffer); + if (!has_obs || shots[shot_index].obs_mask == obs_predicted[shot_index]) { + for (size_t ei : decoder.predicted_errors_buffer) { + ++error_use[ei]; + } } } }, @@ -597,20 +729,20 @@ int main(int argc, char* argv[]) { } else if (has_obs && obs_predicted[shot_index] != shots[shot_index].obs_mask) { ++num_errors; } - total_time_seconds += decoding_time_seconds[shot_index]; + total_time_seconds = total_time_seconds + decoding_time_seconds[shot_index]; if (args.print_stats) { std::cout << "num_shots = " << (shot_index + 1) - << " num_low_confidence = " << num_low_confidence; + << " num_low_confidence = " << num_low_confidence.load(); if (has_obs) { - std::cout << " num_errors = " << num_errors; + std::cout << " num_errors = " << num_errors.load(); } - std::cout << " total_time_seconds = " << total_time_seconds << std::endl; + std::cout << " total_time_seconds = " << total_time_seconds.load() << std::endl; std::cout << "cost = " << cost_predicted[shot_index] << std::endl; std::cout.flush(); } // Disable early termination due to \`--max-errors\` when we don't have the ground-truth // observables - return !has_obs || num_errors < args.max_errors; + return !has_obs || num_errors.load() < args.max_errors; }); std::vector error_use_totals(original_dem.count_errors()); @@ -625,7 +757,7 @@ int main(int argc, char* argv[]) { size_t num_usage_dem_shots = shot; if (has_obs) { // When we know the obs, we only count non-error shots. - num_usage_dem_shots -= num_errors; + num_usage_dem_shots -= num_errors.load(); } stim::DetectorErrorModel est_dem = common::dem_from_counts(original_dem, counts, num_usage_dem_shots); @@ -667,11 +799,14 @@ int main(int argc, char* argv[]) { {"pqlimit", args.pqlimit}, {"num_det_orders", args.num_det_orders}, {"det_order_seed", args.det_order_seed}, - {"total_time_seconds", total_time_seconds}, - {"num_errors", has_obs ? nlohmann::json(num_errors) : nullptr}, - {"num_low_confidence", num_low_confidence}, + {"total_time_seconds", total_time_seconds.load()}, + {"num_errors", has_obs ? nlohmann::json(num_errors.load()) : nullptr}, + {"num_low_confidence", num_low_confidence.load()}, {"num_shots", shot}, {"num_threads", args.num_threads}, + {"multipass", args.multipass}, + {"strategy", args.multipass_strategy}, + {"num_passes", args.num_passes}, {"sample_num_shots", args.sample_num_shots}, {"sparsify_errors", args.sparsify_errors}, {"sparsify_base_degree", args.sparsify_base_degree}, @@ -688,11 +823,11 @@ int main(int argc, char* argv[]) { } if (print_final_stats) { std::cout << "num_shots = " << shot; - std::cout << " num_low_confidence = " << num_low_confidence; + std::cout << " num_low_confidence = " << num_low_confidence.load(); if (has_obs) { - std::cout << " num_errors = " << num_errors; + std::cout << " num_errors = " << num_errors.load(); } - std::cout << " total_time_seconds = " << total_time_seconds; + std::cout << " total_time_seconds = " << total_time_seconds.load(); std::cout << std::endl; } }