⬆️ Update mqt-core submodule

apps/noise_aware.cpp

@@ -118,7 +118,7 @@ int main(int argc, char** argv) { // NOLINT(bugprone-exception-escape)
 
         auto t1 = std::chrono::steady_clock::now();
 
-        const std::map<std::string, double> measurementResults = ddsim->deterministicSimulate();
+        const auto measurementResults = ddsim->deterministicSimulate();
 
         auto t2 = std::chrono::steady_clock::now();
 
@@ -146,6 +146,6 @@ int main(int argc, char** argv) { // NOLINT(bugprone-exception-escape)
         if (vm.count("pm") > 0) {
             outputObj["measurement_results"] = measurementResults;
         }
-        std::cout << std::setw(2) << outputObj << std::endl;
+        std::cout << std::setw(2) << outputObj << "\n";
     }
 }

apps/simple.cpp

@@ -19,6 +19,18 @@
 
 namespace nl = nlohmann;
 
+namespace std {
+    template<class T>
+    void to_json(nl::json& j, const std::complex<T>& p) { // NOLINT(readability-identifier-naming)
+        j = nl::json{p.real(), p.imag()};
+    }
+    template<class T>
+    void from_json(const nl::json& j, std::complex<T>& p) { // NOLINT(readability-identifier-naming)
+        p.real(j.at(0));
+        p.imag(j.at(1));
+    }
+} /* namespace std */
+
 int main(int argc, char** argv) { // NOLINT(bugprone-exception-escape)
     cxxopts::Options options("MQT DDSIM", "for more information see https://www.cda.cit.tum.de/");
     // clang-format off
@@ -211,9 +223,9 @@ int main(int argc, char** argv) { // NOLINT(bugprone-exception-escape)
 
     if (vm.count("pv") > 0) {
         if (auto* hsfSim = dynamic_cast<HybridSchrodingerFeynmanSimulator<>*>(ddsim.get())) {
-            outputObj["state_vector"] = hsfSim->getVectorFromHybridSimulation<std::pair<dd::fp, dd::fp>>();
+            outputObj["state_vector"] = hsfSim->getVectorFromHybridSimulation();
         } else {
-            outputObj["state_vector"] = ddsim->getVector<std::pair<dd::fp, dd::fp>>();
+            outputObj["state_vector"] = ddsim->getVector();
         }
     }
 
@@ -245,5 +257,5 @@ int main(int argc, char** argv) { // NOLINT(bugprone-exception-escape)
         dd::exportEdgeWeights(ddsim->rootEdge, ostream);
     }
 
-    std::cout << std::setw(2) << outputObj << std::endl;
+    std::cout << std::setw(2) << outputObj << "\n";
 }

include/DeterministicNoiseSimulator.hpp

@@ -36,27 +36,27 @@ class DeterministicNoiseSimulator: public Simulator<Config> {
         return sampleFromProbabilityMap(deterministicSimulate(), shots);
     };
 
-    std::map<std::string, dd::fp> deterministicSimulate();
+    dd::SparsePVecStrKeys deterministicSimulate();
 
-    std::map<std::string, std::size_t> sampleFromProbabilityMap(const std::map<std::string, dd::fp>& resultProbabilityMap, std::size_t shots);
+    std::map<std::string, std::size_t> sampleFromProbabilityMap(const dd::SparsePVecStrKeys& resultProbabilityMap, std::size_t shots);
 
     [[nodiscard]] std::size_t getNumberOfQubits() const override { return qc->getNqubits(); };
 
     [[nodiscard]] std::size_t getNumberOfOps() const override { return qc->getNops(); };
 
     [[nodiscard]] std::string getName() const override { return qc->getName(); };
 
-    [[nodiscard]] std::size_t getActiveNodeCount() const override { return Simulator<Config>::dd->template getUniqueTable<dd::dNode>().getStats().activeEntryCount; }
-    [[nodiscard]] std::size_t getMaxNodeCount() const override { return Simulator<Config>::dd->template getUniqueTable<dd::dNode>().getStats().peakActiveEntryCount; }
+    [[nodiscard]] std::size_t getActiveNodeCount() const override { return Simulator<Config>::dd->template getUniqueTable<dd::dNode>().getNumActiveEntries(); }
+    [[nodiscard]] std::size_t getMaxNodeCount() const override { return Simulator<Config>::dd->template getUniqueTable<dd::dNode>().getPeakNumActiveEntries(); }
 
     [[nodiscard]] std::size_t countNodesFromRoot() override {
         if (useDensityMatrixType) {
             qc::DensityMatrixDD::alignDensityEdge(rootEdge);
-            const std::size_t tmp = Simulator<Config>::dd->size(rootEdge);
+            const std::size_t tmp = rootEdge.size();
             qc::DensityMatrixDD::setDensityMatrixTrue(rootEdge);
             return tmp;
         }
-        return Simulator<Config>::dd->size(rootEdge);
+        return rootEdge.size();
     }
 
     qc::DensityMatrixDD rootEdge{};

include/HybridSchrodingerFeynmanSimulator.hpp

@@ -49,21 +49,15 @@ class HybridSchrodingerFeynmanSimulator: public CircuitSimulator<Config> {
 
     Mode mode = Mode::Amplitude;
 
-    template<class ReturnType = dd::ComplexValue>
-    [[nodiscard]] std::vector<ReturnType> getVectorFromHybridSimulation() const {
+    [[nodiscard]] dd::CVec getVectorFromHybridSimulation() const {
         if (CircuitSimulator<Config>::getNumberOfQubits() >= 60) {
             // On 64bit system the vector can hold up to (2^60)-1 elements, if memory permits
             throw std::range_error("getVector only supports less than 60 qubits.");
         }
         if (getMode() == Mode::Amplitude) {
-            if constexpr (std::is_same_v<ReturnType, decltype(finalAmplitudes)>) {
-                return finalAmplitudes;
-            }
-            std::vector<ReturnType> amplitudes;
-            std::transform(finalAmplitudes.begin(), finalAmplitudes.end(), std::back_inserter(amplitudes), [](std::complex<dd::fp> x) -> ReturnType { return {x.real(), x.imag()}; });
-            return amplitudes;
+            return finalAmplitudes;
         }
-        return CircuitSimulator<Config>::template getVector<ReturnType>();
+        return CircuitSimulator<Config>::getVector();
     }
 
     //  Get # of decisions for given split_qubit, so that lower slice: q0 < i < qubit; upper slice: qubit <= i < nqubits
@@ -72,8 +66,8 @@ class HybridSchrodingerFeynmanSimulator: public CircuitSimulator<Config> {
     [[nodiscard]] Mode getMode() const { return mode; }
 
 private:
-    std::size_t                       nthreads = 2;
-    std::vector<std::complex<dd::fp>> finalAmplitudes{};
+    std::size_t nthreads = 2;
+    dd::CVec    finalAmplitudes{};
 
     void simulateHybridTaskflow(qc::Qubit splitQubit);
     void simulateHybridAmplitudes(qc::Qubit splitQubit);

include/Simulator.hpp

@@ -66,33 +66,23 @@ class Simulator {
 
     std::map<std::string, std::size_t> sampleFromAmplitudeVectorInPlace(std::vector<std::complex<dd::fp>>& amplitudes, std::size_t shots);
 
-    template<class ReturnType = dd::ComplexValue>
-    [[nodiscard]] std::vector<ReturnType> getVector() const {
+    [[nodiscard]] dd::CVec getVector() const {
         if (getNumberOfQubits() >= 60) {
             // On 64bit system the vector can hold up to (2^60)-1 elements, if memory permits
             throw std::range_error("getVector only supports less than 60 qubits.");
         }
-        std::string             path(getNumberOfQubits(), '0');
-        std::vector<ReturnType> results;
-        results.resize(1ULL << getNumberOfQubits());
-        for (std::size_t i = 0; i < 1ULL << getNumberOfQubits(); ++i) {
-            const std::string      correctedPath{path.rbegin(), path.rend()};
-            const dd::ComplexValue cv = dd->getValueByPath(rootEdge, correctedPath);
-            results[i]                = {cv.r, cv.i};
-            nextPath(path);
-        }
-        return results;
+        return rootEdge.getVector();
     }
 
-    [[nodiscard]] virtual std::size_t getActiveNodeCount() const { return dd->template getUniqueTable<dd::vNode>().getStats().activeEntryCount; }
+    [[nodiscard]] virtual std::size_t getActiveNodeCount() const { return dd->template getUniqueTable<dd::vNode>().getNumActiveEntries(); }
 
-    [[nodiscard]] virtual std::size_t getMaxNodeCount() const { return dd->template getUniqueTable<dd::vNode>().getStats().peakActiveEntryCount; }
+    [[nodiscard]] virtual std::size_t getMaxNodeCount() const { return dd->template getUniqueTable<dd::vNode>().getPeakNumActiveEntries(); }
 
-    [[nodiscard]] virtual std::size_t getMaxMatrixNodeCount() const { return dd->template getUniqueTable<dd::mNode>().getStats().activeEntryCount; }
+    [[nodiscard]] virtual std::size_t getMaxMatrixNodeCount() const { return dd->template getUniqueTable<dd::mNode>().getPeakNumActiveEntries(); }
 
-    [[nodiscard]] virtual std::size_t getMatrixActiveNodeCount() const { return dd->template getUniqueTable<dd::mNode>().getStats().activeEntryCount; }
+    [[nodiscard]] virtual std::size_t getMatrixActiveNodeCount() const { return dd->template getUniqueTable<dd::mNode>().getNumActiveEntries(); }
 
-    [[nodiscard]] virtual std::size_t countNodesFromRoot() { return dd->size(rootEdge); }
+    [[nodiscard]] virtual std::size_t countNodesFromRoot() { return rootEdge.size(); }
 
     [[nodiscard]] std::pair<dd::ComplexValue, std::string> getPathOfLeastResistance() const;
 
@@ -104,16 +94,6 @@ class Simulator {
 
     [[nodiscard]] virtual std::string getName() const = 0;
 
-    [[nodiscard]] static inline std::string toBinaryString(const std::size_t value, const std::size_t numberOfQubits) {
-        std::string binary(numberOfQubits, '0');
-        for (std::size_t j = 0; j < numberOfQubits; ++j) {
-            if ((value & (1U << j)) != 0U) {
-                binary[numberOfQubits - 1 - j] = '1';
-            }
-        }
-        return binary;
-    }
-
     void setTolerance(const dd::fp tolerance) {
         dd->cn.setTolerance(tolerance);
     }
@@ -144,8 +124,6 @@ class Simulator {
     std::uint64_t seed = 0;
     bool          hasFixedSeed;
     dd::fp        epsilon = 0.001;
-
-    static void nextPath(std::string& s);
 };
 
 struct StochasticNoiseSimulatorDDPackageConfig: public dd::DDPackageConfig {

include/UnitarySimulator.hpp

@@ -45,8 +45,8 @@ class UnitarySimulator: public CircuitSimulator<Config> {
     [[nodiscard]] Mode         getMode() const { return mode; }
     [[nodiscard]] qc::MatrixDD getConstructedDD() const { return e; }
     [[nodiscard]] double       getConstructionTime() const { return constructionTime; }
-    [[nodiscard]] std::size_t  getFinalNodeCount() const { return Simulator<Config>::dd->size(e); }
-    [[nodiscard]] std::size_t  getMaxNodeCount() const override { return Simulator<Config>::dd->template getUniqueTable<dd::mNode>().getStats().peakActiveEntryCount; }
+    [[nodiscard]] std::size_t  getFinalNodeCount() const { return e.size(); }
+    [[nodiscard]] std::size_t  getMaxNodeCount() const override { return Simulator<Config>::dd->template getUniqueTable<dd::mNode>().getPeakNumActiveEntries(); }
 
 private:
     qc::MatrixDD e{};

src/CircuitSimulator.cpp

@@ -185,7 +185,7 @@ std::map<std::size_t, bool> CircuitSimulator<Config>::singleShot(const bool igno
             if (approximationInfo.stepNumber > 0 && approximationInfo.stepFidelity < 1.0) {
                 if (approximationInfo.strategy == ApproximationInfo::FidelityDriven && (opNum + 1) % approxMod == 0 &&
                     approximationRuns < approximationInfo.stepNumber) {
-                    [[maybe_unused]] const auto sizeBefore = Simulator<Config>::dd->size(Simulator<Config>::rootEdge);
+                    [[maybe_unused]] const auto sizeBefore = Simulator<Config>::rootEdge.size();
                     const auto                  apFid      = Simulator<Config>::approximateByFidelity(approximationInfo.stepFidelity, false, true);
                     approximationRuns++;
                     finalFidelity *= static_cast<long double>(apFid);
@@ -197,7 +197,7 @@ std::map<std::size_t, bool> CircuitSimulator<Config>::singleShot(const bool igno
                               << "; #runs=" << approximation_runs
                               << "\n";//*/
                 } else if (approximationInfo.strategy == ApproximationInfo::MemoryDriven) {
-                    [[maybe_unused]] const auto sizeBefore = Simulator<Config>::dd->size(Simulator<Config>::rootEdge);
+                    [[maybe_unused]] const auto sizeBefore = Simulator<Config>::rootEdge.size();
                     if (Simulator<Config>::dd->template getUniqueTable<dd::vNode>().possiblyNeedsCollection()) {
                         const auto apFid = Simulator<Config>::approximateByFidelity(approximationInfo.stepFidelity, false, true);
                         approximationRuns++;

src/DeterministicNoiseSimulator.cpp

@@ -5,7 +5,7 @@
 using CN = dd::ComplexNumbers;
 
 template<class Config>
-std::map<std::string, double> DeterministicNoiseSimulator<Config>::deterministicSimulate() {
+dd::SparsePVecStrKeys DeterministicNoiseSimulator<Config>::deterministicSimulate() {
     rootEdge = Simulator<Config>::dd->makeZeroDensityOperator(static_cast<dd::Qubit>(qc->getNqubits()));
     Simulator<Config>::dd->incRef(rootEdge);
 
@@ -42,12 +42,11 @@ std::map<std::string, double> DeterministicNoiseSimulator<Config>::deterministic
 
         deterministicNoiseFunctionality.applyNoiseEffects(rootEdge, op);
     }
-    return Simulator<Config>::dd->getProbVectorFromDensityMatrix(rootEdge, measurementThreshold);
+    return rootEdge.getSparseProbabilityVectorStrKeys(measurementThreshold);
 }
 
 template<class Config>
-std::map<std::string, std::size_t> DeterministicNoiseSimulator<Config>::sampleFromProbabilityMap(const std::map<std::string, dd::fp>& resultProbabilityMap, std::size_t shots) {
-    // Create probability distribution from measure probabilities
+std::map<std::string, std::size_t> DeterministicNoiseSimulator<Config>::sampleFromProbabilityMap(const dd::SparsePVecStrKeys& resultProbabilityMap, std::size_t shots) {
     std::vector<dd::fp> weights;
     weights.reserve(resultProbabilityMap.size());
 

src/HybridSchrodingerFeynmanSimulator.cpp

@@ -244,11 +244,11 @@ void HybridSchrodingerFeynmanSimulator<Config>::simulateHybridAmplitudes(qc::Qub
     const auto requiredVectors     = static_cast<std::size_t>(std::ceil(static_cast<double>(maxControl) / static_cast<double>(nslicesOnOneCpu)));
     Simulator<Config>::rootEdge    = qc::VectorDD::zero;
 
-    std::vector<std::vector<std::complex<dd::fp>>> amplitudes(requiredVectors, std::vector<std::complex<dd::fp>>(1U << nqubits, {0, 0}));
+    std::vector<dd::CVec> amplitudes(requiredVectors, dd::CVec(1ULL << nqubits, 0));
 
     tf::Executor executor(actuallyUsedThreads);
     for (std::size_t control = 0, i = 0; control < maxControl; control += nslicesOnOneCpu, i++) {
-        executor.silent_async([this, i, &amplitudes, nslicesOnOneCpu, control, nqubits, splitQubit, maxControl]() {
+        executor.silent_async([this, i, &amplitudes, nslicesOnOneCpu, control, splitQubit, maxControl]() {
             const auto                         currentThread    = i;
             std::vector<std::complex<dd::fp>>& threadAmplitudes = amplitudes.at(currentThread);
 
@@ -259,7 +259,7 @@ void HybridSchrodingerFeynmanSimulator<Config>::simulateHybridAmplitudes(qc::Qub
                 }
                 std::unique_ptr<dd::Package<Config>> sliceDD = std::make_unique<dd::Package<Config>>(CircuitSimulator<Config>::getNumberOfQubits());
                 auto                                 result  = simulateSlicing(sliceDD, splitQubit, totalControl);
-                sliceDD->addAmplitudes(result, threadAmplitudes, static_cast<dd::Qubit>(nqubits));
+                result.addToVector(threadAmplitudes);
             }
         });
     }

src/ShorFastSimulator.cpp

@@ -63,7 +63,7 @@ std::map<std::string, std::size_t> ShorFastSimulator<Config>::simulate([[maybe_u
         uAEmulate2(as[i]);
 
         if (verbose) {
-            std::clog << "[ " << i + 1 << "/" << 2 * requiredBits << " ] QFT Pass. dd size=" << Simulator<Config>::dd->size(Simulator<Config>::rootEdge)
+            std::clog << "[ " << i + 1 << "/" << 2 * requiredBits << " ] QFT Pass. dd size=" << Simulator<Config>::rootEdge.size()
                       << "\n";
         }
 

src/ShorSimulator.cpp

@@ -68,13 +68,13 @@ std::map<std::string, std::size_t> ShorSimulator<Config>::simulate([[maybe_unuse
     }
 
     if (verbose) {
-        std::clog << "Nodes before QFT: " << Simulator<Config>::dd->size(Simulator<Config>::rootEdge) << "\n";
+        std::clog << "Nodes before QFT: " << Simulator<Config>::rootEdge.size() << "\n";
     }
 
     //EXACT QFT
     for (std::int32_t i = 0; i < static_cast<std::int32_t>(2 * requiredBits); i++) {
         if (verbose) {
-            std::clog << "[ " << i + 1 << "/" << 2 * requiredBits << " ] QFT Pass. dd size=" << Simulator<Config>::dd->size(Simulator<Config>::rootEdge)
+            std::clog << "[ " << i + 1 << "/" << 2 * requiredBits << " ] QFT Pass. dd size=" << Simulator<Config>::rootEdge.size()
                       << "\n";
         }
         double q = 2;

src/Simulator.cpp

@@ -28,30 +28,12 @@ std::map<std::string, std::size_t> Simulator<Config>::sampleFromAmplitudeVectorI
         auto m   = std::distance(amplitudes.begin(), mit);
 
         // construct basis state string
-        auto basisState = toBinaryString(static_cast<std::size_t>(m), getNumberOfQubits());
+        auto basisState = dd::intToBinaryString(static_cast<std::size_t>(m), getNumberOfQubits());
         results[basisState]++;
     }
     return results;
 }
 
-template<class Config>
-void Simulator<Config>::nextPath(std::string& s) {
-    std::string::reverse_iterator       iter = s.rbegin();
-    const std::string::reverse_iterator end  = s.rend();
-
-    int carry = 1;
-
-    while ((carry != 0) && iter != end) {
-        const int value = (*iter - '0') + carry;
-        carry           = (value / 2);
-        *iter           = static_cast<char>('0' + (value % 2));
-        ++iter;
-    }
-    if (carry != 0) {
-        s.insert(0, "1");
-    }
-}
-
 /**
  * Calculate the contributions of each node and return as vector of priority queues (each queue corresponds to a level in the decision diagram)
  * @tparam Config Configuration for the underlying DD package
@@ -169,8 +151,8 @@ double Simulator<Config>::approximateByFidelity(std::unique_ptr<dd::Package<Conf
     }
 
     if (verbose) {
-        const auto sizeBefore = localDD->size(edge);
-        const auto sizeAfter  = localDD->size(newEdge);
+        const auto sizeBefore = edge.size();
+        const auto sizeAfter  = newEdge.size();
         std::cout
                 << getName() << ","
                 << +getNumberOfQubits() << "," // unary plus for int promotion
@@ -258,8 +240,8 @@ double Simulator<Config>::approximateBySampling(std::unique_ptr<dd::Package<Conf
     }
 
     if (verbose) {
-        const auto sizeAfter  = localDD->size(newEdge);
-        const auto sizeBefore = localDD->size(edge);
+        const auto sizeAfter  = newEdge.size();
+        const auto sizeBefore = edge.size();
         std::cout
                 << getName() << ","
                 << +getNumberOfQubits() << "," // unary plus for int promotion

src/StochasticNoiseSimulator.cpp

@@ -243,10 +243,12 @@ void StochasticNoiseSimulator<Config>::runStochSimulationForId(std::size_t
                     controls = op->getControls();
 
                     if (targets.size() == 1 && controls.empty()) {
-                        operation = localDD->stochasticNoiseOperationCache.lookup(op->getType(), static_cast<dd::Qubit>(targets.front()));
-                        if (operation.p == nullptr) {
+                        auto* oper = localDD->stochasticNoiseOperationCache.lookup(op->getType(), static_cast<dd::Qubit>(targets.front()));
+                        if (oper == nullptr) {
                             operation = dd::getDD(op.get(), localDD);
                             localDD->stochasticNoiseOperationCache.insert(op->getType(), static_cast<dd::Qubit>(targets.front()), operation);
+                        } else {
+                            operation = *oper;
                         }
                     } else {
                         operation = dd::getDD(op.get(), localDD);
@@ -293,9 +295,8 @@ void StochasticNoiseSimulator<Config>::runStochSimulationForId(std::size_t
             } else {
                 // extract amplitude for state
                 const auto basisVector = recordedPropertiesList[i].second;
-                const auto amplitude   = localDD->getValueByPath(localRootEdge, basisVector);
-                const auto prob        = amplitude.r * amplitude.r + amplitude.i * amplitude.i;
-                recordedPropertiesStorage[i] += prob;
+                const auto amplitude   = localRootEdge.getValueByPath(basisVector);
+                recordedPropertiesStorage[i] += std::norm(amplitude);
             }
         }
     }