Improvements for the Paper "Depth-Optimal Synthesis of Clifford Circuits with SAT Solvers" (#304)

## Description

This PR introduces a couple of changes simultaneously that were
introduced during the work culminating in [this
paper](https://arxiv.org/abs/2305.01674). These are:

- A parallel heuristic based on vertical circuit decomposition.
- Various improvements to the symmetry breaking constraints
- Allowing for setting of solver parameters from the Python side
- Some smaller fixes and improvements

Fixes # (issue)

## Checklist:

<!---
This checklist serves as a reminder of a couple of things that ensure
your pull request will be merged swiftly.
-->

- [x] The pull request only contains commits that are related to it.
- [x] I have added appropriate tests and documentation.
- [x] I have made sure that all CI jobs on GitHub pass.
- [x] The pull request introduces no new warnings and follows the
project's style guidelines.

Author: burgholzer

docs/source/Publications.rst

@@ -4,12 +4,12 @@ Publications
 *QMAP* is academic software. Thus, many of its built-in algorithms have been published as scientific papers.
 See :cite:labelpar:`wille2023qmap` for a general overview of *QMAP* and its features.
 
-If you use *QMAP* in your work, we would appreciate if you cited :cite:labelpar:`zulehnerEfficientMethodologyMapping2019` when using the heuristic mapper, :cite:labelpar:`willeMappingQuantumCircuits2019` when using the exact mapper, and :cite:labelpar:`schneider2023satEncodingOptimalClifford` when using the Clifford circuit synthesis approach.
+If you use *QMAP* in your work, we would appreciate if you cited :cite:labelpar:`zulehnerEfficientMethodologyMapping2019` when using the heuristic mapper, :cite:labelpar:`willeMappingQuantumCircuits2019` when using the exact mapper, and :cite:labelpar:`schneider2023satEncodingOptimalClifford`, :cite:labelpar:`peham2023DepthOptimalSynthesis` when using the Clifford circuit synthesis approach.
 Furthermore, if you use any of the particular algorithms such as
 
 - the heuristic mapping scheme using teleportation :cite:labelpar:`hillmichExlpoitingQuantumTeleportation2021`
 - the search space limitation techniques of the exact mapper (some of which are enabled per default) :cite:labelpar:`burgholzer2022limitingSearchSpace`
-- the method for finding (near-)optimal subarchitectures :cite:labelpar:`peham2022OptimalSubarchitectures`
+- the method for finding (near-)optimal subarchitectures :cite:labelpar:`peham2023OptimalSubarchitectures`
 
 please consider citing their respective papers as well. A full list of related papers is given below.
 

docs/source/Synthesis.ipynb

@@ -191,7 +191,75 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "The `include_destabilizers` flag guarantees that the unitary of the circuit is preserved during optimization."
+    "The `include_destabilizers` flag guarantees that the unitary of the circuit is preserved during optimization.\n",
+    "\n",
+    "By default *QMAP* generates optimal Clifford circuits with respect to the target metric. This might lead to runtime problems when trying to optimize larger circuits. When optimizing for depth, *QMAP* provides a heuristic that splits the circuits into several independent parts and optimizes them separately. This allows to optimize larger circuits while not guaranteeing that the depth-optimal circuit is synthesized.\n",
+    "\n",
+    "The heuristic synthesizer can be used as follows:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from qiskit import QuantumCircuit\n",
+    "from mqt import qmap\n",
+    "\n",
+    "qc = QuantumCircuit(2)\n",
+    "qc.x(0)\n",
+    "qc.cx(0, 1)\n",
+    "qc.x(0)\n",
+    "qc.s(1)\n",
+    "qc.x(1)\n",
+    "qc.cx(1, 0)\n",
+    "qc.x(1)\n",
+    "\n",
+    "qc_opt, results = qmap.optimize_clifford(\n",
+    "    circuit=qc, heuristic=True, split_size=3, include_destabilizers=True, target_metric=\"depth\"\n",
+    ")\n",
+    "\n",
+    "qc_opt.draw(output=\"mpl\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The parameter `split_size` determines how many layers of the circuit are optimized individually.\n",
+    "\n",
+    "In this example the synthesized circuit does not have optimal depth as can be checked by running the optimal synthesis method:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from qiskit import QuantumCircuit\n",
+    "from mqt import qmap\n",
+    "\n",
+    "qc = QuantumCircuit(2)\n",
+    "qc.x(0)\n",
+    "qc.cx(0, 1)\n",
+    "qc.x(0)\n",
+    "qc.s(1)\n",
+    "qc.x(1)\n",
+    "qc.cx(1, 0)\n",
+    "qc.x(1)\n",
+    "\n",
+    "qc_opt, results = qmap.optimize_clifford(circuit=qc, heuristic=False, include_destabilizers=True, target_metric=\"depth\")\n",
+    "\n",
+    "qc_opt.draw(output=\"mpl\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "However, the heuristic still gives a good depth reduction in many cases."
    ]
   },
   {

docs/source/library/Subarchitectures.rst

@@ -1,7 +1,7 @@
 Optimal Subarchitectures
 ========================
 
-To compute (near-)optimal subarchitectures of quantum computing architectures with restricted connectivity as described in :cite:labelpar:`peham2022OptimalSubarchitectures` the :code:`SubarchitectureOrder` class is provided. This class has functionality to compute the quasi-order that allows for fast computation of optimal subarchitectures.
+To compute (near-)optimal subarchitectures of quantum computing architectures with restricted connectivity as described in :cite:labelpar:`peham2023OptimalSubarchitectures` the :code:`SubarchitectureOrder` class is provided. This class has functionality to compute the quasi-order that allows for fast computation of optimal subarchitectures.
 
 Note that the initial construction of the ordering might take a while for larger architectures.
 

docs/source/refs.bib

@@ -51,12 +51,16 @@ @inproceedings{boteaComplexityQuantumCircuit2018
   year = {2018},
 }
 
-@inproceedings{peham2022OptimalSubarchitectures,
-  title = {On Optimal Subarchitectures for Quantum Circuit Mapping},
-  author = {Peham, Tom and Burgholzer, Lukas and Wille, Robert},
-  booktitle = {arXiv:2210.09321},
-  year = {2022},
-  url = {https://arxiv.org/pdf/2210.09321.pdf},
+@article{peham2023OptimalSubarchitectures,
+author = {Peham, Tom and Burgholzer, Lukas and Wille, Robert},
+title = {On Optimal Subarchitectures for Quantum Circuit Mapping},
+year = {2023},
+publisher = {Association for Computing Machinery},
+address = {New York, NY, USA},
+issn = {2643-6809},
+url = {https://doi.org/10.1145/3593594},
+doi = {10.1145/3593594},
+journal = {ACM Transactions on Quantum Computing},
 }
 
 @INPROCEEDINGS{schneider2023satEncodingOptimalClifford,
@@ -74,3 +78,11 @@ @inproceedings{wille2023qmap
   year = {2023},
   url = {https://www.cda.cit.tum.de/files/eda/2023_ispd_mqt_qmap_efficient_quantum_circuit_mapping.pdf}
 }
+
+@inproceedings{peham2023DepthOptimalSynthesis,
+  title = {Depth-Optimal Synthesis of Clifford Circuits with SAT Solvers},
+  author = {Peham, Tom and Burgholzer, Lukas and Wille, Robert},
+  booktitle = {arXiv:2305.01674},
+  year = {2023},
+  url = {https://arxiv.org/pdf/2305.01674.pdf},
+}

include/cliffordsynthesis/CliffordSynthesizer.hpp

@@ -31,24 +31,31 @@ class CliffordSynthesizer {
       : initialTableau(std::move(initial)),
         targetTableau(qc, 0, std::numeric_limits<std::size_t>::max(),
                       initialTableau.hasDestabilizers()),
+        initialCircuit(std::make_shared<qc::QuantumComputation>(qc.clone())),
         results(qc, targetTableau) {}
   explicit CliffordSynthesizer(qc::QuantumComputation& qc,
                                const bool              useDestabilizers = false)
       : initialTableau(qc.getNqubits(), useDestabilizers),
         targetTableau(qc, 0, std::numeric_limits<std::size_t>::max(),
                       useDestabilizers),
+        initialCircuit(std::make_shared<qc::QuantumComputation>(qc.clone())),
         results(qc, targetTableau) {}
 
   virtual ~CliffordSynthesizer() = default;
 
   void synthesize(const Configuration& config = {});
 
-  [[nodiscard]] Results&                getResults() { return results; };
-  [[nodiscard]] qc::QuantumComputation& getResultCircuit() {
+  [[nodiscard]] Results& getResults() { return results; };
+
+  void initResultCircuitFromResults() {
     std::stringstream ss;
     ss << results.getResultCircuit();
     resultCircuit = std::make_unique<qc::QuantumComputation>();
     resultCircuit->import(ss, qc::Format::OpenQASM);
+  }
+
+  [[nodiscard]] qc::QuantumComputation& getResultCircuit() {
+    initResultCircuitFromResults();
     return *resultCircuit;
   };
   [[nodiscard]] Tableau& getResultTableau() {
@@ -59,8 +66,9 @@ class CliffordSynthesizer {
   }
 
 protected:
-  Tableau initialTableau{};
-  Tableau targetTableau{};
+  Tableau                                 initialTableau{};
+  Tableau                                 targetTableau{};
+  std::shared_ptr<qc::QuantumComputation> initialCircuit{};
 
   Configuration configuration{};
 
@@ -84,6 +92,7 @@ class CliffordSynthesizer {
                             std::size_t upper);
   void depthOptimalSynthesis(EncoderConfig config, std::size_t lower,
                              std::size_t upper);
+  void depthHeuristicSynthesis();
   void twoQubitGateOptimalSynthesis(EncoderConfig config, std::size_t lower,
                                     std::size_t upper);
 
@@ -114,6 +123,33 @@ class CliffordSynthesizer {
     INFO() << "Found optimum: " << lowerBound;
   }
 
+  template <typename T>
+  void runLinearSearch(T& value, T lowerBound, T upperBound,
+                       const EncoderConfig& config) {
+    INFO() << "Running linear search in range [" << lowerBound << ", "
+           << upperBound << ")";
+
+    if (upperBound == 0U) {
+      upperBound = std::numeric_limits<std::size_t>::max();
+    }
+    for (value = lowerBound; value < upperBound; ++value) {
+      INFO() << "Trying value " << value << " in range [" << lowerBound << ", "
+             << upperBound << ")";
+      const auto r = callSolver(config);
+      updateResults(configuration, r, results);
+      if (r.sat()) {
+        INFO() << "Found optimum " << value;
+        return;
+      }
+      INFO() << "No solution found. Trying next value.";
+    }
+    INFO() << "No solution found in given interval.";
+  }
+
+  static std::shared_ptr<qc::QuantumComputation>
+  synthesizeSubcircuit(const std::shared_ptr<qc::QuantumComputation>& qc,
+                       std::size_t begin, std::size_t end,
+                       const Configuration& config);
   static void updateResults(const Configuration& config,
                             const Results& newResults, Results& currentResults);
 };

include/cliffordsynthesis/Configuration.hpp

@@ -9,23 +9,31 @@
 #include "nlohmann/json.hpp"
 
 #include <plog/Log.h>
+#include <thread>
+#include <unordered_map>
+#include <variant>
 
 namespace cs {
+
+using SolverParameter = std::variant<bool, std::uint32_t, double, std::string>;
+using SolverParameterMap = std::unordered_map<std::string, SolverParameter>;
+
 struct Configuration {
   Configuration() = default;
 
   /// General configuration for the synthesis algorithm
   std::size_t    initialTimestepLimit    = 0U;
   std::size_t    minimalTimesteps        = 0U;
   bool           useMaxSAT               = false;
+  bool           linearSearch            = false;
   TargetMetric   target                  = TargetMetric::Gates;
   bool           useSymmetryBreaking     = true;
   bool           dumpIntermediateResults = false;
   std::string    intermediateResultsPath = "./";
   plog::Severity verbosity               = plog::Severity::warning;
 
   /// Settings for the SAT solver
-  std::size_t nThreads = 1U;
+  SolverParameterMap solverParameters = {};
 
   /// Settings for depth-optimal synthesis
   bool minimizeGatesAfterDepthOptimization = false;
@@ -35,22 +43,40 @@ struct Configuration {
   double gateLimitFactor                               = 1.1;
   bool   minimizeGatesAfterTwoQubitGateOptimization    = false;
 
+  // Settings for the heuristic solver
+  bool        heuristic         = false;
+  std::size_t splitSize         = 5U;
+  std::size_t nThreadsHeuristic = std::thread::hardware_concurrency();
+
   [[nodiscard]] nlohmann::json json() const {
     nlohmann::json j;
     j["initial_timestep_limit"] = initialTimestepLimit;
     j["minimal_timesteps"]      = minimalTimesteps;
     j["use_max_sat"]            = useMaxSAT;
+    j["linear_search"]          = linearSearch;
     j["target_metric"]          = toString(target);
     j["use_symmetry_breaking"]  = useSymmetryBreaking;
-    j["n_threads"]              = nThreads;
     j["minimize_gates_after_depth_optimization"] =
         minimizeGatesAfterDepthOptimization;
     j["try_higher_gate_limit_for_two_qubit_gate_optimization"] =
         tryHigherGateLimitForTwoQubitGateOptimization;
     j["gate_limit_factor"] = gateLimitFactor;
     j["minimize_gates_after_two_qubit_gate_optimization"] =
         minimizeGatesAfterTwoQubitGateOptimization;
-
+    j["heuristic"]           = heuristic;
+    j["split_size"]          = splitSize;
+    j["n_threads_heuristic"] = nThreadsHeuristic;
+    if (!solverParameters.empty()) {
+      nlohmann::json solverParametersJson;
+      for (const auto& entry : solverParameters) {
+        std::visit(
+            [&solverParametersJson, &entry](const auto& v) {
+              solverParametersJson[entry.first] = v;
+            },
+            entry.second);
+      }
+      j["solver_parameters"] = solverParametersJson;
+    }
     return j;
   }
 

include/cliffordsynthesis/Tableau.hpp

@@ -121,6 +121,8 @@ class Tableau {
     return !(lhs == rhs);
   }
 
+  [[nodiscard]] bool isIdentityTableau() const;
+
   void createDiagonalTableau(std::size_t nQ, bool includeDestabilizers = false);
 
   friend std::ostream& operator<<(std::ostream& os, const Tableau& dt) {

include/cliffordsynthesis/encoding/GateEncoder.hpp

@@ -59,6 +59,7 @@ class GateEncoder {
   static constexpr std::array<qc::OpType, 7> SINGLE_QUBIT_GATES = {
       qc::OpType::None, qc::OpType::X, qc::OpType::Y,   qc::OpType::Z,
       qc::OpType::H,    qc::OpType::S, qc::OpType::Sdag};
+
   [[nodiscard]] static constexpr std::size_t
   gateToIndex(const qc::OpType type) {
     for (std::size_t i = 0; i < SINGLE_QUBIT_GATES.size(); ++i) {
@@ -69,6 +70,36 @@ class GateEncoder {
     return 0;
   }
 
+  template <qc::OpType Gate>
+  [[nodiscard]] static constexpr bool containsGate() {
+    for (const auto& g : // NOLINT(readability-use-anyofallof)
+         SINGLE_QUBIT_GATES) {
+      if (g == Gate) {
+        return true;
+      }
+    }
+    return false;
+  }
+
+  [[nodiscard]] static constexpr bool containsX() {
+    return containsGate<qc::OpType::X>();
+  }
+  [[nodiscard]] static constexpr bool containsY() {
+    return containsGate<qc::OpType::Y>();
+  }
+  [[nodiscard]] static constexpr bool containsZ() {
+    return containsGate<qc::OpType::Z>();
+  }
+  [[nodiscard]] static constexpr bool containsH() {
+    return containsGate<qc::OpType::H>();
+  }
+  [[nodiscard]] static constexpr bool containsS() {
+    return containsGate<qc::OpType::S>();
+  }
+  [[nodiscard]] static constexpr bool containsSdag() {
+    return containsGate<qc::OpType::Sdag>();
+  }
+
 protected:
   // number of qubits N
   std::size_t N{}; // NOLINT (readability-identifier-naming)

include/cliffordsynthesis/encoding/MultiGateEncoder.hpp

@@ -17,7 +17,8 @@ class MultiGateEncoder : public GateEncoder {
   using GateEncoder::GateEncoder;
 
 protected:
-  logicbase::LogicTerm rChanges{};
+  logicbase::LogicTerm   rChanges{};
+  logicbase::LogicMatrix xorHelpers{};
 
   void assertConsistency() const override;
   void assertGateConstraints() override;
@@ -32,6 +33,7 @@ class MultiGateEncoder : public GateEncoder {
                                              std::size_t qubit) override;
   void assertTwoQubitGateOrderConstraints(std::size_t pos, std::size_t ctrl,
                                           std::size_t trgt) override;
+  void splitXorR(const logicbase::LogicTerm& changes, std::size_t pos);
 };
 
 } // namespace cs::encoding

include/cliffordsynthesis/encoding/ObjectiveEncoder.hpp

@@ -37,7 +37,7 @@ class ObjectiveEncoder {
 
   void optimizeGateCount(bool includeSingleQubitGates = true) const;
 
-  void optimizeDepth(bool includeSingleQubitGates = true) const;
+  void optimizeDepth() const;
 
 protected:
   // number of qubits N