From 1959cf4a7c057ff273eefb687a073b933f477fe6 Mon Sep 17 00:00:00 2001 From: Oscar Higgott Date: Tue, 19 Aug 2025 01:37:27 +0000 Subject: [PATCH 01/10] Add more exhaustive tests of correlations. Fix a bug where edge weight rewrites in the search graph weren't always undone. --- .../sparse_blossom/flooder/graph.cc | 60 +++++++++++++++++++ src/pymatching/sparse_blossom/flooder/graph.h | 1 + 2 files changed, 61 insertions(+) diff --git a/src/pymatching/sparse_blossom/flooder/graph.cc b/src/pymatching/sparse_blossom/flooder/graph.cc index 5c1b8850d..c0e690392 100644 --- a/src/pymatching/sparse_blossom/flooder/graph.cc +++ b/src/pymatching/sparse_blossom/flooder/graph.cc @@ -138,6 +138,66 @@ MatchingGraph::MatchingGraph(MatchingGraph&& graph) noexcept loaded_from_dem_without_correlations(graph.loaded_from_dem_without_correlations) { } +bool MatchingGraph::graph_structure_equal(const MatchingGraph& other) const { + if ((this->negative_weight_detection_events_set != other.negative_weight_detection_events_set) || + (this->negative_weight_observables_set != other.negative_weight_observables_set) || + (this->negative_weight_sum != other.negative_weight_sum) || + (this->is_user_graph_boundary_node != other.is_user_graph_boundary_node) || + (this->num_nodes != other.num_nodes) || (this->num_observables != other.num_observables) || + (this->normalising_constant != other.normalising_constant) || + (this->previous_weights != other.previous_weights) || + (this->edges_to_implied_weights_unconverted != other.edges_to_implied_weights_unconverted) || + (this->loaded_from_dem_without_correlations != other.loaded_from_dem_without_correlations)) { + return false; + } + if (this->nodes.size() != other.nodes.size()) { + return false; + } + for (size_t i = 0; i < nodes.size(); i++) { + if ((nodes[i].neighbors.size() != other.nodes[i].neighbors.size()) || + (nodes[i].neighbor_weights != other.nodes[i].neighbor_weights) || + (nodes[i].neighbor_implied_weights.size() != other.nodes[i].neighbor_implied_weights.size())) { + return false; + } + for (size_t j = 0; j < nodes[i].neighbors.size(); j++) { + if (j == 0 && nodes[i].neighbors[0] == nullptr && other.nodes[i].neighbors[0] == nullptr) { + continue; + } + auto node_i_j = nodes[i].neighbors[j] - nodes.data(); + auto other_node_i_j = other.nodes[i].neighbors[j] - other.nodes.data(); + if (node_i_j != other_node_i_j) { + // std::cout << " i: " << i << " j: " << j << " node_i_j: " << node_i_j << " other_node_i_j: " << + // other_node_i_j << std::endl; + return false; + } + } + for (size_t j = 0; j < nodes[i].neighbor_implied_weights.size(); j++) { + if (nodes[i].neighbor_implied_weights[j].size() != other.nodes[i].neighbor_implied_weights[j].size()) { + return false; + } + for (size_t k = 0; k < nodes[i].neighbor_implied_weights[j].size(); k++) { + if (*nodes[i].neighbor_implied_weights[j][k].edge0_ptr != + *other.nodes[i].neighbor_implied_weights[j][k].edge0_ptr) { + return false; + } + if (nodes[i].neighbor_implied_weights[j][k].implied_weight != + other.nodes[i].neighbor_implied_weights[j][k].implied_weight) { + return false; + } + if ((nodes[i].neighbor_implied_weights[j][k].edge1_ptr == nullptr) && + (other.nodes[i].neighbor_implied_weights[j][k].edge1_ptr == nullptr)) { + continue; + } + if (*nodes[i].neighbor_implied_weights[j][k].edge1_ptr != + *other.nodes[i].neighbor_implied_weights[j][k].edge1_ptr) { + return false; + } + } + } + } + return true; +} + MatchingGraph::MatchingGraph() : negative_weight_sum(0), num_nodes(0), num_observables(0), normalising_constant(0) { } diff --git a/src/pymatching/sparse_blossom/flooder/graph.h b/src/pymatching/sparse_blossom/flooder/graph.h index abfddfc25..7a0e0603b 100644 --- a/src/pymatching/sparse_blossom/flooder/graph.h +++ b/src/pymatching/sparse_blossom/flooder/graph.h @@ -67,6 +67,7 @@ class MatchingGraph { MatchingGraph(size_t num_nodes, size_t num_observables); MatchingGraph(size_t num_nodes, size_t num_observables, double normalising_constant); MatchingGraph(MatchingGraph&& graph) noexcept; + bool graph_structure_equal(const MatchingGraph& other) const; void add_edge( size_t u, size_t v, From 424cd8904bec770a461e652506732f48945d3e11 Mon Sep 17 00:00:00 2001 From: Oscar Higgott Date: Tue, 19 Aug 2025 01:41:08 +0000 Subject: [PATCH 02/10] remove graph structure comparison method from MatchingGraph --- .../sparse_blossom/flooder/graph.cc | 60 ------------------- src/pymatching/sparse_blossom/flooder/graph.h | 1 - 2 files changed, 61 deletions(-) diff --git a/src/pymatching/sparse_blossom/flooder/graph.cc b/src/pymatching/sparse_blossom/flooder/graph.cc index c0e690392..5c1b8850d 100644 --- a/src/pymatching/sparse_blossom/flooder/graph.cc +++ b/src/pymatching/sparse_blossom/flooder/graph.cc @@ -138,66 +138,6 @@ MatchingGraph::MatchingGraph(MatchingGraph&& graph) noexcept loaded_from_dem_without_correlations(graph.loaded_from_dem_without_correlations) { } -bool MatchingGraph::graph_structure_equal(const MatchingGraph& other) const { - if ((this->negative_weight_detection_events_set != other.negative_weight_detection_events_set) || - (this->negative_weight_observables_set != other.negative_weight_observables_set) || - (this->negative_weight_sum != other.negative_weight_sum) || - (this->is_user_graph_boundary_node != other.is_user_graph_boundary_node) || - (this->num_nodes != other.num_nodes) || (this->num_observables != other.num_observables) || - (this->normalising_constant != other.normalising_constant) || - (this->previous_weights != other.previous_weights) || - (this->edges_to_implied_weights_unconverted != other.edges_to_implied_weights_unconverted) || - (this->loaded_from_dem_without_correlations != other.loaded_from_dem_without_correlations)) { - return false; - } - if (this->nodes.size() != other.nodes.size()) { - return false; - } - for (size_t i = 0; i < nodes.size(); i++) { - if ((nodes[i].neighbors.size() != other.nodes[i].neighbors.size()) || - (nodes[i].neighbor_weights != other.nodes[i].neighbor_weights) || - (nodes[i].neighbor_implied_weights.size() != other.nodes[i].neighbor_implied_weights.size())) { - return false; - } - for (size_t j = 0; j < nodes[i].neighbors.size(); j++) { - if (j == 0 && nodes[i].neighbors[0] == nullptr && other.nodes[i].neighbors[0] == nullptr) { - continue; - } - auto node_i_j = nodes[i].neighbors[j] - nodes.data(); - auto other_node_i_j = other.nodes[i].neighbors[j] - other.nodes.data(); - if (node_i_j != other_node_i_j) { - // std::cout << " i: " << i << " j: " << j << " node_i_j: " << node_i_j << " other_node_i_j: " << - // other_node_i_j << std::endl; - return false; - } - } - for (size_t j = 0; j < nodes[i].neighbor_implied_weights.size(); j++) { - if (nodes[i].neighbor_implied_weights[j].size() != other.nodes[i].neighbor_implied_weights[j].size()) { - return false; - } - for (size_t k = 0; k < nodes[i].neighbor_implied_weights[j].size(); k++) { - if (*nodes[i].neighbor_implied_weights[j][k].edge0_ptr != - *other.nodes[i].neighbor_implied_weights[j][k].edge0_ptr) { - return false; - } - if (nodes[i].neighbor_implied_weights[j][k].implied_weight != - other.nodes[i].neighbor_implied_weights[j][k].implied_weight) { - return false; - } - if ((nodes[i].neighbor_implied_weights[j][k].edge1_ptr == nullptr) && - (other.nodes[i].neighbor_implied_weights[j][k].edge1_ptr == nullptr)) { - continue; - } - if (*nodes[i].neighbor_implied_weights[j][k].edge1_ptr != - *other.nodes[i].neighbor_implied_weights[j][k].edge1_ptr) { - return false; - } - } - } - } - return true; -} - MatchingGraph::MatchingGraph() : negative_weight_sum(0), num_nodes(0), num_observables(0), normalising_constant(0) { } diff --git a/src/pymatching/sparse_blossom/flooder/graph.h b/src/pymatching/sparse_blossom/flooder/graph.h index 7a0e0603b..abfddfc25 100644 --- a/src/pymatching/sparse_blossom/flooder/graph.h +++ b/src/pymatching/sparse_blossom/flooder/graph.h @@ -67,7 +67,6 @@ class MatchingGraph { MatchingGraph(size_t num_nodes, size_t num_observables); MatchingGraph(size_t num_nodes, size_t num_observables, double normalising_constant); MatchingGraph(MatchingGraph&& graph) noexcept; - bool graph_structure_equal(const MatchingGraph& other) const; void add_edge( size_t u, size_t v, From 92cc8a63d70dd39ad216d6127a6f8e3e82dab5e6 Mon Sep 17 00:00:00 2001 From: Oscar Higgott Date: Tue, 19 Aug 2025 06:43:31 +0000 Subject: [PATCH 03/10] use -log(p) weights and only add edges to the graph if they are not a component in a decomposed error with more than one component --- .../driver/mwpm_decoding.test.cc | 2 +- .../sparse_blossom/driver/user_graph.cc | 40 +++++++++++------- .../sparse_blossom/driver/user_graph.h | 41 ++++++++----------- .../sparse_blossom/driver/user_graph.test.cc | 29 ++++++------- 4 files changed, 56 insertions(+), 56 deletions(-) diff --git a/src/pymatching/sparse_blossom/driver/mwpm_decoding.test.cc b/src/pymatching/sparse_blossom/driver/mwpm_decoding.test.cc index a40ca1fbf..3ac4bff34 100644 --- a/src/pymatching/sparse_blossom/driver/mwpm_decoding.test.cc +++ b/src/pymatching/sparse_blossom/driver/mwpm_decoding.test.cc @@ -330,7 +330,7 @@ TEST(MwpmCorrelatedDecoding, BetterLogicalErrorRateThanUncorrelated) { size_t num_mistakes_uncorrelated = 0; size_t num_shots = 0; size_t max_shots = 300; - size_t expected_mistakes_correlated = 5; + size_t expected_mistakes_correlated = 4; size_t expected_mistakes_uncorrelated = 11; pm::ExtendedMatchingResult res(mwpm_correlated.flooder.graph.num_observables); std::vector edges; diff --git a/src/pymatching/sparse_blossom/driver/user_graph.cc b/src/pymatching/sparse_blossom/driver/user_graph.cc index 19734730a..5b7eca37f 100644 --- a/src/pymatching/sparse_blossom/driver/user_graph.cc +++ b/src/pymatching/sparse_blossom/driver/user_graph.cc @@ -17,6 +17,19 @@ #include "pymatching/rand/rand_gen.h" #include "pymatching/sparse_blossom/driver/implied_weights.h" +namespace { + +double bernoulli_xor(double p1, double p2) { + return p1 * (1 - p2) + p2 * (1 - p1); +} + +} // namespace + + +double pm::to_weight_for_correlations(double probability) { + return -std::log(probability); +} + double pm::merge_weights(double a, double b) { auto sgn = std::copysign(1, a) * std::copysign(1, b); auto signed_min = sgn * std::min(std::abs(a), std::abs(b)); @@ -342,6 +355,15 @@ void pm::UserGraph::handle_dem_instruction( } } +void pm::UserGraph::handle_dem_instruction_include_correlations( + double p, const std::vector& detectors, const std::vector& observables) { + if (detectors.size() == 2) { + add_or_merge_edge(detectors[0], detectors[1], observables, pm::to_weight_for_correlations(p), p, INDEPENDENT); + } else if (detectors.size() == 1) { + add_or_merge_boundary_edge(detectors[0], observables, pm::to_weight_for_correlations(p), p, INDEPENDENT); + } +} + void pm::UserGraph::get_nodes_on_shortest_path_from_source(size_t src, size_t dst, std::vector& out_nodes) { auto& mwpm = get_mwpm_with_search_graph(); bool src_is_boundary = is_boundary_node(src); @@ -444,18 +466,6 @@ double pm::UserGraph::get_edge_weight_normalising_constant(size_t max_num_distin } } -namespace { - -double bernoulli_xor(double p1, double p2) { - return p1 * (1 - p2) + p2 * (1 - p1); -} - -double to_weight(double probability) { - return std::log((1 - probability) / probability); -} - -} // namespace - void pm::add_decomposed_error_to_joint_probabilities( DecomposedDemError& error, std::map, std::map, double>>& joint_probabilites) { @@ -490,7 +500,7 @@ pm::UserGraph pm::detector_error_model_to_user_graph( pm::iter_dem_instructions_include_correlations( detector_error_model, [&](double p, const std::vector& detectors, std::vector& observables) { - user_graph.handle_dem_instruction(p, detectors, observables); + user_graph.handle_dem_instruction_include_correlations(p, detectors, observables); }, joint_probabilites); @@ -525,8 +535,8 @@ void pm::UserGraph::populate_implied_edge_weights( // error, would lead to a negatively weighted error. We do not support this (yet), and use a // minimum of 0.5 as an implied probability for an edge to be reweighted. double implied_probability_for_other_edge = - std::min(0.5, affected_edge_and_probability.second / marginal_probability); - double w = to_weight(implied_probability_for_other_edge); + std::min(0.9, affected_edge_and_probability.second / marginal_probability); + double w = pm::to_weight_for_correlations(implied_probability_for_other_edge); ImpliedWeightUnconverted implied{affected_edge.first, affected_edge.second, w}; edge.implied_weights_for_other_edges.push_back(implied); } diff --git a/src/pymatching/sparse_blossom/driver/user_graph.h b/src/pymatching/sparse_blossom/driver/user_graph.h index ab945cb4b..98fb7d513 100644 --- a/src/pymatching/sparse_blossom/driver/user_graph.h +++ b/src/pymatching/sparse_blossom/driver/user_graph.h @@ -120,6 +120,7 @@ class UserGraph { Mwpm& get_mwpm(); Mwpm& get_mwpm_with_search_graph(); void handle_dem_instruction(double p, const std::vector& detectors, const std::vector& observables); + void handle_dem_instruction_include_correlations(double p, const std::vector& detectors, const std::vector& observables); void get_nodes_on_shortest_path_from_source(size_t src, size_t dst, std::vector& out_nodes); void populate_implied_edge_weights( std::map, std::map, double>>& joint_probabilites); @@ -131,6 +132,8 @@ class UserGraph { bool _all_edges_have_error_probabilities; }; +double to_weight_for_correlations(double probability); + template inline double UserGraph::iter_discretized_edges( pm::weight_int num_distinct_weights, @@ -311,35 +314,22 @@ void iter_dem_instructions_include_correlations( component->observable_indices.push_back(target.val()); } else if (target.is_separator()) { instruction_contains_separator = true; - // If the previous error in the decomposition had 3 or more detectors, we throw an exception. - if (num_component_detectors > 2) { - throw std::invalid_argument( - "Encountered a decomposed error instruction with a hyperedge component (3 or more detectors). " - "This is not supported."); - } else if (num_component_detectors == 0) { + // We cannot have num_component_detectors > 2 at this point, or we would have already thrown an exception + if (num_component_detectors == 0) { throw std::invalid_argument( "Encountered a decomposed error instruction with an undetectable component (0 detectors). " "This is not supported."); - } else if (num_component_detectors > 0) { - // If the previous error in the decomposition had 1 or 2 detectors, we handle it - handle_dem_error(p, {component->node1, component->node2}, component->observable_indices); - decomposed_err.components.push_back({}); - component = &decomposed_err.components.back(); - component->node1 = SIZE_MAX; - component->node2 = SIZE_MAX; - num_component_detectors = 0; } + // The previous error in the decomposition must have 1 or 2 detectors + decomposed_err.components.push_back({}); + component = &decomposed_err.components.back(); + component->node1 = SIZE_MAX; + component->node2 = SIZE_MAX; + num_component_detectors = 0; } } - if (num_component_detectors > 2) { - // Undecomposed hyperedges are not supported - throw std::invalid_argument( - "Encountered an undecomposed error instruction with 3 or mode detectors. " - "This is not supported when using `enable_correlations=True`. " - "Did you forget to set `decompose_errors=True` when " - "converting the stim circuit to a detector error model?"); - } else if (num_component_detectors == 0) { + if (num_component_detectors == 0) { if (instruction_contains_separator) { throw std::invalid_argument( "Encountered a decomposed error instruction with an undetectable component (0 detectors). " @@ -348,8 +338,11 @@ void iter_dem_instructions_include_correlations( // Ignore errors that are undetectable, provided they are not a component of a decomposed error return; } - - } else if (num_component_detectors > 0) { + } + + // Only add the edge into the graph if it is not a component in a decomposed error with more than one component + if (decomposed_err.components.size() == 1) { + component = &decomposed_err.components.back(); if (component->node2 == SIZE_MAX) { handle_dem_error(p, {component->node1}, component->observable_indices); } else { diff --git a/src/pymatching/sparse_blossom/driver/user_graph.test.cc b/src/pymatching/sparse_blossom/driver/user_graph.test.cc index cd08bf18e..4b4f6bed3 100644 --- a/src/pymatching/sparse_blossom/driver/user_graph.test.cc +++ b/src/pymatching/sparse_blossom/driver/user_graph.test.cc @@ -320,9 +320,10 @@ TEST(IterDemInstructionsTest, DecomposedError) { pm::iter_dem_instructions_include_correlations(dem, handler, joint_probabilities); // Check handler calls - ASSERT_EQ(handler.handled_errors.size(), 3); - std::vector expected_handled = {{0.1, 0, 1, {}}, {0.1, 2, 3, {0}}, {0.1, 4, SIZE_MAX, {}}}; - EXPECT_EQ(handler.handled_errors, expected_handled); + ASSERT_EQ(handler.handled_errors.size(), 0); + // ASSERT_EQ(handler.handled_errors.size(), 3); + // std::vector expected_handled = {{0.1, 0, 1, {}}, {0.1, 2, 3, {0}}, {0.1, 4, SIZE_MAX, {}}}; + // EXPECT_EQ(handler.handled_errors, expected_handled); // Check joint probabilities std::pair key01 = {0, 1}; @@ -383,11 +384,14 @@ TEST(IterDemInstructionsTest, CombinedComplexDem) { std::map, std::map, double>> joint_probabilities; pm::iter_dem_instructions_include_correlations(dem, handler, joint_probabilities); - ASSERT_EQ(handler.handled_errors.size(), 4); - + ASSERT_EQ(handler.handled_errors.size(), 2); std::vector expected = { - {0.1, 0, SIZE_MAX, {}}, {0.2, 1, 2, {0}}, {0.4, 8, SIZE_MAX, {}}, {0.4, 9, SIZE_MAX, {1}}}; + {0.1, 0, SIZE_MAX, {}}, {0.2, 1, 2, {0}}}; EXPECT_EQ(handler.handled_errors, expected); + // ASSERT_EQ(handler.handled_errors.size(), 4); + // std::vector expected = { + // {0.1, 0, SIZE_MAX, {}}, {0.2, 1, 2, {0}}, {0.4, 8, SIZE_MAX, {}}, {0.4, 9, SIZE_MAX, {1}}}; + // EXPECT_EQ(handler.handled_errors, expected); // Check joint probabilities std::pair key0B = {0, SIZE_MAX}; @@ -485,14 +489,6 @@ TEST(UserGraph, PopulateImpliedEdgeWeights) { graph.populate_implied_edge_weights(joint_probabilities); - auto to_weight = [](double p) { - if (p == 1.0) - return -std::numeric_limits::infinity(); - if (p == 0.0) - return std::numeric_limits::infinity(); - return std::log((1 - p) / p); - }; - auto it_01 = std::find_if(graph.edges.begin(), graph.edges.end(), [](const pm::UserEdge& edge) { return edge.node1 == 0 && edge.node2 == 1; }); @@ -503,7 +499,7 @@ TEST(UserGraph, PopulateImpliedEdgeWeights) { ASSERT_EQ(implied_01.node2, 3); double p_01 = 0.1 / 0.26; - double w_01 = to_weight(p_01); + double w_01 = pm::to_weight_for_correlations(p_01); ASSERT_EQ(implied_01.implied_weight, w_01); auto it_23 = std::find_if(graph.edges.begin(), graph.edges.end(), [](const pm::UserEdge& edge) { @@ -515,7 +511,8 @@ TEST(UserGraph, PopulateImpliedEdgeWeights) { const auto& implied_23 = it_23->implied_weights_for_other_edges[0]; ASSERT_EQ(implied_23.node1, 0); ASSERT_EQ(implied_23.node2, 1); - ASSERT_EQ(implied_23.implied_weight, 0); + // ASSERT_EQ(implied_23.implied_weight, 0); + ASSERT_NEAR(implied_23.implied_weight, 0.105361, 0.00001); } TEST(UserGraph, ConvertImpliedWeights) { From 13bbea6575edf8605e1b9ca0045d9d6e0669b66b Mon Sep 17 00:00:00 2001 From: Oscar Higgott Date: Mon, 25 Aug 2025 07:51:38 +0000 Subject: [PATCH 04/10] handle boundary edge component and refactor --- .../driver/mwpm_decoding.test.cc | 2 +- .../sparse_blossom/driver/user_graph.cc | 4 +- .../sparse_blossom/driver/user_graph.h | 155 +++++++++--------- .../sparse_blossom/driver/user_graph.test.cc | 2 +- 4 files changed, 84 insertions(+), 79 deletions(-) diff --git a/src/pymatching/sparse_blossom/driver/mwpm_decoding.test.cc b/src/pymatching/sparse_blossom/driver/mwpm_decoding.test.cc index 3ac4bff34..a5f66f363 100644 --- a/src/pymatching/sparse_blossom/driver/mwpm_decoding.test.cc +++ b/src/pymatching/sparse_blossom/driver/mwpm_decoding.test.cc @@ -330,7 +330,7 @@ TEST(MwpmCorrelatedDecoding, BetterLogicalErrorRateThanUncorrelated) { size_t num_mistakes_uncorrelated = 0; size_t num_shots = 0; size_t max_shots = 300; - size_t expected_mistakes_correlated = 4; + size_t expected_mistakes_correlated = 3; size_t expected_mistakes_uncorrelated = 11; pm::ExtendedMatchingResult res(mwpm_correlated.flooder.graph.num_observables); std::vector edges; diff --git a/src/pymatching/sparse_blossom/driver/user_graph.cc b/src/pymatching/sparse_blossom/driver/user_graph.cc index 5b7eca37f..313343978 100644 --- a/src/pymatching/sparse_blossom/driver/user_graph.cc +++ b/src/pymatching/sparse_blossom/driver/user_graph.cc @@ -27,7 +27,7 @@ double bernoulli_xor(double p1, double p2) { double pm::to_weight_for_correlations(double probability) { - return -std::log(probability); + return std::log((1 - probability) / probability); } double pm::merge_weights(double a, double b) { @@ -535,7 +535,7 @@ void pm::UserGraph::populate_implied_edge_weights( // error, would lead to a negatively weighted error. We do not support this (yet), and use a // minimum of 0.5 as an implied probability for an edge to be reweighted. double implied_probability_for_other_edge = - std::min(0.9, affected_edge_and_probability.second / marginal_probability); + std::min(0.5, affected_edge_and_probability.second / marginal_probability); double w = pm::to_weight_for_correlations(implied_probability_for_other_edge); ImpliedWeightUnconverted implied{affected_edge.first, affected_edge.second, w}; edge.implied_weights_for_other_edges.push_back(implied); diff --git a/src/pymatching/sparse_blossom/driver/user_graph.h b/src/pymatching/sparse_blossom/driver/user_graph.h index 98fb7d513..0c3b88a4c 100644 --- a/src/pymatching/sparse_blossom/driver/user_graph.h +++ b/src/pymatching/sparse_blossom/driver/user_graph.h @@ -120,7 +120,8 @@ class UserGraph { Mwpm& get_mwpm(); Mwpm& get_mwpm_with_search_graph(); void handle_dem_instruction(double p, const std::vector& detectors, const std::vector& observables); - void handle_dem_instruction_include_correlations(double p, const std::vector& detectors, const std::vector& observables); + void handle_dem_instruction_include_correlations( + double p, const std::vector& detectors, const std::vector& observables); void get_nodes_on_shortest_path_from_source(size_t src, size_t dst, std::vector& out_nodes); void populate_implied_edge_weights( std::map, std::map, double>>& joint_probabilites); @@ -267,91 +268,95 @@ void iter_dem_instructions_include_correlations( const stim::DetectorErrorModel& detector_error_model, const Handler& handle_dem_error, std::map, std::map, double>>& joint_probabilites) { - detector_error_model.iter_flatten_error_instructions([&](const stim::DemInstruction& instruction) { - double p = instruction.arg_data[0]; - pm::DecomposedDemError decomposed_err; - decomposed_err.probability = p; - if (p > 0.5) { - throw ::std::invalid_argument( - "Errors with probability greater than 0.5 are not supported with correlations enabled"); - } - if (p == 0) { - // Ignore errors with no error probability. - return; - } - decomposed_err.components = {}; - decomposed_err.components.push_back({}); - UserEdge* component = &decomposed_err.components.back(); - // Mark component as empty to begin with. - component->node1 = SIZE_MAX; - component->node2 = SIZE_MAX; - size_t num_component_detectors = 0; - bool instruction_contains_separator = false; - for (auto& target : instruction.target_data) { - // Decompose error - if (target.is_relative_detector_id()) { - num_component_detectors++; - if (num_component_detectors == 1) { - const size_t& d1 = target.raw_id(); - component->node1 = d1; - } else if (num_component_detectors == 2) { - // Maintain invariant that node1 <= node2. - if (component->node1 <= target.raw_id()) { - component->node2 = target.raw_id(); + detector_error_model.iter_flatten_error_instructions( + [&](const stim::DemInstruction& instruction, bool include_decomposed_error_components_in_edge_weights = false) { + double p = instruction.arg_data[0]; + pm::DecomposedDemError decomposed_err; + decomposed_err.probability = p; + if (p > 0.5) { + throw ::std::invalid_argument( + "Errors with probability greater than 0.5 are not supported with correlations enabled"); + } + if (p == 0) { + // Ignore errors with no error probability. + return; + } + decomposed_err.components = {}; + decomposed_err.components.push_back({}); + UserEdge* component = &decomposed_err.components.back(); + // Mark component as empty to begin with. + component->node1 = SIZE_MAX; + component->node2 = SIZE_MAX; + size_t num_component_detectors = 0; + bool instruction_contains_separator = false; + for (auto& target : instruction.target_data) { + // Decompose error + if (target.is_relative_detector_id()) { + num_component_detectors++; + if (num_component_detectors == 1) { + const size_t& d1 = target.raw_id(); + component->node1 = d1; + } else if (num_component_detectors == 2) { + // Maintain invariant that node1 <= node2. + if (component->node1 <= target.raw_id()) { + component->node2 = target.raw_id(); + } else { + component->node2 = component->node1; + component->node1 = target.raw_id(); + } } else { - component->node2 = component->node1; - component->node1 = target.raw_id(); + // Undecomposed hyperedges are not supported + throw std::invalid_argument( + "Encountered an undecomposed error instruction with 3 or mode detectors. " + "This is not supported when using `enable_correlations=True`. " + "Did you forget to set `decompose_errors=True` when " + "converting the stim circuit to a detector error model?"); } - } else { - // Undecomposed hyperedges are not supported - throw std::invalid_argument( - "Encountered an undecomposed error instruction with 3 or mode detectors. " - "This is not supported when using `enable_correlations=True`. " - "Did you forget to set `decompose_errors=True` when " - "converting the stim circuit to a detector error model?"); + } else if (target.is_observable_id()) { + component->observable_indices.push_back(target.val()); + } else if (target.is_separator()) { + instruction_contains_separator = true; + // We cannot have num_component_detectors > 2 at this point, or we would have already thrown an + // exception + if (num_component_detectors == 0) { + throw std::invalid_argument( + "Encountered a decomposed error instruction with an undetectable component (0 detectors). " + "This is not supported."); + } + // The previous error in the decomposition must have 1 or 2 detectors + decomposed_err.components.push_back({}); + component = &decomposed_err.components.back(); + component->node1 = SIZE_MAX; + component->node2 = SIZE_MAX; + num_component_detectors = 0; } - } else if (target.is_observable_id()) { - component->observable_indices.push_back(target.val()); - } else if (target.is_separator()) { - instruction_contains_separator = true; - // We cannot have num_component_detectors > 2 at this point, or we would have already thrown an exception - if (num_component_detectors == 0) { + } + + if (num_component_detectors == 0) { + if (instruction_contains_separator) { throw std::invalid_argument( "Encountered a decomposed error instruction with an undetectable component (0 detectors). " "This is not supported."); + } else { + // Ignore errors that are undetectable, provided they are not a component of a decomposed error + return; } - // The previous error in the decomposition must have 1 or 2 detectors - decomposed_err.components.push_back({}); - component = &decomposed_err.components.back(); - component->node1 = SIZE_MAX; - component->node2 = SIZE_MAX; - num_component_detectors = 0; } - } - if (num_component_detectors == 0) { - if (instruction_contains_separator) { - throw std::invalid_argument( - "Encountered a decomposed error instruction with an undetectable component (0 detectors). " - "This is not supported."); - } else { - // Ignore errors that are undetectable, provided they are not a component of a decomposed error - return; - } - } - - // Only add the edge into the graph if it is not a component in a decomposed error with more than one component - if (decomposed_err.components.size() == 1) { - component = &decomposed_err.components.back(); - if (component->node2 == SIZE_MAX) { - handle_dem_error(p, {component->node1}, component->observable_indices); - } else { - handle_dem_error(p, {component->node1, component->node2}, component->observable_indices); + // If include_decomposed_error_components_in_edge_weights, then only add the edge into the graph if + // it is not a component in a decomposed error with more than one component + if (include_decomposed_error_components_in_edge_weights || decomposed_err.components.size() == 1) { + for (pm::UserEdge& component : decomposed_err.components) { + if (component.node2 == SIZE_MAX) { + handle_dem_error(p, {component.node1}, component.observable_indices); + } else { + handle_dem_error(p, {component.node1, component.node2}, component.observable_indices); + } + } } - } - add_decomposed_error_to_joint_probabilities(decomposed_err, joint_probabilites); - }); + add_decomposed_error_to_joint_probabilities(decomposed_err, joint_probabilites); + }); } } // namespace pm diff --git a/src/pymatching/sparse_blossom/driver/user_graph.test.cc b/src/pymatching/sparse_blossom/driver/user_graph.test.cc index 4b4f6bed3..3bdf6d81a 100644 --- a/src/pymatching/sparse_blossom/driver/user_graph.test.cc +++ b/src/pymatching/sparse_blossom/driver/user_graph.test.cc @@ -512,7 +512,7 @@ TEST(UserGraph, PopulateImpliedEdgeWeights) { ASSERT_EQ(implied_23.node1, 0); ASSERT_EQ(implied_23.node2, 1); // ASSERT_EQ(implied_23.implied_weight, 0); - ASSERT_NEAR(implied_23.implied_weight, 0.105361, 0.00001); + ASSERT_NEAR(implied_23.implied_weight, 0.0, 0.00001); } TEST(UserGraph, ConvertImpliedWeights) { From 02cb7b1bd7e96279829b4759e54dd7b70995ffc0 Mon Sep 17 00:00:00 2001 From: Oscar Higgott Date: Mon, 25 Aug 2025 08:25:43 +0000 Subject: [PATCH 05/10] fix default choice to include components, fix tests --- .../sparse_blossom/driver/user_graph.h | 156 +++++++++--------- .../sparse_blossom/driver/user_graph.test.cc | 15 +- 2 files changed, 83 insertions(+), 88 deletions(-) diff --git a/src/pymatching/sparse_blossom/driver/user_graph.h b/src/pymatching/sparse_blossom/driver/user_graph.h index 0c3b88a4c..a8c6ff8ff 100644 --- a/src/pymatching/sparse_blossom/driver/user_graph.h +++ b/src/pymatching/sparse_blossom/driver/user_graph.h @@ -267,96 +267,96 @@ template void iter_dem_instructions_include_correlations( const stim::DetectorErrorModel& detector_error_model, const Handler& handle_dem_error, - std::map, std::map, double>>& joint_probabilites) { - detector_error_model.iter_flatten_error_instructions( - [&](const stim::DemInstruction& instruction, bool include_decomposed_error_components_in_edge_weights = false) { - double p = instruction.arg_data[0]; - pm::DecomposedDemError decomposed_err; - decomposed_err.probability = p; - if (p > 0.5) { - throw ::std::invalid_argument( - "Errors with probability greater than 0.5 are not supported with correlations enabled"); - } - if (p == 0) { - // Ignore errors with no error probability. - return; - } - decomposed_err.components = {}; - decomposed_err.components.push_back({}); - UserEdge* component = &decomposed_err.components.back(); - // Mark component as empty to begin with. - component->node1 = SIZE_MAX; - component->node2 = SIZE_MAX; - size_t num_component_detectors = 0; - bool instruction_contains_separator = false; - for (auto& target : instruction.target_data) { - // Decompose error - if (target.is_relative_detector_id()) { - num_component_detectors++; - if (num_component_detectors == 1) { - const size_t& d1 = target.raw_id(); - component->node1 = d1; - } else if (num_component_detectors == 2) { - // Maintain invariant that node1 <= node2. - if (component->node1 <= target.raw_id()) { - component->node2 = target.raw_id(); - } else { - component->node2 = component->node1; - component->node1 = target.raw_id(); - } + std::map, std::map, double>>& joint_probabilites, + bool include_decomposed_error_components_in_edge_weights = true) { + detector_error_model.iter_flatten_error_instructions([&](const stim::DemInstruction& instruction) { + double p = instruction.arg_data[0]; + pm::DecomposedDemError decomposed_err; + decomposed_err.probability = p; + if (p > 0.5) { + throw ::std::invalid_argument( + "Errors with probability greater than 0.5 are not supported with correlations enabled"); + } + if (p == 0) { + // Ignore errors with no error probability. + return; + } + decomposed_err.components = {}; + decomposed_err.components.push_back({}); + UserEdge* component = &decomposed_err.components.back(); + // Mark component as empty to begin with. + component->node1 = SIZE_MAX; + component->node2 = SIZE_MAX; + size_t num_component_detectors = 0; + bool instruction_contains_separator = false; + for (auto& target : instruction.target_data) { + // Decompose error + if (target.is_relative_detector_id()) { + num_component_detectors++; + if (num_component_detectors == 1) { + const size_t& d1 = target.raw_id(); + component->node1 = d1; + } else if (num_component_detectors == 2) { + // Maintain invariant that node1 <= node2. + if (component->node1 <= target.raw_id()) { + component->node2 = target.raw_id(); } else { - // Undecomposed hyperedges are not supported - throw std::invalid_argument( - "Encountered an undecomposed error instruction with 3 or mode detectors. " - "This is not supported when using `enable_correlations=True`. " - "Did you forget to set `decompose_errors=True` when " - "converting the stim circuit to a detector error model?"); + component->node2 = component->node1; + component->node1 = target.raw_id(); } - } else if (target.is_observable_id()) { - component->observable_indices.push_back(target.val()); - } else if (target.is_separator()) { - instruction_contains_separator = true; - // We cannot have num_component_detectors > 2 at this point, or we would have already thrown an - // exception - if (num_component_detectors == 0) { - throw std::invalid_argument( - "Encountered a decomposed error instruction with an undetectable component (0 detectors). " - "This is not supported."); - } - // The previous error in the decomposition must have 1 or 2 detectors - decomposed_err.components.push_back({}); - component = &decomposed_err.components.back(); - component->node1 = SIZE_MAX; - component->node2 = SIZE_MAX; - num_component_detectors = 0; + } else { + // Undecomposed hyperedges are not supported + throw std::invalid_argument( + "Encountered an undecomposed error instruction with 3 or mode detectors. " + "This is not supported when using `enable_correlations=True`. " + "Did you forget to set `decompose_errors=True` when " + "converting the stim circuit to a detector error model?"); } - } - - if (num_component_detectors == 0) { - if (instruction_contains_separator) { + } else if (target.is_observable_id()) { + component->observable_indices.push_back(target.val()); + } else if (target.is_separator()) { + instruction_contains_separator = true; + // We cannot have num_component_detectors > 2 at this point, or we would have already thrown an + // exception + if (num_component_detectors == 0) { throw std::invalid_argument( "Encountered a decomposed error instruction with an undetectable component (0 detectors). " "This is not supported."); - } else { - // Ignore errors that are undetectable, provided they are not a component of a decomposed error - return; } + // The previous error in the decomposition must have 1 or 2 detectors + decomposed_err.components.push_back({}); + component = &decomposed_err.components.back(); + component->node1 = SIZE_MAX; + component->node2 = SIZE_MAX; + num_component_detectors = 0; } + } - // If include_decomposed_error_components_in_edge_weights, then only add the edge into the graph if - // it is not a component in a decomposed error with more than one component - if (include_decomposed_error_components_in_edge_weights || decomposed_err.components.size() == 1) { - for (pm::UserEdge& component : decomposed_err.components) { - if (component.node2 == SIZE_MAX) { - handle_dem_error(p, {component.node1}, component.observable_indices); - } else { - handle_dem_error(p, {component.node1, component.node2}, component.observable_indices); - } + if (num_component_detectors == 0) { + if (instruction_contains_separator) { + throw std::invalid_argument( + "Encountered a decomposed error instruction with an undetectable component (0 detectors). " + "This is not supported."); + } else { + // Ignore errors that are undetectable, provided they are not a component of a decomposed error + return; + } + } + + // If include_decomposed_error_components_in_edge_weights, then only add the edge into the graph if + // it is not a component in a decomposed error with more than one component + if (include_decomposed_error_components_in_edge_weights || decomposed_err.components.size() == 1) { + for (pm::UserEdge& component : decomposed_err.components) { + if (component.node2 == SIZE_MAX) { + handle_dem_error(p, {component.node1}, component.observable_indices); + } else { + handle_dem_error(p, {component.node1, component.node2}, component.observable_indices); } } + } - add_decomposed_error_to_joint_probabilities(decomposed_err, joint_probabilites); - }); + add_decomposed_error_to_joint_probabilities(decomposed_err, joint_probabilites); + }); } } // namespace pm diff --git a/src/pymatching/sparse_blossom/driver/user_graph.test.cc b/src/pymatching/sparse_blossom/driver/user_graph.test.cc index 3bdf6d81a..c38879f9d 100644 --- a/src/pymatching/sparse_blossom/driver/user_graph.test.cc +++ b/src/pymatching/sparse_blossom/driver/user_graph.test.cc @@ -320,10 +320,9 @@ TEST(IterDemInstructionsTest, DecomposedError) { pm::iter_dem_instructions_include_correlations(dem, handler, joint_probabilities); // Check handler calls - ASSERT_EQ(handler.handled_errors.size(), 0); - // ASSERT_EQ(handler.handled_errors.size(), 3); - // std::vector expected_handled = {{0.1, 0, 1, {}}, {0.1, 2, 3, {0}}, {0.1, 4, SIZE_MAX, {}}}; - // EXPECT_EQ(handler.handled_errors, expected_handled); + ASSERT_EQ(handler.handled_errors.size(), 3); + std::vector expected_handled = {{0.1, 0, 1, {}}, {0.1, 2, 3, {0}}, {0.1, 4, SIZE_MAX, {}}}; + EXPECT_EQ(handler.handled_errors, expected_handled); // Check joint probabilities std::pair key01 = {0, 1}; @@ -384,14 +383,10 @@ TEST(IterDemInstructionsTest, CombinedComplexDem) { std::map, std::map, double>> joint_probabilities; pm::iter_dem_instructions_include_correlations(dem, handler, joint_probabilities); - ASSERT_EQ(handler.handled_errors.size(), 2); + ASSERT_EQ(handler.handled_errors.size(), 4); std::vector expected = { - {0.1, 0, SIZE_MAX, {}}, {0.2, 1, 2, {0}}}; + {0.1, 0, SIZE_MAX, {}}, {0.2, 1, 2, {0}}, {0.4, 8, SIZE_MAX, {}}, {0.4, 9, SIZE_MAX, {1}}}; EXPECT_EQ(handler.handled_errors, expected); - // ASSERT_EQ(handler.handled_errors.size(), 4); - // std::vector expected = { - // {0.1, 0, SIZE_MAX, {}}, {0.2, 1, 2, {0}}, {0.4, 8, SIZE_MAX, {}}, {0.4, 9, SIZE_MAX, {1}}}; - // EXPECT_EQ(handler.handled_errors, expected); // Check joint probabilities std::pair key0B = {0, SIZE_MAX}; From 0dc38d29ee8d18dfc4e4de2ae6e3749f56b86b67 Mon Sep 17 00:00:00 2001 From: Oscar Higgott Date: Mon, 25 Aug 2025 08:50:10 +0000 Subject: [PATCH 06/10] fix test --- src/pymatching/sparse_blossom/driver/mwpm_decoding.test.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pymatching/sparse_blossom/driver/mwpm_decoding.test.cc b/src/pymatching/sparse_blossom/driver/mwpm_decoding.test.cc index a5f66f363..a40ca1fbf 100644 --- a/src/pymatching/sparse_blossom/driver/mwpm_decoding.test.cc +++ b/src/pymatching/sparse_blossom/driver/mwpm_decoding.test.cc @@ -330,7 +330,7 @@ TEST(MwpmCorrelatedDecoding, BetterLogicalErrorRateThanUncorrelated) { size_t num_mistakes_uncorrelated = 0; size_t num_shots = 0; size_t max_shots = 300; - size_t expected_mistakes_correlated = 3; + size_t expected_mistakes_correlated = 5; size_t expected_mistakes_uncorrelated = 11; pm::ExtendedMatchingResult res(mwpm_correlated.flooder.graph.num_observables); std::vector edges; From e49dcecbd6eea7f7fafc9c7c231821db12a7faa7 Mon Sep 17 00:00:00 2001 From: Oscar Higgott Date: Mon, 25 Aug 2025 08:52:20 +0000 Subject: [PATCH 07/10] format tests --- src/pymatching/sparse_blossom/driver/user_graph.test.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pymatching/sparse_blossom/driver/user_graph.test.cc b/src/pymatching/sparse_blossom/driver/user_graph.test.cc index c38879f9d..5d1e08bfe 100644 --- a/src/pymatching/sparse_blossom/driver/user_graph.test.cc +++ b/src/pymatching/sparse_blossom/driver/user_graph.test.cc @@ -384,6 +384,7 @@ TEST(IterDemInstructionsTest, CombinedComplexDem) { pm::iter_dem_instructions_include_correlations(dem, handler, joint_probabilities); ASSERT_EQ(handler.handled_errors.size(), 4); + std::vector expected = { {0.1, 0, SIZE_MAX, {}}, {0.2, 1, 2, {0}}, {0.4, 8, SIZE_MAX, {}}, {0.4, 9, SIZE_MAX, {1}}}; EXPECT_EQ(handler.handled_errors, expected); @@ -506,7 +507,6 @@ TEST(UserGraph, PopulateImpliedEdgeWeights) { const auto& implied_23 = it_23->implied_weights_for_other_edges[0]; ASSERT_EQ(implied_23.node1, 0); ASSERT_EQ(implied_23.node2, 1); - // ASSERT_EQ(implied_23.implied_weight, 0); ASSERT_NEAR(implied_23.implied_weight, 0.0, 0.00001); } From 73a25f65d6668b21762725f80cdbab4ea3ec50de Mon Sep 17 00:00:00 2001 From: Oscar Higgott Date: Thu, 25 Sep 2025 17:10:46 +0000 Subject: [PATCH 08/10] fix docstring typo --- src/pymatching/sparse_blossom/driver/user_graph.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/pymatching/sparse_blossom/driver/user_graph.h b/src/pymatching/sparse_blossom/driver/user_graph.h index a8c6ff8ff..4fedb9961 100644 --- a/src/pymatching/sparse_blossom/driver/user_graph.h +++ b/src/pymatching/sparse_blossom/driver/user_graph.h @@ -343,8 +343,8 @@ void iter_dem_instructions_include_correlations( } } - // If include_decomposed_error_components_in_edge_weights, then only add the edge into the graph if - // it is not a component in a decomposed error with more than one component + // If include_decomposed_error_components_in_edge_weights is False, then only add the edge into + // the graph if it is not a component in a decomposed error with more than one component if (include_decomposed_error_components_in_edge_weights || decomposed_err.components.size() == 1) { for (pm::UserEdge& component : decomposed_err.components) { if (component.node2 == SIZE_MAX) { From 23b3be993cf9d8c1f7063c30d8775e068282b3cd Mon Sep 17 00:00:00 2001 From: Oscar Higgott Date: Thu, 25 Sep 2025 17:51:07 +0000 Subject: [PATCH 09/10] remove support for cp38-macosx_arm64 --- .github/workflows/ci.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 7f76a5b8b..56790b93b 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -51,7 +51,7 @@ jobs: {os: macos-latest, dist: cp312-macosx_x86_64}, {os: macos-latest, dist: cp313-macosx_x86_64}, # macosx arm64 - {os: macos-latest, dist: cp38-macosx_arm64}, + # {os: macos-latest, dist: cp38-macosx_arm64}, {os: macos-latest, dist: cp39-macosx_arm64}, {os: macos-latest, dist: cp310-macosx_arm64}, {os: macos-latest, dist: cp311-macosx_arm64}, From 7fef1785cc157380b1db87e7eef985c9b29def1e Mon Sep 17 00:00:00 2001 From: Noah Shutty Date: Thu, 25 Sep 2025 10:53:07 -0700 Subject: [PATCH 10/10] Add test to check for correlated matching boundary component handling (#176) Co-authored-by: oscarhiggott <29460323+oscarhiggott@users.noreply.github.com> --- tests/matching/decode_test.py | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/tests/matching/decode_test.py b/tests/matching/decode_test.py index 6f37fe56c..3e2bcd513 100644 --- a/tests/matching/decode_test.py +++ b/tests/matching/decode_test.py @@ -400,6 +400,26 @@ def test_decode_to_edges_with_correlations(): assert np.array_equal(edges, expected_edges) +def test_correlated_matching_handles_single_detector_components(): + stim = pytest.importorskip("stim") + p = 0.1 + circuit = stim.Circuit.generated( + code_task="surface_code:rotated_memory_x", + distance=5, + rounds=5, + before_round_data_depolarization=p, + ) + circ_str = str(circuit).replace( + f"DEPOLARIZE1({p})", f"PAULI_CHANNEL_1(0, {p}, 0)" + ) + noisy_circuit = stim.Circuit(circ_str) + dem = noisy_circuit.detector_error_model( + decompose_errors=True, approximate_disjoint_errors=True + ) + m = Matching.from_detector_error_model(dem, enable_correlations=True) + assert m.num_detectors > 0 + + def test_load_from_circuit_with_correlations(): stim = pytest.importorskip("stim") circuit = stim.Circuit.generated(