✨ Stripping DDs of their identity (#361)

## Description

This PR pulls in the changes from cda-tum/mqt-core#358, which
considerably changes the way matrix decision diagrams are represented.
Particularly, any node resembling the identity is now eliminated and
only implicitly represented.
This further compacts the representation of quantum gates and makes the
identity the most compact it can be---a single terminal node.
The overall performance improvements are still to be evaluated. 

## 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.

---------

Signed-off-by: burgholzer <burgholzer@me.com>

include/DeterministicNoiseSimulator.hpp

@@ -62,7 +62,7 @@ class DeterministicNoiseSimulator: public CircuitSimulator<dd::DensityMatrixSimu
     }
 
     std::map<std::string, std::size_t> measureAllNonCollapsing(std::size_t shots) override {
-        return sampleFromProbabilityMap(rootEdge.getSparseProbabilityVectorStrKeys(measurementThreshold), shots);
+        return sampleFromProbabilityMap(rootEdge.getSparseProbabilityVectorStrKeys(getNumberOfQubits(), measurementThreshold), shots);
     }
 
     void initializeSimulation(std::size_t nQubits) override;

src/HybridSchrodingerFeynmanSimulator.cpp

@@ -50,7 +50,7 @@ qc::VectorDD HybridSchrodingerFeynmanSimulator<Config>::simulateSlicing(std::uni
         sliceDD->garbageCollect();
     }
 
-    auto result = sliceDD->kronecker(upper.edge, lower.edge, false);
+    auto result = sliceDD->kronecker(upper.edge, lower.edge, lower.nqubits, false);
     sliceDD->incRef(result);
 
     return result;
@@ -104,16 +104,19 @@ bool HybridSchrodingerFeynmanSimulator<Config>::Slice::apply(std::unique_ptr<dd:
             isSplitOp          = true;
             const bool control = getNextControl();
             for (const auto& c: opControls) {
-                sliceDD->decRef(edge); // TODO incref and decref could be integrated in delete edge
-                edge = sliceDD->deleteEdge(edge, static_cast<dd::Qubit>(c.qubit), control ? (c.type == qc::Control::Type::Pos ? 0 : 1) : (c.type == qc::Control::Type::Pos ? 1 : 0));
+                auto tmp = edge;
+                edge     = sliceDD->deleteEdge(edge, static_cast<dd::Qubit>(c.qubit), control ? (c.type == qc::Control::Type::Pos ? 0 : 1) : (c.type == qc::Control::Type::Pos ? 1 : 0));
+                // TODO incref and decref could be integrated in delete edge
                 sliceDD->incRef(edge);
+                sliceDD->decRef(tmp);
             }
         } else if (targetInSplit) { // target slice for split or operation in split
             const auto&           param = op->getParameter();
             qc::StandardOperation newOp(nqubits, opControls, opTargets, op->getType(), param, start);
-            sliceDD->decRef(edge);
-            edge = sliceDD->multiply(dd::getDD(&newOp, *sliceDD), edge, static_cast<dd::Qubit>(start));
+            auto                  tmp = edge;
+            edge                      = sliceDD->multiply(dd::getDD(&newOp, *sliceDD), edge, static_cast<dd::Qubit>(start));
             sliceDD->incRef(edge);
+            sliceDD->decRef(tmp);
         }
     } else {
         throw std::invalid_argument("Only StandardOperations are supported for now.");
@@ -130,7 +133,7 @@ std::map<std::string, std::size_t> HybridSchrodingerFeynmanSimulator<Config>::si
     auto splitQubit = static_cast<qc::Qubit>(nqubits / 2);
     if (mode == Mode::DD) {
         simulateHybridTaskflow(splitQubit);
-        return Simulator<Config>::measureAllNonCollapsing(shots);
+        return CircuitSimulator<Config>::measureAllNonCollapsing(shots);
     }
     simulateHybridAmplitudes(splitQubit);
 

src/ShorFastSimulator.cpp

@@ -223,7 +223,7 @@ dd::mEdge ShorFastSimulator<Config>::limitTo(std::uint64_t a) {
 
     for (std::uint32_t p = 1; p < requiredBits + 1; p++) {
         if (((a >> p) & 1U) > 0) {
-            edges[0] = Simulator<Config>::dd->makeIdent(0, static_cast<dd::Qubit>(p - 1));
+            edges[0] = Simulator<Config>::dd->makeIdent();
             edges[3] = f;
         } else {
             edges[0] = f;

src/ShorSimulator.cpp

@@ -255,7 +255,7 @@ dd::mEdge ShorSimulator<Config>::limitTo(std::uint64_t a) {
 
     for (std::uint32_t p = 1; p < requiredBits + 1; p++) {
         if (((a >> p) & 1U) > 0) {
-            edges[0] = Simulator<Config>::dd->makeIdent(0, static_cast<dd::Qubit>(p - 1));
+            edges[0] = Simulator<Config>::dd->makeIdent();
             edges[3] = f;
         } else {
             edges[0] = f;
@@ -349,7 +349,7 @@ dd::mEdge ShorSimulator<Config>::addConstMod(std::uint64_t a) {
 
 template<class Config>
 void ShorSimulator<Config>::uAEmulate(std::uint64_t a, std::int32_t q) {
-    const dd::mEdge limit = Simulator<Config>::dd->makeIdent(0, static_cast<dd::Qubit>(requiredBits - 1));
+    const dd::mEdge limit = Simulator<Config>::dd->makeIdent();
 
     dd::mEdge                f = dd::mEdge::one();
     std::array<dd::mEdge, 4> edges{
@@ -411,7 +411,7 @@ void ShorSimulator<Config>::uAEmulate(std::uint64_t a, std::int32_t q) {
     for (auto i = static_cast<std::int32_t>(2 * requiredBits - 1); i >= 0; --i) {
         if (i == q) {
             edges[1] = edges[2] = dd::mEdge::zero();
-            edges[0]            = Simulator<Config>::dd->makeIdent(0, static_cast<dd::Qubit>(nQubits - static_cast<std::size_t>(i) - 2));
+            edges[0]            = Simulator<Config>::dd->makeIdent();
             edges[3]            = e;
             e                   = Simulator<Config>::dd->makeDDNode(static_cast<dd::Qubit>(nQubits - 1 - static_cast<std::size_t>(i)), edges, false);
         } else {

src/python/bindings.cpp

@@ -69,32 +69,31 @@ std::unique_ptr<Simulator> constructSimulatorWithoutSeed(const py::object& circ,
     return constructSimulator<Simulator>(circ, 1., 1, "fidelity", -1, std::forward<Args>(args)...);
 }
 
-void getNumPyMatrixRec(const qc::MatrixDD& e, const std::complex<dd::fp>& amp, std::size_t i, std::size_t j, std::size_t dim, std::complex<dd::fp>* mat) {
+void getNumPyMatrixRec(const qc::MatrixDD& e, const std::complex<dd::fp>& amp, std::size_t i, std::size_t j, std::size_t dim, std::complex<dd::fp>* mat, std::size_t level) {
     // calculate new accumulated amplitude
-    auto w = std::complex<dd::fp>{dd::RealNumber::val(e.w.r), dd::RealNumber::val(e.w.i)};
-    auto c = amp * w;
+    const auto c = amp * static_cast<std::complex<dd::fp>>(e.w);
 
     // base case
-    if (e.isTerminal()) {
+    if (level == 0U) {
         mat[i * dim + j] = c; // NOLINT(cppcoreguidelines-pro-bounds-pointer-arithmetic)
         return;
     }
 
-    const std::size_t x = i | (1 << e.p->v);
-    const std::size_t y = j | (1 << e.p->v);
-
-    // recursive case
-    if (!e.p->e[0].w.approximatelyZero()) {
-        getNumPyMatrixRec(e.p->e[0], c, i, j, dim, mat);
-    }
-    if (!e.p->e[1].w.approximatelyZero()) {
-        getNumPyMatrixRec(e.p->e[1], c, i, y, dim, mat);
-    }
-    if (!e.p->e[2].w.approximatelyZero()) {
-        getNumPyMatrixRec(e.p->e[2], c, x, j, dim, mat);
+    const auto        nextLevel = static_cast<dd::Qubit>(level - 1U);
+    const std::size_t x         = i | (1 << nextLevel);
+    const std::size_t y         = j | (1 << nextLevel);
+    if (e.isTerminal() || e.p->v < nextLevel) {
+        getNumPyMatrixRec(e, c, i, j, dim, mat, nextLevel);
+        getNumPyMatrixRec(e, c, x, y, dim, mat, nextLevel);
+        return;
     }
-    if (!e.p->e[3].w.approximatelyZero()) {
-        getNumPyMatrixRec(e.p->e[3], c, x, y, dim, mat);
+
+    const auto  coords = {std::pair{i, j}, {i, y}, {x, j}, {x, y}};
+    std::size_t k      = 0U;
+    for (const auto& [a, b]: coords) {
+        if (auto& f = e.p->e[k++]; !f.w.exactlyZero()) {
+            getNumPyMatrixRec(f, c, a, b, dim, mat, nextLevel);
+        }
     }
 }
 
@@ -109,12 +108,12 @@ void getNumPyMatrix(UnitarySimulator<Config>& sim, py::array_t<std::complex<dd::
         throw std::runtime_error("Provided matrix is not a square matrix.");
     }
 
-    const std::size_t dim = 1 << (e.p->v + 1);
+    const std::size_t dim = 1ULL << sim.getNumberOfQubits();
     if (static_cast<std::size_t>(rows) != dim) {
         throw std::runtime_error("Provided matrix does not have the right size.");
     }
 
-    getNumPyMatrixRec(e, std::complex<dd::fp>{1.0, 0.0}, 0, 0, dim, dataPtr);
+    getNumPyMatrixRec(e, std::complex<dd::fp>{1.0, 0.0}, 0, 0, dim, dataPtr, sim.getNumberOfQubits());
 }
 
 void dumpTensorNetwork(const py::object& circ, const std::string& filename) {