Fix for #20 with --det-order-bfs, and update visualizer (#30)

Fix #20 by adding
`--det-order-bfs` argument.
Apparently it can improve performance for the same LER as well:
```
Before change:
bazel build src:all && time ./bazel-bin/src/tesseract --pqlimit 200000 --beam 15 --beam-climbing --num-det-orders 16  --sample-num-shots 10000 --det-order-seed 13267562 --circuit testdata/colorcodes/r\=9\,d\=9\,p\=0.002\,noise\=si1000\,c\=superdense_color_code_X\,q\=121\,gates\=cz.stim --sample-seed 717347  --threads 60  --print-stats
num_shots = 10000 num_low_confidence = 0 num_errors = 52 total_time_seconds = 64829.16896500005

After change:
bazel build src:all && time ./bazel-bin/src/tesseract --pqlimit 200000 --beam 5 --num-det-orders 20  --sample-num-shots 10000 --det-order-seed 13267562 --circuit testdata/colorcodes/r\=9\,d\=9\,p\=0.002\,noise\=si1000\,c\=superdense_color_code_X\,q\=121\,gates\=cz.stim --sample-seed 717347  --threads 60 --det-order-bfs
num_shots = 10000 num_low_confidence = 2 num_errors = 26 total_time_seconds = 46809.56727799983
```

Plus a few updates for the visualizer, now it should load a default
logfile automatically!

Author: noajshu

README.md

@@ -1,6 +1,6 @@
 <div align="center">
 
-# Tesseract Decoder
+# [Tesseract Decoder](https://quantumlib.github.io/tesseract-decoder)
 
 A Search-Based Decoder for Quantum Error Correction.
 
@@ -48,9 +48,10 @@ We tested the Tesseract decoder for:
 *   **Detailed Statistics:** provides comprehensive statistics output, including shot counts, error
     counts, and processing times.
 *   **Heuristics**: includes flexible heuristic options: `--beam`, `--det-penalty`,
-    `--beam-climbing`, `--no-revisit-dets`, `--at-most-two-errors-per-detector` and `--pqlimit` to
+    `--beam-climbing`, `--no-revisit-dets`, `--at-most-two-errors-per-detector`, `--det-order-bfs` and `--pqlimit` to
     improve performance while maintaining a low logical error rate. To learn more about these
     options, use `./bazel-bin/src/tesseract --help`
+*   **Visualization tool:** open the [viz directory](viz/) in your browser to view decoding results. See [viz/README.md](viz/README.md) for instructions on generating the visualization JSON.
 
 ## Installation
 

src/tesseract_main.cc

@@ -15,6 +15,9 @@
 #include <argparse/argparse.hpp>
 #include <atomic>
 #include <fstream>
+#include <queue>
+#include <algorithm>
+#include <numeric>
 #include <nlohmann/json.hpp>
 #include <thread>
 
@@ -30,6 +33,7 @@ struct Args {
   // Manifold orientation options
   uint64_t det_order_seed;
   size_t num_det_orders = 10;
+  bool det_order_bfs = false;
 
   // Sampling options
   size_t sample_num_shots = 0;
@@ -192,46 +196,88 @@ struct Args {
         }
       }
 
-      std::vector<double> inner_products(config.dem.count_detectors());
-
-      if (!detector_coords.size() or !detector_coords.at(0).size()) {
-        // If there are no detector coordinates, just use the standard ordering
-        // of the indices.
+      if (det_order_bfs) {
+        auto graph = build_detector_graph(config.dem);
+        std::uniform_int_distribution<size_t> dist_det(0, graph.size() - 1);
         for (size_t det_order = 0; det_order < num_det_orders; ++det_order) {
-          config.det_orders.emplace_back();
-          std::iota(config.det_orders.back().begin(),
-                    config.det_orders.front().end(), 0);
-        }
-      } else {
-        // Use the coordinates to order the detectors based on a random
-        // orientation
-        for (size_t det_order = 0; det_order < num_det_orders; ++det_order) {
-          // Sample a direction
-          std::vector<double> orientation_vector;
-          for (size_t i = 0; i < detector_coords.at(0).size(); ++i) {
-            orientation_vector.push_back(dist(rng));
-          }
-
-          for (size_t i = 0; i < detector_coords.size(); ++i) {
-            inner_products[i] = 0;
-            for (size_t j = 0; j < orientation_vector.size(); ++j) {
-              inner_products[i] +=
-                  detector_coords[i][j] * orientation_vector[j];
+          std::vector<size_t> perm;
+          perm.reserve(graph.size());
+          std::vector<char> visited(graph.size(), false);
+          std::queue<size_t> q;
+          size_t start = dist_det(rng);
+          while (perm.size() < graph.size()) {
+            if (!visited[start]) {
+              visited[start] = true;
+              q.push(start);
+              perm.push_back(start);
+            }
+            while (!q.empty()) {
+              size_t cur = q.front();
+              q.pop();
+              auto neigh = graph[cur];
+              std::shuffle(neigh.begin(), neigh.end(), rng);
+              for (size_t n : neigh) {
+                if (!visited[n]) {
+                  visited[n] = true;
+                  q.push(n);
+                  perm.push_back(n);
+                }
+              }
+            }
+            if (perm.size() < graph.size()) {
+              do {
+                start = dist_det(rng);
+              } while (visited[start]);
             }
           }
-          std::vector<size_t> perm(config.dem.count_detectors());
-          std::iota(perm.begin(), perm.end(), 0);
-          std::sort(perm.begin(), perm.end(),
-                    [&](const size_t& i, const size_t& j) {
-                      return inner_products[i] > inner_products[j];
-                    });
-          // Invert the permutation
-          std::vector<size_t> inv_perm(config.dem.count_detectors());
+          std::vector<size_t> inv_perm(graph.size());
           for (size_t i = 0; i < perm.size(); ++i) {
             inv_perm[perm[i]] = i;
           }
           config.det_orders[det_order] = inv_perm;
         }
+      } else {
+        std::vector<double> inner_products(config.dem.count_detectors());
+
+        if (!detector_coords.size() || !detector_coords.at(0).size()) {
+          // If there are no detector coordinates, just use the standard ordering
+          // of the indices.
+          for (size_t det_order = 0; det_order < num_det_orders; ++det_order) {
+            config.det_orders.emplace_back();
+            std::iota(config.det_orders.back().begin(),
+                      config.det_orders.front().end(), 0);
+          }
+        } else {
+          // Use the coordinates to order the detectors based on a random
+          // orientation
+          for (size_t det_order = 0; det_order < num_det_orders; ++det_order) {
+            // Sample a direction
+            std::vector<double> orientation_vector;
+            for (size_t i = 0; i < detector_coords.at(0).size(); ++i) {
+              orientation_vector.push_back(dist(rng));
+            }
+
+            for (size_t i = 0; i < detector_coords.size(); ++i) {
+              inner_products[i] = 0;
+              for (size_t j = 0; j < orientation_vector.size(); ++j) {
+                inner_products[i] +=
+                    detector_coords[i][j] * orientation_vector[j];
+              }
+            }
+            std::vector<size_t> perm(config.dem.count_detectors());
+            std::iota(perm.begin(), perm.end(), 0);
+            std::sort(perm.begin(), perm.end(),
+                      [&](const size_t& i, const size_t& j) {
+                        return inner_products[i] > inner_products[j];
+                      });
+            // Invert the permutation
+            std::vector<size_t> inv_perm(config.dem.count_detectors());
+            for (size_t i = 0; i < perm.size(); ++i) {
+              inv_perm[perm[i]] = i;
+            }
+            config.det_orders[det_order] = inv_perm;
+          }
+        }
       }
     }
 
@@ -359,6 +405,10 @@ int main(int argc, char* argv[]) {
       .metavar("N")
       .default_value(size_t(1))
       .store_into(args.num_det_orders);
+  program.add_argument("--det-order-bfs")
+      .help("Use BFS-based detector ordering instead of geometric orientation")
+      .flag()
+      .store_into(args.det_order_bfs);
   program.add_argument("--det-order-seed")
       .help(
           "Seed used when initializing the random detector traversal "

src/utils.cc

@@ -51,6 +51,36 @@ std::vector<std::vector<double>> get_detector_coords(
   return detector_coords;
 }
 
+std::vector<std::vector<size_t>> build_detector_graph(
+    const stim::DetectorErrorModel& dem) {
+  size_t num_detectors = dem.count_detectors();
+  std::vector<std::vector<size_t>> neighbors(num_detectors);
+  for (const stim::DemInstruction& instruction : dem.flattened().instructions) {
+    if (instruction.type != stim::DemInstructionType::DEM_ERROR) {
+      continue;
+    }
+    std::vector<int> dets;
+    for (const stim::DemTarget& target : instruction.target_data) {
+      if (target.is_relative_detector_id()) {
+        dets.push_back(target.val());
+      }
+    }
+    for (size_t i = 0; i < dets.size(); ++i) {
+      for (size_t j = i + 1; j < dets.size(); ++j) {
+        size_t a = dets[i];
+        size_t b = dets[j];
+        neighbors[a].push_back(b);
+        neighbors[b].push_back(a);
+      }
+    }
+  }
+  for (auto& neigh : neighbors) {
+    std::sort(neigh.begin(), neigh.end());
+    neigh.erase(std::unique(neigh.begin(), neigh.end()), neigh.end());
+  }
+  return neighbors;
+}
+
 bool sampling_from_dem(uint64_t seed, size_t num_shots,
                        stim::DetectorErrorModel dem,
                        std::vector<stim::SparseShot>& shots) {

src/utils.h

@@ -31,6 +31,11 @@ constexpr const double EPSILON = 1e-7;
 std::vector<std::vector<double>> get_detector_coords(
     stim::DetectorErrorModel& dem);
 
+// Builds an adjacency list graph where two detectors share an edge iff an error
+// in the model activates them both.
+std::vector<std::vector<size_t>> build_detector_graph(
+    const stim::DetectorErrorModel& dem);
+
 const double INF = std::numeric_limits<double>::infinity();
 
 bool sampling_from_dem(uint64_t seed, size_t num_shots,

viz/README.md

@@ -0,0 +1,34 @@
+# Visualization
+
+This tool displays the detectors and errors from a Tesseract decoding run in 3D.
+
+## Generating the JSON
+
+Use the `--verbose` flag when running the decoder and filter the output to
+include only the lines used by the converter script:
+
+```bash
+bazel build src:all && \
+./bazel-bin/src/tesseract \
+  --sample-num-shots 1 --det-order-seed 13267562 --pqlimit 10000 --beam 1 --num-det-orders 20 --det-order-bfs \
+  --circuit testdata/colorcodes/r\=9\,d\=9\,p\=0.002\,noise\=si1000\,c\=superdense_color_code_X\,q\=121\,gates\=cz.stim \
+  --sample-seed 717347 --threads 1 --verbose | \
+  grep -E 'Error|Detector|activated_errors|activated_dets' > logfile.txt
+
+python viz/to_json.py logfile.txt -o logfile.json
+```
+
+
+The `--det-order-bfs` flag is compatible with visualization logs. Just make
+sure `--verbose` is enabled so the detector coordinates are printed for
+`to_json.py` to parse.
+
+The `to_json.py` script produces `logfile.json`, which contains the detector
+coordinates and animation frames for the viewer.
+
+## Viewing
+
+Open `viz/index.html` in a modern browser. It will automatically try to load
+`logfile.json` from the same directory. If the file picker is used, any JSON
+produced by `to_json.py` can be visualized.
+

viz/index.html

@@ -314,6 +314,18 @@
   if (playing && frames.length) restartAnimationTimer();
   else clearInterval(animationTimer);
 }
+
+// Try to load a default logfile when the page is opened.
+window.addEventListener('load', () => {
+  fetch('logfile.json')
+    .then(r => r.ok ? r.json() : Promise.reject())
+    .then(data => initializeScene(
+      data.detectorCoords,
+      data.errorCoords,
+      data.frames,
+      data.errorToDetectors))
+    .catch(() => {/* ignore if not present */});
+});
 </script>
 </body>
 </html>