Repository: cda-tum/mqt-qmap (language: python)

Issue:
🐛 Missmatches in current implmentation to publication
### mqt.qmap version

mqt-qmap 2.1.2 via pip

### OS

Ubuntu

### Python version

_No response_

### C++ compiler

_No response_

### Additional environment information

_No response_

### Description

Please find below an (incomplete) list of benchmarks where the current implementation of the **exact** mapping algorithm shows differences when reproducing the results of [Wille et al.](https://doi.org/10.1145/3316781.3317859) (IBM QX4 arch):

- 3_17_13.qasm
- 4gt11_83.qasm
- 4mod5-v0_19.qasm
- alu-v3_34.qasm



### Expected behavior

I would expect the same results as mentioned in the paper.

### How to Reproduce

Run listed benchmarks above with exact mapper for arch graph resembling IBM QX4 (asymmetric) connectivity and compare with results of publication.
🐛 Missmatches in current implmentation to publication
### mqt.qmap version

mqt-qmap 2.1.2 via pip

### OS

Ubuntu

### Python version

_No response_

### C++ compiler

_No response_

### Additional environment information

_No response_

### Description

Please find below an (incomplete) list of benchmarks where the current implementation of the **exact** mapping algorithm shows differences when reproducing the results of [Wille et al.](https://doi.org/10.1145/3316781.3317859) (IBM QX4 arch):

- 3_17_13.qasm
- 4gt11_83.qasm
- 4mod5-v0_19.qasm
- alu-v3_34.qasm



### Expected behavior

I would expect the same results as mentioned in the paper.

### How to Reproduce

Run listed benchmarks above with exact mapper for arch graph resembling IBM QX4 (asymmetric) connectivity and compare with results of publication.


Discussion / hints:
Thanks for opening this issue. I'll try to look into it in more detail over the next couple of days.
A note on that right away: Since the original publication, QMAP as a tool has evolved quite a bit. So I wouldn't be too surprised, if the results it produces changed over time.
If the results got better, I'd even go as far and say that's desirable.

Could you provide me with a little more context on the difference in results you are seeing?
Is it final gate count, number of SWAPs, number of direction reverses? Is it better or worse?
I checked against the original publication, and I find the following differences based on the benchmarks shown on the publication.
I count the total number of gates after mapping (as in the publication).
The format below is: 
{benchmark name}: {optimal results from publication} vs. {qmap_exact}

3_17_13.qasm: 59 vs. 60
4gt11_83.qasm: 49 vs. 50
4mod5-v0_19.qasm: 64 vs. 69
alu-v3_34.qasm: 90 vs. 91
mod5mils_65.qasm: 64 vs. 65
rd32-v0_66.qasm: 63 vs. 73
rd32-v1_68.qasm: 65 vs. 75

I currently see no reason to argue that the results in the publication are not correct. Maybe there is a bug in the permutation function that SWAPs some qubits even though that is not necessary for physical connectivity reasons?

Thanks for the details. Now that I see the numbers, I think I remember the cause for the difference.
IIRC, it comes down to a very subtle difference in the coupling map shown in the paper and the _real_ coupling map of the IBM QX4 device (as it can be seen in https://arxiv.org/abs/1808.05661 and exposed through QMAP).

In the paper, the coupling map for the IBM QX4 device shown in Figure 2 is:
```
[[1, 0], [2, 0], [2, 1], [3, 2], [3, 4], [4, 2]]
```
whereas `AvailableArchitecture::IbmQx4` in QMAP is defined as:
```
[[1, 0], [2, 0], [2, 1], [3, 2], [3, 4], [2, 4]]
```
Note the reversed direction of the last edge of the coupling map.

The results in the paper were generated based on the coupling map shown in the paper. This can also be seen in the original code repository at https://github.com/iic-jku/minimal_ibm_qx_mapping/blob/e41de077a1c1096be74a9bcdd9ced35228bd69d5/minimal_ibm_qx_mapping/main.cpp#L79-L86.
The coupling map error has since been corrected in QMAP (see above), but there is no way we could update the publication.

If you are using the Python interface of QMAP, using the old Architecture is as easy as
```python
from mqt import qmap

circ = ...
arch = qmap.Architecture(
    5,
    {
        (1, 0),
        (2, 0),
        (2, 1),
        (3, 2),
        (3, 4),
        (4, 2)
    },
)
qc_mapped, res = qmap.compile(qc, arch, method="exact")
```

Please let us know if you can replicate the results based on the above suggestions.
I already included that "mismatch", however, results are still not the optimal ones as in the publication.
Note that is actually used 
`qmap.compile(circ, arch, method="exact", post_mapping_optimizations=False)`
I will have a look at the _rd32-*_ circuits, as here the difference seems to be the largest. Frome these circuits, I will try to derive a smaller circuit as a MVE and potentially create an example where qmap (in _excact_ mode) does not find the minimal solution (and provide this minimal solution independently). Hopefully this will help to narrow down if there really is an issue or not...
Thanks for opening this issue. I'll try to look into it in more detail over the next couple of days.
A note on that right away: Since the original publication, QMAP as a tool has evolved quite a bit. So I wouldn't be too surprised, if the results it produces changed over time.
If the results got better, I'd even go as far and say that's desirable.

Could you provide me with a little more context on the difference in results you are seeing?
Is it final gate count, number of SWAPs, number of direction reverses? Is it better or worse?
I checked against the original publication, and I find the following differences based on the benchmarks shown on the publication.
I count the total number of gates after mapping (as in the publication).
The format below is: 
{benchmark name}: {optimal results from publication} vs. {qmap_exact}

3_17_13.qasm: 59 vs. 60
4gt11_83.qasm: 49 vs. 50
4mod5-v0_19.qasm: 64 vs. 69
alu-v3_34.qasm: 90 vs. 91
mod5mils_65.qasm: 64 vs. 65
rd32-v0_66.qasm: 63 vs. 73
rd32-v1_68.qasm: 65 vs. 75

I currently see no reason to argue that the results in the publication are not correct. Maybe there is a bug in the permutation function that SWAPs some qubits even though that is not necessary for physical connectivity reasons?

Thanks for the details. Now that I see the numbers, I think I remember the cause for the difference.
IIRC, it comes down to a very subtle difference in the coupling map shown in the paper and the _real_ coupling map of the IBM QX4 device (as it can be seen in https://arxiv.org/abs/1808.05661 and exposed through QMAP).

In the paper, the coupling map for the IBM QX4 device shown in Figure 2 is:
```
[[1, 0], [2, 0], [2, 1], [3, 2], [3, 4], [4, 2]]
```
whereas `AvailableArchitecture::IbmQx4` in QMAP is defined as:
```
[[1, 0], [2, 0], [2, 1], [3, 2], [3, 4], [2, 4]]
```
Note the reversed direction of the last edge of the coupling map.

The results in the paper were generated based on the coupling map shown in the paper. This can also be seen in the original code repository at https://github.com/iic-jku/minimal_ibm_qx_mapping/blob/e41de077a1c1096be74a9bcdd9ced35228bd69d5/minimal_ibm_qx_mapping/main.cpp#L79-L86.
The coupling map error has since been corrected in QMAP (see above), but there is no way we could update the publication.

If you are using the Python interface of QMAP, using the old Architecture is as easy as
```python
from mqt import qmap

circ = ...
arch = qmap.Architecture(
    5,
    {
        (1, 0),
        (2, 0),
        (2, 1),
        (3, 2),
        (3, 4),
        (4, 2)
    },
)
qc_mapped, res = qmap.compile(qc, arch, method="exact")
```

Please let us know if you can replicate the results based on the above suggestions.
I already included that "mismatch", however, results are still not the optimal ones as in the publication.
Note that is actually used 
`qmap.compile(circ, arch, method="exact", post_mapping_optimizations=False)`
I will have a look at the _rd32-*_ circuits, as here the difference seems to be the largest. Frome these circuits, I will try to derive a smaller circuit as a MVE and potentially create an example where qmap (in _excact_ mode) does not find the minimal solution (and provide this minimal solution independently). Hopefully this will help to narrow down if there really is an issue or not...

Fix patch:
diff --git a/include/Architecture.hpp b/include/Architecture.hpp
index 2487063fd..d00c4d969 100644
--- a/include/Architecture.hpp
+++ b/include/Architecture.hpp
@@ -15,6 +15,7 @@
 #include <map>
 #include <regex>
 #include <unordered_map>
+#include <unordered_set>
 #include <vector>
 
 constexpr std::uint8_t GATES_OF_BIDIRECTIONAL_SWAP  = 3U;
@@ -390,10 +391,10 @@ class Architecture {
   static std::size_t findCouplingLimit(const CouplingMap& cm,
                                        std::uint16_t      nQubits);
   static std::size_t
-  findCouplingLimit(const CouplingMap& cm, std::uint16_t nQubits,
-                    const std::set<std::uint16_t>& qubitChoice);
-  static void
-  findCouplingLimit(std::uint16_t node, std::uint16_t curSum,
-                    const std::vector<std::vector<std::uint16_t>>& connections,
-                    std::vector<std::uint16_t>& d, std::vector<bool>& visited);
+              findCouplingLimit(const CouplingMap& cm, std::uint16_t nQubits,
+                                const std::set<std::uint16_t>& qubitChoice);
+  static void findCouplingLimit(
+      std::uint16_t node, std::uint16_t curSum,
+      const std::vector<std::unordered_set<std::uint16_t>>& connections,
+      std::vector<std::uint16_t>& d, std::vector<bool>& visited);
 };
diff --git a/include/configuration/Configuration.hpp b/include/configuration/Configuration.hpp
index bfd7c1fa0..ada24f871 100644
--- a/include/configuration/Configuration.hpp
+++ b/include/configuration/Configuration.hpp
@@ -73,5 +73,8 @@ struct Configuration {
   [[nodiscard]] nlohmann::json json() const;
   [[nodiscard]] std::string    toString() const { return json().dump(2); }
 
-  void setTimeout(const std::size_t sec) { timeout = sec; }
+  void               setTimeout(const std::size_t sec) { timeout = sec; }
+  [[nodiscard]] bool swapLimitsEnabled() const {
+    return (swapReduction != SwapReduction::None) && enableSwapLimits;
+  }
 };
diff --git a/src/Architecture.cpp b/src/Architecture.cpp
index 4b352ba6d..89c065d61 100644
--- a/src/Architecture.cpp
+++ b/src/Architecture.cpp
@@ -486,13 +486,15 @@ std::uint64_t Architecture::bfs(const std::uint16_t   start,
 
 std::size_t Architecture::findCouplingLimit(const CouplingMap&  cm,
                                             const std::uint16_t nQubits) {
-  std::vector<std::vector<std::uint16_t>> connections;
-  std::vector<std::uint16_t>              d;
-  std::vector<bool>                       visited;
+  std::vector<std::unordered_set<std::uint16_t>> connections;
+  std::vector<std::uint16_t>                     d;
+  std::vector<bool>                              visited;
   connections.resize(nQubits);
   std::uint16_t maxSum = 0;
   for (const auto& edge : cm) {
-    connections.at(edge.first).emplace_back(edge.second);
+    connections.at(edge.first).emplace(edge.second);
+    // make sure that the connections are bidirectional
+    connections.at(edge.second).emplace(edge.first);
   }
   for (std::uint16_t q = 0; q < nQubits; ++q) {
     d.clear();
@@ -513,15 +515,17 @@ std::size_t Architecture::findCouplingLimit(const CouplingMap&  cm,
 std::size_t Architecture::findCouplingLimit(const CouplingMap&  cm,
                                             const std::uint16_t nQubits,
                                             const QubitSubset&  qubitChoice) {
-  std::vector<std::vector<std::uint16_t>> connections;
-  std::vector<std::uint16_t>              d;
-  std::vector<bool>                       visited;
+  std::vector<std::unordered_set<std::uint16_t>> connections;
+  std::vector<std::uint16_t>                     d;
+  std::vector<bool>                              visited;
   connections.resize(nQubits);
   std::uint16_t maxSum = 0;
   for (const auto& edge : cm) {
     if ((qubitChoice.count(edge.first) != 0U) &&
         (qubitChoice.count(edge.second) != 0U)) {
-      connections.at(edge.first).emplace_back(edge.second);
+      connections.at(edge.first).emplace(edge.second);
+      // make sure that the connections are bidirectional
+      connections.at(edge.second).emplace(edge.first);
     }
   }
   for (std::uint16_t q = 0; q < nQubits; ++q) {
@@ -545,7 +549,7 @@ std::size_t Architecture::findCouplingLimit(const CouplingMap&  cm,
 
 void Architecture::findCouplingLimit(
     const std::uint16_t node, const std::uint16_t curSum,
-    const std::vector<std::vector<std::uint16_t>>& connections,
+    const std::vector<std::unordered_set<std::uint16_t>>& connections,
     std::vector<std::uint16_t>& d, std::vector<bool>& visited) {
   if (visited.at(node)) {
     return;
diff --git a/src/exact/ExactMapper.cpp b/src/exact/ExactMapper.cpp
index 3896e6705..28fdf9e5d 100644
--- a/src/exact/ExactMapper.cpp
+++ b/src/exact/ExactMapper.cpp
@@ -104,6 +104,12 @@ void ExactMapper::map(const Configuration& settings) {
       maxLimit = architecture.getCouplingLimit() - 1U;
     }
     if (config.swapReduction == SwapReduction::CouplingLimit) {
+      if (!architecture.bidirectional()) {
+        // on a directed architecture, one more SWAP might be needed overall
+        // due to the directionality of the edges and direction reversal not
+        // being possible for every gate.
+        maxLimit += 1U;
+      }
       limit = maxLimit;
     } else if (config.swapReduction == SwapReduction::Increasing) {
       limit = 0U;
@@ -401,7 +407,7 @@ void ExactMapper::coreMappingRoutine(
   //////////////////////////////////////////
   /// 	Check necessary permutations	//
   //////////////////////////////////////////
-  if (config.enableSwapLimits && !config.useBDD) {
+  if (config.swapLimitsEnabled() && !config.useBDD) {
     do {
       auto picost = architecture.minimumNumberOfSwaps(
           pi, static_cast<std::int64_t>(limit));
@@ -449,7 +455,7 @@ number of variables: (|L|-1) * m!
     y.emplace_back();
     piCount = 0;
     do {
-      if (skippedPi.count(piCount) == 0 || !config.enableSwapLimits) {
+      if (skippedPi.count(piCount) == 0 || !config.swapLimitsEnabled()) {
         yName.str("");
         yName << "y_" << k << '_' << piCount;
         y.back().emplace_back(lb->makeVariable(yName.str(), CType::BOOL));
@@ -621,7 +627,7 @@ number of variables: (|L|-1) * m!
     auto& i         = x[k - 1];
     auto& j         = x[k];
     do {
-      if (skippedPi.count(piCount) == 0 || !config.enableSwapLimits) {
+      if (skippedPi.count(piCount) == 0 || !config.swapLimitsEnabled()) {
         auto equal = LogicTerm(true);
         for (const auto qubit : qubitChoice) {
           for (std::size_t q = 0; q < qc.getNqubits(); ++q) {
@@ -645,7 +651,7 @@ number of variables: (|L|-1) * m!
       piCount         = 0;
       internalPiCount = 0;
       do {
-        if (skippedPi.count(piCount) == 0 || !config.enableSwapLimits) {
+        if (skippedPi.count(piCount) == 0 || !config.swapLimitsEnabled()) {
           onlyOne = onlyOne + LogicTerm::ite(y[k - 1][internalPiCount],
                                              LogicTerm(1), LogicTerm(0));
           ++internalPiCount;
@@ -660,7 +666,7 @@ number of variables: (|L|-1) * m!
       piCount         = 0;
       internalPiCount = 0;
       do {
-        if (skippedPi.count(piCount) == 0 || !config.enableSwapLimits) {
+        if (skippedPi.count(piCount) == 0 || !config.swapLimitsEnabled()) {
           varIDs.push_back(y[k - 1][internalPiCount]);
           ++internalPiCount;
         }
@@ -694,7 +700,7 @@ number of variables: (|L|-1) * m!
       reducedLayerIndices.size());
   auto cost = LogicTerm(0);
   do {
-    if (skippedPi.count(piCount) == 0 || !config.enableSwapLimits) {
+    if (skippedPi.count(piCount) == 0 || !config.swapLimitsEnabled()) {
       auto picost = architecture.minimumNumberOfSwaps(pi);
       if (architecture.bidirectional()) {
         picost *= GATES_OF_BIDIRECTIONAL_SWAP;
@@ -702,9 +708,8 @@ number of variables: (|L|-1) * m!
         picost *= GATES_OF_UNIDIRECTIONAL_SWAP;
       }
       for (std::size_t k = 1; k < reducedLayerIndices.size(); ++k) {
-        cost = cost + LogicTerm::ite(y[k - 1][internalPiCount],
-                                     LogicTerm(static_cast<int>(picost)),
-                                     LogicTerm(0));
+        lb->weightedTerm(y[k - 1][internalPiCount],
+                         static_cast<double>(picost));
         if (config.useBDD) {
           weightedVars[k].emplace(y[k - 1][internalPiCount],
                                   static_cast<int>(picost));
@@ -714,8 +719,7 @@ number of variables: (|L|-1) * m!
     }
     ++piCount;
   } while (std::next_permutation(pi.begin(), pi.end()));
-  lb->minimize(cost);
-  if (config.enableSwapLimits && config.useBDD) {
+  if (config.swapLimitsEnabled() && config.useBDD) {
     for (std::size_t k = 1; k < reducedLayerIndices.size(); ++k) {
       lb->assertFormula(BuildBDD(weightedVars[k], y[k - 1],
                                  static_cast<int>(limit), lb.get()));
@@ -724,7 +728,6 @@ number of variables: (|L|-1) * m!
 
   // cost for reversed directions
   if (!architecture.bidirectional()) {
-    cost                 = LogicTerm(0);
     const auto numLayers = reducedLayerIndices.size();
     for (std::size_t k = 0; k < numLayers; ++k) {
       for (const auto& gate : layers.at(reducedLayerIndices.at(k))) {
@@ -739,13 +742,11 @@ number of variables: (|L|-1) * m!
                                 [static_cast<std::size_t>(gate.control)];
           reverse = reverse || (indexFT && indexSC);
         }
-        cost = cost + LogicTerm::ite(reverse,
-                                     LogicTerm(::GATES_OF_DIRECTION_REVERSE),
-                                     LogicTerm(0));
+        lb->weightedTerm(reverse, ::GATES_OF_DIRECTION_REVERSE);
       }
     }
-    lb->minimize(cost);
   }
+  lb->makeMinimize();
 
   if (config.includeWCNF) {
     choiceResults.wcnf = lb->dumpInternalSolver();
@@ -807,7 +808,7 @@ number of variables: (|L|-1) * m!
         // sort the permutation of the qubits to start fresh
         std::sort(pi.begin(), pi.end());
         do {
-          if (skippedPi.count(piCount) == 0 || !config.enableSwapLimits) {
+          if (skippedPi.count(piCount) == 0 || !config.swapLimitsEnabled()) {
             if (m->getBoolValue(y[k - 1][internalPiCount], lb.get())) {
               break;
             }


Test patch:
diff --git a/test/test_exact.cpp b/test/test_exact.cpp
index 45133f408..d3ae23dc1 100644
--- a/test/test_exact.cpp
+++ b/test/test_exact.cpp
@@ -534,7 +534,7 @@ TEST_F(ExactTest, Test4QCircuitThatUsesAll5Q) {
   EXPECT_EQ(results.output.swaps, 1);
 }
 
-TEST_F(ExactTest, Test) {
+TEST_F(ExactTest, RegressionTestDirectionReverseCost) {
   // Regression test for https://github.com/cda-tum/qmap/issues/251
   using namespace qc::literals;
 
@@ -556,3 +556,76 @@ TEST_F(ExactTest, Test) {
   EXPECT_EQ(mapper.getResults().output.swaps, 0);
   EXPECT_EQ(mapper.getResults().output.directionReverse, 2);
 }
+
+TEST_F(ExactTest, RegressionTestExactMapperPerformance) {
+  // Regression test for https://github.com/cda-tum/qmap/issues/256
+  std::stringstream ss{"OPENQASM 2.0;\n"
+                       "include \"qelib1.inc\";\n"
+                       "qreg q[3];\n"
+                       "cx q[0],q[2];\n"
+                       "cx q[2],q[1];\n"
+                       "cx q[2],q[1];\n"
+                       "cx q[0],q[2];\n"
+                       "cx q[1],q[0];\n"
+                       "cx q[1],q[2];\n"
+                       "cx q[0],q[2];\n"
+                       "cx q[1],q[0];\n"
+                       "cx q[2],q[1];\n"
+                       "cx q[1],q[0];\n"
+                       "cx q[2],q[1];\n"
+                       "cx q[0],q[2];\n"
+                       "cx q[0],q[1];\n"
+                       "cx q[2],q[1];\n"
+                       "cx q[0],q[2];\n"
+                       "cx q[1],q[0];\n"
+                       "cx q[1],q[2];\n"};
+
+  Architecture      arch;
+  const CouplingMap cm = {{1, 0}, {2, 0}, {2, 1}, {3, 2}, {3, 4}, {4, 2}};
+  arch.loadCouplingMap(5, cm);
+  qc.import(ss, qc::Format::OpenQASM);
+
+  auto mapper            = ExactMapper(qc, arch);
+  settings.swapReduction = SwapReduction::CouplingLimit;
+  mapper.map(settings);
+  EXPECT_EQ(mapper.getResults().output.swaps, 1);
+  EXPECT_EQ(mapper.getResults().output.directionReverse, 4);
+
+  auto mapper2           = ExactMapper(qc, arch);
+  settings.swapReduction = SwapReduction::None;
+  mapper2.map(settings);
+  EXPECT_EQ(mapper2.getResults().output.swaps, 1);
+  EXPECT_EQ(mapper2.getResults().output.directionReverse, 4);
+}
+
+TEST_F(ExactTest, RegressionTestExactMapperPerformance2) {
+  // Regression test for https://github.com/cda-tum/qmap/issues/256
+  std::stringstream ss{"OPENQASM 2.0;\n"
+                       "include \"qelib1.inc\";\n"
+                       "qreg q[4];\n"
+                       "cx q[0],q[1];\n"
+                       "cx q[3],q[0];\n"
+                       "cx q[1],q[3];\n"
+                       "cx q[1],q[0];\n"
+                       "cx q[3],q[0];\n"
+                       "cx q[1],q[3];\n"
+                       "cx q[0],q[1];\n"
+                       "cx q[1],q[2];\n"};
+
+  Architecture      arch;
+  const CouplingMap cm = {{1, 0}, {2, 0}, {2, 1}, {3, 2}, {3, 4}, {4, 2}};
+  arch.loadCouplingMap(5, cm);
+  qc.import(ss, qc::Format::OpenQASM);
+
+  auto mapper            = ExactMapper(qc, arch);
+  settings.swapReduction = SwapReduction::CouplingLimit;
+  mapper.map(settings);
+  EXPECT_EQ(mapper.getResults().output.swaps, 1);
+  EXPECT_EQ(mapper.getResults().output.directionReverse, 1);
+
+  auto mapper2           = ExactMapper(qc, arch);
+  settings.swapReduction = SwapReduction::None;
+  mapper2.map(settings);
+  EXPECT_EQ(mapper2.getResults().output.swaps, 1);
+  EXPECT_EQ(mapper2.getResults().output.directionReverse, 1);
+}