fix: resolve C++ detector-to-component correlation analyzer bug

TAG=agy

Author: aria-googler

src/error_correlations.cc

@@ -1,123 +1,123 @@
 #include "error_correlations.h"
-#include <sstream>
+
 #include <iostream>
+#include <sstream>
 
 namespace tesseract {
 
 std::string ImpliedProbability::str() const {
-    std::stringstream ss;
-    ss << "ImpliedProbability(affected={";
-    for (size_t i = 0; i < affected_hyperedge.size(); ++i) {
-        ss << affected_hyperedge[i] << (i == affected_hyperedge.size() - 1 ? "" : ",");
-    }
-    ss << "}, prob=" << probability << ")";
-    return ss.str();
+  std::stringstream ss;
+  ss << "ImpliedProbability(affected={";
+  for (size_t i = 0; i < affected_hyperedge.size(); ++i) {
+    ss << affected_hyperedge[i] << (i == affected_hyperedge.size() - 1 ? "" : ",");
+  }
+  ss << "}, prob=" << probability << ")";
+  return ss.str();
 }
 
 bool ImpliedProbability::operator==(const ImpliedProbability& other) const {
-    return affected_hyperedge == other.affected_hyperedge && 
-           std::abs(probability - other.probability) < 1e-12;
+  return affected_hyperedge == other.affected_hyperedge &&
+         std::abs(probability - other.probability) < 1e-12;
 }
 
 bool ImpliedProbability::operator<(const ImpliedProbability& other) const {
-    if (affected_hyperedge != other.affected_hyperedge) {
-        return affected_hyperedge < other.affected_hyperedge;
-    }
-    return probability < other.probability;
+  if (affected_hyperedge != other.affected_hyperedge) {
+    return affected_hyperedge < other.affected_hyperedge;
+  }
+  return probability < other.probability;
 }
 
-JointProbsMap get_hyperedge_joint_probabilities(const stim::DetectorErrorModel& dem) {
-    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<Hyperedge> components;
-        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()) {
-                Hyperedge hyperedge;
-                for (size_t j = group_start; j < k; ++j) {
-                    const auto& target = inst.target_data[j];
-                    if (target.is_relative_detector_id()) {
-                        hyperedge.push_back(target.val());
-                    }
-                }
-                if (!hyperedge.empty()) {
-                    std::sort(hyperedge.begin(), hyperedge.end());
-                    components.push_back(hyperedge);
-                }
-                group_start = k + 1;
-            }
-        }
+JointProbsMap get_hyperedge_joint_probabilities(const stim::DetectorErrorModel& dem,
+                                                const std::vector<int>& global_det_to_comp_id) {
+  JointProbsMap joint_probs;
+  auto flattened = dem.flattened();
 
-        // 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]);
-        }
+  for (const auto& inst : flattened.instructions) {
+    if (inst.type != stim::DemInstructionType::DEM_ERROR) continue;
+
+    double p = inst.arg_data[0];
+
+    std::map<int, Hyperedge> comp_targets;
+    for (const auto& target : inst.target_data) {
+      if (target.is_relative_detector_id()) {
+        int d = target.val();
+        int cid =
+            (d >= 0 && (size_t)d < global_det_to_comp_id.size()) ? global_det_to_comp_id[d] : -1;
+        if (cid != -1) comp_targets[cid].push_back(d);
+      }
+    }
+
+    std::vector<Hyperedge> components;
+    for (auto& [cid, h] : comp_targets) {
+      std::sort(h.begin(), h.end());
+      components.push_back(h);
+    }
+
+    // 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]);
-                }
-            }
+    // 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;
+  return joint_probs;
 }
 
 ImpliedProbsMap get_implied_hyperedge_probabilities(const JointProbsMap& joint_probs) {
-    ImpliedProbsMap implied_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;
-        }
+  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;
+    if (p_causal <= 0 || p_causal >= 1.0) continue;
 
-        for (const auto& [affected, p_joint] : affected_map) {
-            if (causal == affected) 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;
+      // Conditional Probability P(affected | causal) = P(affected and causal) / P(causal)
+      double p_conditional = p_joint / p_causal;
 
-            implied_probs[causal].push_back({affected, p_conditional});
-        }
+      // 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;
+  return implied_probs;
 }
 
-ImpliedProbsMap process_dem_correlations(const stim::DetectorErrorModel& dem) {
-    auto joint = get_hyperedge_joint_probabilities(dem);
-    return get_implied_hyperedge_probabilities(joint);
+ImpliedProbsMap process_dem_correlations(const stim::DetectorErrorModel& dem,
+                                         const std::vector<int>& 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
+}  // namespace tesseract

src/error_correlations.h

@@ -1,12 +1,12 @@
 #ifndef ERROR_CORRELATIONS_H
 #define ERROR_CORRELATIONS_H
 
+#include <algorithm>
 #include <cmath>
 #include <map>
 #include <numeric>
-#include <vector>
 #include <string>
-#include <algorithm>
+#include <vector>
 
 #include "stim.h"
 
@@ -16,27 +16,29 @@ namespace tesseract {
  * Represents a probability adjustment for an affected hyperedge given a causal hyperedge.
  */
 struct ImpliedProbability {
-    std::vector<int> affected_hyperedge;
-    double probability; // Represents the conditional probability P(affected | causal)
+  std::vector<int> affected_hyperedge;
+  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;
+  std::string str() const;
+  bool operator==(const ImpliedProbability& other) const;
+  bool operator<(const ImpliedProbability& other) const;
 };
 
 // Type alias for hyperedge (sorted detector indices)
 using Hyperedge = std::vector<int>;
 // Type alias for joint probabilities map: causal_hyperedge -> {affected_hyperedge -> joint_prob}
 using JointProbsMap = std::map<Hyperedge, std::map<Hyperedge, double>>;
-// Type alias for implied probabilities map: causal_hyperedge -> list of conditional probability updates
+// Type alias for implied probabilities map: causal_hyperedge -> list of conditional probability
+// updates
 using ImpliedProbsMap = std::map<Hyperedge, std::vector<ImpliedProbability>>;
 
 /**
  * Calculates marginal and joint probabilities for hyperedges in a DEM.
- * Note: Assumes the input DEM has NOT been decomposed yet, as we need bridging errors 
+ * Note: Assumes the input DEM has NOT been decomposed yet, as we need bridging errors
  * to find joint probabilities.
  */
-JointProbsMap get_hyperedge_joint_probabilities(const stim::DetectorErrorModel& dem);
+JointProbsMap get_hyperedge_joint_probabilities(const stim::DetectorErrorModel& dem,
+                                                const std::vector<int>& global_det_to_comp_id);
 
 /**
  * Calculates conditional probabilities from joint probabilities.
@@ -46,8 +48,9 @@ ImpliedProbsMap get_implied_hyperedge_probabilities(const JointProbsMap& joint_p
 /**
  * Complete workflow for analyzing correlations within a stim::DetectorErrorModel.
  */
-ImpliedProbsMap process_dem_correlations(const stim::DetectorErrorModel& dem);
+ImpliedProbsMap process_dem_correlations(const stim::DetectorErrorModel& dem,
+                                         const std::vector<int>& global_det_to_comp_id);
 
-} // namespace tesseract
+}  // namespace tesseract
 
-#endif // ERROR_CORRELATIONS_H
+#endif  // ERROR_CORRELATIONS_H

src/multi_pass_tesseract_decoder.cc

@@ -53,7 +53,6 @@ void MultiPassTesseractDecoder::initialize(const stim::DetectorErrorModel& dem,
   // std::cout << "DEBUG decomposed:\n" << decomposed << std::endl;
   stim::DetectorErrorModel merged = merge_indistinguishable_errors(decomposed);
   // std::cout << "DEBUG merged:\n" << merged << std::endl;
-  ImpliedProbsMap raw_correlations = process_dem_correlations(merged);
 
   std::set<int> unique_classes;
   for (int c : detector_classes)
@@ -82,6 +81,8 @@ void MultiPassTesseractDecoder::initialize(const stim::DetectorErrorModel& dem,
     }
   }
 
+  ImpliedProbsMap raw_correlations = process_dem_correlations(flattened, global_det_to_comp_id);
+
   auto component_dems_raw = split_dem_by_component(merged, [&](int d) {
     return (d >= 0 && (size_t)d < total_global_detectors) ? global_det_to_comp_id[d] : -1;
   });

src/multi_pass_tesseract_decoder.test.cc

@@ -269,3 +269,29 @@ TEST(MultiPassTesseractDecoderTest, BoundaryConditionAndCappingTest) {
     std::vector<uint64_t> hits = {0};
     ASSERT_NO_THROW(decoder.decode(hits));
 }
+
+TEST(MultiPassTesseractDecoderTest, IntermediatePassLeakageTest) {
+    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<double>& coords, const std::string& tag) -> int {
+        return index;
+    };
+
+    TesseractConfig config;
+    config.dem = dem;
+
+    MultiPassTesseractDecoder decoder(dem, 3, classifier, config, 1, DetOrder::DetIndex, 12345, SchedulingStrategy::Causal);
+
+    std::vector<uint64_t> 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);
+}