Heuristic Benchmarking (#278)

docs/source/library/MappingResults.rst

@@ -1,9 +1,10 @@
 Results
 =======
 
-This class captures additional results from the :code:`compile` method.
-
     .. currentmodule:: mqt.qmap
+
+This class captures additional results from the :func:`compile` method.
+
     .. autoclass:: MappingResults
 
 In addition to the mapped circuit
@@ -20,10 +21,18 @@ and information about the input and output circuits.
     .. autoattribute:: MappingResults.input
     .. autoattribute:: MappingResults.output
 
-If specified in the call to the :code:`compile` method, the `weighted MaxSAT formula <http://www.maxhs.org/docs/wdimacs.html>`_ is also tracked.
+If specified in the call to the :func:`compile` method, the `weighted MaxSAT formula <http://www.maxhs.org/docs/wdimacs.html>`_ is also tracked.
 
     .. autoattribute:: MappingResults.wcnf
 
+If the heuristic mapper is used and :code:`debug` is set to :code:`True`, some benchmark information is provided for the whole circuit
+
+    .. autoattribute:: MappingResults.heuristic_benchmark
+
+and for each layer separately.
+
+    .. autoattribute:: MappingResults.layer_heuristic_benchmark
+
 In addition, the class provides methods to export to other formats.
 
     .. automethod:: MappingResults.json

include/MappingResults.hpp

@@ -28,6 +28,15 @@ struct MappingResults {
     std::size_t teleportations   = 0;
   };
 
+  struct HeuristicBenchmarkInfo {
+    std::size_t expandedNodes            = 0;
+    std::size_t generatedNodes           = 0;
+    std::size_t solutionDepth            = 0;
+    double      timePerNode              = 0.;
+    double      averageBranchingFactor   = 0.;
+    double      effectiveBranchingFactor = 0.;
+  };
+
   CircuitInfo input{};
 
   std::string   architecture{};
@@ -41,15 +50,20 @@ struct MappingResults {
 
   std::string wcnf{};
 
+  HeuristicBenchmarkInfo              heuristicBenchmark{};
+  std::vector<HeuristicBenchmarkInfo> layerHeuristicBenchmark{};
+
   MappingResults()          = default;
   virtual ~MappingResults() = default;
 
   virtual void copyInput(const MappingResults& mappingResults) {
-    input        = mappingResults.input;
-    architecture = mappingResults.architecture;
-    config       = mappingResults.config;
-    output       = mappingResults.output;
-    wcnf         = mappingResults.wcnf;
+    input                   = mappingResults.input;
+    architecture            = mappingResults.architecture;
+    config                  = mappingResults.config;
+    output                  = mappingResults.output;
+    wcnf                    = mappingResults.wcnf;
+    heuristicBenchmark      = mappingResults.heuristicBenchmark;
+    layerHeuristicBenchmark = mappingResults.layerHeuristicBenchmark;
   }
 
   [[nodiscard]] std::string toString() const { return json().dump(2); }
@@ -88,7 +102,15 @@ struct MappingResults {
         stats["WCNF"] = wcnf;
       }
     } else if (config.method == Method::Heuristic) {
-      stats["teleportations"] = output.teleportations;
+      stats["teleportations"]      = output.teleportations;
+      auto& benchmark              = stats["benchmark"];
+      benchmark["expanded_nodes"]  = heuristicBenchmark.expandedNodes;
+      benchmark["generated_nodes"] = heuristicBenchmark.generatedNodes;
+      benchmark["time_per_node"]   = heuristicBenchmark.timePerNode;
+      benchmark["average_branching_factor"] =
+          heuristicBenchmark.averageBranchingFactor;
+      benchmark["effective_branching_factor"] =
+          heuristicBenchmark.effectiveBranchingFactor;
     }
     stats["additional_gates"] =
         static_cast<std::make_signed_t<decltype(output.gates)>>(output.gates) -

include/configuration/Configuration.hpp

@@ -27,6 +27,7 @@ struct Configuration {
   bool addMeasurementsToMappedCircuit = true;
 
   bool verbose = false;
+  bool debug   = false;
 
   // map to particular subgraph of architecture (in exact mapper)
   std::set<std::uint16_t> subgraph{};

include/heuristic/HeuristicMapper.hpp

@@ -6,6 +6,8 @@
 #include "Mapper.hpp"
 #include "heuristic/UniquePriorityQueue.hpp"
 
+#include <cmath>
+
 #pragma once
 
 /**
@@ -25,6 +27,8 @@ class HeuristicMapper : public Mapper {
 public:
   using Mapper::Mapper; // import constructors from parent class
 
+  static constexpr double EFFECTIVE_BRANCH_RATE_TOLERANCE = 1e-10;
+
   /**
    * @brief map the circuit passed at initialization to the architecture
    *
@@ -70,13 +74,15 @@ class HeuristicMapper : public Mapper {
     /** number of swaps used to get from mapping after last layer to the current
      * mapping */
     std::size_t nswaps = 0;
+    /** depth in search tree (starting with 0 at the root) */
+    std::size_t depth = 0;
 
     Node() = default;
     Node(const std::array<std::int16_t, MAX_DEVICE_QUBITS>& q,
          const std::array<std::int16_t, MAX_DEVICE_QUBITS>& loc,
-         const std::vector<std::vector<Exchange>>&          sw            = {},
-         const double                                       initCostFixed = 0)
-        : costFixed(initCostFixed) {
+         const std::vector<std::vector<Exchange>>&          sw = {},
+         const double initCostFixed = 0, const std::size_t searchDepth = 0)
+        : costFixed(initCostFixed), depth(searchDepth) {
       std::copy(q.begin(), q.end(), qubits.begin());
       std::copy(loc.begin(), loc.end(), locations.begin());
       std::copy(sw.begin(), sw.end(), std::back_inserter(swaps));
@@ -263,6 +269,32 @@ class HeuristicMapper : public Mapper {
     }
     return currentCost + newCost;
   }
+
+  static double computeEffectiveBranchingRate(std::size_t       nodesProcessed,
+                                              const std::size_t solutionDepth) {
+    // N = (b*)^d + (b*)^(d-1) + ... + (b*)^2 + b* + 1
+    // no closed-form solution for b*, so we use approximation via binary search
+    if (solutionDepth == 0) {
+      return 0.;
+    }
+    --nodesProcessed; // N - 1 = (b*)^d + (b*)^(d-1) + ... + (b*)^2 + b*
+    double upper = std::pow(static_cast<double>(nodesProcessed),
+                            1.0 / static_cast<double>(solutionDepth));
+    double lower = upper / static_cast<double>(solutionDepth);
+    while (upper - lower > 2 * EFFECTIVE_BRANCH_RATE_TOLERANCE) {
+      const double mid = (lower + upper) / 2.0;
+      double       sum = 0.0;
+      for (std::size_t i = 1; i <= solutionDepth; ++i) {
+        sum += std::pow(mid, i);
+      }
+      if (sum < static_cast<double>(nodesProcessed)) {
+        lower = mid;
+      } else {
+        upper = mid;
+      }
+    }
+    return (lower + upper) / 2.0;
+  }
 };
 
 inline bool operator<(const HeuristicMapper::Node& x,

mqt/qmap/bindings.cpp

@@ -182,6 +182,7 @@ PYBIND11_MODULE(pyqmap, m) {
       .def(py::init<>())
       .def_readwrite("method", &Configuration::method)
       .def_readwrite("verbose", &Configuration::verbose)
+      .def_readwrite("debug", &Configuration::debug)
       .def_readwrite("layering", &Configuration::layering)
       .def_readwrite("initial_layout", &Configuration::initialLayout)
       .def_readwrite("lookahead", &Configuration::lookahead)
@@ -227,6 +228,9 @@ PYBIND11_MODULE(pyqmap, m) {
       .def_readwrite("time", &MappingResults::time)
       .def_readwrite("timeout", &MappingResults::timeout)
       .def_readwrite("mapped_circuit", &MappingResults::mappedCircuit)
+      .def_readwrite("heuristic_benchmark", &MappingResults::heuristicBenchmark)
+      .def_readwrite("layer_heuristic_benchmark",
+                     &MappingResults::layerHeuristicBenchmark)
       .def_readwrite("wcnf", &MappingResults::wcnf)
       .def("json", &MappingResults::json)
       .def("csv", &MappingResults::csv)
@@ -249,6 +253,25 @@ PYBIND11_MODULE(pyqmap, m) {
       .def_readwrite("teleportations",
                      &MappingResults::CircuitInfo::teleportations);
 
+  // Heuristic benchmark information
+  py::class_<MappingResults::HeuristicBenchmarkInfo>(
+      m, "HeuristicBenchmarkInfo", "Heuristic benchmark information")
+      .def(py::init<>())
+      .def_readwrite("expanded_nodes",
+                     &MappingResults::HeuristicBenchmarkInfo::expandedNodes)
+      .def_readwrite("generated_nodes",
+                     &MappingResults::HeuristicBenchmarkInfo::generatedNodes)
+      .def_readwrite("solution_depth",
+                     &MappingResults::HeuristicBenchmarkInfo::solutionDepth)
+      .def_readwrite("time_per_node",
+                     &MappingResults::HeuristicBenchmarkInfo::timePerNode)
+      .def_readwrite(
+          "average_branching_factor",
+          &MappingResults::HeuristicBenchmarkInfo::averageBranchingFactor)
+      .def_readwrite(
+          "effective_branching_factor",
+          &MappingResults::HeuristicBenchmarkInfo::effectiveBranchingFactor);
+
   auto arch = py::class_<Architecture>(
       m, "Architecture", "Class representing device/backend information");
   auto properties = py::class_<Architecture::Properties>(

mqt/qmap/compile.py

@@ -76,6 +76,7 @@ def compile(  # noqa: A001
     post_mapping_optimizations: bool = True,
     add_measurements_to_mapped_circuit: bool = True,
     verbose: bool = False,
+    debug: bool = False,
 ) -> tuple[QuantumCircuit, MappingResults]:
     """Interface to the MQT QMAP tool for mapping quantum circuits.
 
@@ -101,6 +102,7 @@ def compile(  # noqa: A001
         post_mapping_optimizations: Run post-mapping optimizations. Defaults to True.
         add_measurements_to_mapped_circuit: Whether to add measurements at the end of the mapped circuit. Defaults to True.
         verbose: Print more detailed information during the mapping process. Defaults to False.
+        debug: Gather additional information during the mapping process (e.g. number of generated nodes, branching factors, ...). Defaults to False.
 
     Returns:
         The mapped circuit and the mapping results.
@@ -134,6 +136,7 @@ def compile(  # noqa: A001
     config.post_mapping_optimizations = post_mapping_optimizations
     config.add_measurements_to_mapped_circuit = add_measurements_to_mapped_circuit
     config.verbose = verbose
+    config.debug = debug
 
     results = map(circ, architecture, config)
 

mqt/qmap/pyqmap.pyi

@@ -137,6 +137,7 @@ class Configuration:
     use_subsets: bool
     use_teleportation: bool
     verbose: bool
+    debug: bool
     def __init__(self) -> None: ...
     def json(self) -> dict[str, Any]: ...
 

src/configuration/Configuration.cpp

@@ -16,6 +16,7 @@ nlohmann::json Configuration::json() const {
   config["post_mapping_optimizations"]         = postMappingOptimizations;
   config["add_measurements_to_mapped_circuit"] = addMeasurementsToMappedCircuit;
   config["verbose"]                            = verbose;
+  config["debug"]                              = debug;
 
   if (method == Method::Heuristic) {
     auto& heuristic             = config["settings"];

src/heuristic/HeuristicMapper.cpp

@@ -8,6 +8,7 @@
 #include <chrono>
 
 void HeuristicMapper::map(const Configuration& configuration) {
+  results        = MappingResults{};
   results.config = configuration;
   auto& config   = results.config;
   if (config.layering == Layering::OddGates ||
@@ -138,6 +139,24 @@ void HeuristicMapper::map(const Configuration& configuration) {
     }
   }
 
+  if (config.debug && results.heuristicBenchmark.expandedNodes > 0) {
+    auto& benchmark = results.heuristicBenchmark;
+    benchmark.timePerNode /= static_cast<double>(benchmark.expandedNodes);
+    benchmark.averageBranchingFactor =
+        static_cast<double>(benchmark.generatedNodes - layers.size()) /
+        static_cast<double>(benchmark.expandedNodes);
+    for (const auto& layer : results.layerHeuristicBenchmark) {
+      benchmark.effectiveBranchingFactor +=
+          layer.effectiveBranchingFactor *
+          (static_cast<double>(layer.expandedNodes) /
+           static_cast<double>(benchmark.expandedNodes));
+    }
+    if (benchmark.effectiveBranchingFactor > benchmark.averageBranchingFactor) {
+      throw QMAPException("Something wrong in benchmark tracking: "
+                          "effectiveBranchingFactor > averageBranchingFactor");
+    }
+  }
+
   // infer output permutation from qubit locations
   qcMapped.outputPermutation.clear();
   std::size_t count = 0U;
@@ -439,14 +458,57 @@ HeuristicMapper::Node HeuristicMapper::aStarMap(size_t layer) {
 
   nodes.push(node);
 
+  const auto& debug = results.config.debug;
+  const auto  start = std::chrono::steady_clock::now();
+  if (debug) {
+    results.layerHeuristicBenchmark.emplace_back();
+  }
+  auto& totalExpandedNodes = results.heuristicBenchmark.expandedNodes;
+  auto  layerResultsIt     = results.layerHeuristicBenchmark.rbegin();
+
   while (!nodes.top().done) {
     Node current = nodes.top();
     nodes.pop();
     expandNode(consideredQubits, current, layer, twoQubitGateMultiplicity);
+
+    if (debug) {
+      ++totalExpandedNodes;
+      ++layerResultsIt->expandedNodes;
+    }
   }
 
   Node result = nodes.top();
-  nodes.pop();
+  if (debug) {
+    const auto end = std::chrono::steady_clock::now();
+
+    layerResultsIt->solutionDepth = result.depth;
+
+    const std::chrono::duration<double> diff = end - start;
+    results.heuristicBenchmark.timePerNode += diff.count();
+
+    layerResultsIt->generatedNodes =
+        layerResultsIt->expandedNodes + nodes.size();
+    results.heuristicBenchmark.generatedNodes += layerResultsIt->generatedNodes;
+
+    if (layerResultsIt->expandedNodes > 0) {
+      layerResultsIt->timePerNode =
+          diff.count() / static_cast<double>(layerResultsIt->expandedNodes);
+      layerResultsIt->averageBranchingFactor =
+          static_cast<double>(layerResultsIt->generatedNodes - 1) /
+          static_cast<double>(layerResultsIt->expandedNodes);
+    }
+
+    layerResultsIt->effectiveBranchingFactor = computeEffectiveBranchingRate(
+        layerResultsIt->expandedNodes + 1, result.depth);
+
+    if (layerResultsIt->effectiveBranchingFactor >
+        layerResultsIt->averageBranchingFactor) {
+      throw QMAPException(
+          "Something wrong in benchmark tracking on layer " +
+          std::to_string(layer) +
+          ": effectiveBranchingFactor > averageBranchingFactor");
+    }
+  }
 
   // clear nodes
   while (!nodes.empty()) {
@@ -536,7 +598,8 @@ void HeuristicMapper::expandNodeAddOneSwap(
     const TwoQubitMultiplicity& twoQubitGateMultiplicity) {
   const auto& config = results.config;
 
-  Node newNode = Node(node.qubits, node.locations, node.swaps, node.costFixed);
+  Node newNode = Node(node.qubits, node.locations, node.swaps, node.costFixed,
+                      node.depth + 1);
 
   if (architecture.getCouplingMap().find(swap) !=
           architecture.getCouplingMap().end() ||