From 60d9a239e7742f694c4e5e1952321f05da81c9cc Mon Sep 17 00:00:00 2001 From: Daniel Haag <121057143+denialhaag@users.noreply.github.com> Date: Sun, 2 Aug 2026 10:21:45 +0200 Subject: [PATCH 1/5] =?UTF-8?q?=F0=9F=9A=9A=20Vendor=20density-matrix=20DD?= =?UTF-8?q?=20support?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Assisted-by: Claude Opus 4.8 via Claude Code --- CHANGELOG.md | 3 + include/DensityComputeTable.hpp | 144 +++++++ include/DensityDDPackage.hpp | 182 ++++++++ include/DensityNode.hpp | 330 +++++++++++++++ include/DensityNoiseTable.hpp | 118 ++++++ include/DensityUniqueTable.hpp | 135 ++++++ include/DeterministicNoiseSimulator.hpp | 41 +- include/NoiseFunctionality.hpp | 174 ++++++++ include/StochasticNoiseOperationTable.hpp | 90 ++++ include/StochasticNoiseSimulator.hpp | 22 +- src/DensityDDPackage.cpp | 488 ++++++++++++++++++++++ src/DensityNode.cpp | 400 ++++++++++++++++++ src/DensityUniqueTable.cpp | 163 ++++++++ src/DeterministicNoiseSimulator.cpp | 16 +- src/NoiseFunctionality.cpp | 485 +++++++++++++++++++++ src/StochasticNoiseSimulator.cpp | 15 +- 16 files changed, 2762 insertions(+), 44 deletions(-) create mode 100644 include/DensityComputeTable.hpp create mode 100644 include/DensityDDPackage.hpp create mode 100644 include/DensityNode.hpp create mode 100644 include/DensityNoiseTable.hpp create mode 100644 include/DensityUniqueTable.hpp create mode 100644 include/NoiseFunctionality.hpp create mode 100644 include/StochasticNoiseOperationTable.hpp create mode 100644 src/DensityDDPackage.cpp create mode 100644 src/DensityNode.cpp create mode 100644 src/DensityUniqueTable.cpp create mode 100644 src/NoiseFunctionality.cpp diff --git a/CHANGELOG.md b/CHANGELOG.md index 6d68e5360..8abd46364 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -12,6 +12,8 @@ releases may include breaking changes. ### Changed +- 🚚 Vendor density-matrix DD support used by the noise-aware simulators, as it + will be removed from `mqt-core` in version 4.0.0 ([#940]) ([**@denialhaag**]) - ⬆️ Update `mqt-core` to version 3.8.0 ([#939]) ([**@denialhaag**]) ## [2.4.0] - 2026-07-09 @@ -149,6 +151,7 @@ _📚 Refer to the [GitHub Release Notes] for previous changelogs._ +[#940]: https://github.com/munich-quantum-toolkit/ddsim/pull/940 [#939]: https://github.com/munich-quantum-toolkit/ddsim/pull/939 [#912]: https://github.com/munich-quantum-toolkit/ddsim/pull/912 [#911]: https://github.com/munich-quantum-toolkit/ddsim/pull/911 diff --git a/include/DensityComputeTable.hpp b/include/DensityComputeTable.hpp new file mode 100644 index 000000000..b093542f3 --- /dev/null +++ b/include/DensityComputeTable.hpp @@ -0,0 +1,144 @@ +/* + * Copyright (c) 2023 - 2026 Chair for Design Automation, TUM + * Copyright (c) 2025 - 2026 Munich Quantum Software Company GmbH + * All rights reserved. + * + * SPDX-License-Identifier: MIT + * + * Licensed under the MIT License + */ + +/** + * @file DensityComputeTable.hpp + * @brief Compute table for caching results of density-matrix operations. + * + * @details Self-contained copy of the MQT Core `dd::ComputeTable` that keeps + * the density-matrix specific hashing (accounting for the temporary flags + * encoded in the node pointer) and the `useDensityMatrix` discrimination on + * lookup. Both are required for the correctness of the reduced density-matrix + * representation. + */ + +#pragma once + +#include "DensityNode.hpp" +#include "dd/statistics/TableStatistics.hpp" +#include "ir/Definitions.hpp" + +#include +#include +#include +#include +#include +#include + +namespace dd::ddsim { + +/** + * @brief Data structure for caching computed results of binary operations on + * density-matrix DDs. + */ +template +class DensityComputeTable { +public: + static constexpr std::size_t DEFAULT_NUM_BUCKETS = 16384U; + + explicit DensityComputeTable( + const std::size_t numBuckets = DEFAULT_NUM_BUCKETS) { + if ((numBuckets & (numBuckets - 1)) != 0) { + throw std::invalid_argument("Number of buckets must be a power of two."); + } + stats.entrySize = sizeof(Entry); + stats.numBuckets = numBuckets; + valid = std::vector(numBuckets, false); + table = std::vector(numBuckets); + } + + struct Entry { + LeftOperandType leftOperand; + RightOperandType rightOperand; + ResultType result; + }; + + [[nodiscard]] std::size_t hash(const LeftOperandType& leftOperand, + const RightOperandType& rightOperand) const { + auto h1 = std::hash{}(leftOperand); + if constexpr (std::is_same_v) { + if (!dNode::isTerminal(leftOperand)) { + h1 = qc::combineHash( + h1, dNode::getDensityMatrixTempFlags(leftOperand->flags)); + } + } + auto h2 = std::hash{}(rightOperand); + if constexpr (std::is_same_v) { + if (!dNode::isTerminal(rightOperand)) { + h2 = qc::combineHash( + h2, dNode::getDensityMatrixTempFlags(rightOperand->flags)); + } + } + const auto hash = qc::combineHash(h1, h2); + const auto mask = stats.numBuckets - 1; + return hash & mask; + } + + [[nodiscard]] const auto& getTable() const { return table; } + [[nodiscard]] const auto& getStats() const noexcept { return stats; } + + void insert(const LeftOperandType& leftOperand, + const RightOperandType& rightOperand, const ResultType& result) { + const auto key = hash(leftOperand, rightOperand); + if (valid[key]) { + ++stats.collisions; + } else { + stats.trackInsert(); + valid[key] = true; + } + table[key] = {leftOperand, rightOperand, result}; + } + + ResultType* lookup(const LeftOperandType& leftOperand, + const RightOperandType& rightOperand, + [[maybe_unused]] const bool useDensityMatrix = false) { + ResultType* result = nullptr; + ++stats.lookups; + const auto key = hash(leftOperand, rightOperand); + if (!valid[key]) { + return result; + } + + auto& entry = table[key]; + if (entry.leftOperand != leftOperand) { + return result; + } + if (entry.rightOperand != rightOperand) { + return result; + } + + if constexpr (std::is_same_v || + std::is_same_v) { + // Since density matrices are reduced representations of matrices, a + // density matrix may not be returned when a matrix is required and vice + // versa + if (!dNode::isTerminal(entry.result.p) && + dNode::isDensityMatrixNode(entry.result.p->flags) != + useDensityMatrix) { + return result; + } + } + ++stats.hits; + return &entry.result; + } + + void clear() { valid = std::vector(stats.numBuckets, false); } + + std::ostream& printStatistics(std::ostream& os = std::cout) const { + return os << stats; + } + +private: + std::vector table; + std::vector valid; + dd::TableStatistics stats{}; +}; + +} // namespace dd::ddsim diff --git a/include/DensityDDPackage.hpp b/include/DensityDDPackage.hpp new file mode 100644 index 000000000..282dd9100 --- /dev/null +++ b/include/DensityDDPackage.hpp @@ -0,0 +1,182 @@ +/* + * Copyright (c) 2023 - 2026 Chair for Design Automation, TUM + * Copyright (c) 2025 - 2026 Munich Quantum Software Company GmbH + * All rights reserved. + * + * SPDX-License-Identifier: MIT + * + * Licensed under the MIT License + */ + +/** + * @file DensityDDPackage.hpp + * @brief Self-contained engine for density-matrix decision diagrams. + * + * @details This reconstructs the density-matrix operations that used to be + * member functions of the MQT Core `dd::Package` before they were removed in + * https://github.com/munich-quantum-toolkit/core/pull/1466. The engine owns the + * density-matrix node space (memory manager, unique table, compute tables) and + * borrows a `dd::Package` for the (surviving) matrix/complex-number + * infrastructure it needs (gate construction, conjugate transpose, complex + * number pool). + * + * The density DD reuses matrix nodes and complex numbers from the borrowed + * package (via the layout-compatible `densityFromMatrixEdge` reinterpretation). + * To keep those shared entries alive across the package's own garbage + * collection, every density root is additionally registered in the package's + * matrix root set (again via reinterpretation). The density node space itself + * is collected by this engine's own @ref garbageCollect. + */ + +#pragma once + +#include "DensityComputeTable.hpp" +#include "DensityNode.hpp" +#include "DensityUniqueTable.hpp" +#include "dd/ComplexValue.hpp" +#include "dd/DDDefinitions.hpp" +#include "dd/DDpackageConfig.hpp" +#include "dd/MemoryManager.hpp" +#include "dd/Node.hpp" +#include "dd/Package.hpp" +#include "dd/UnaryComputeTable.hpp" + +#include +#include +#include +#include +#include + +namespace dd::ddsim { + +/// Bucket sizes for the density-matrix compute and unique tables. +struct DensityDDPackageConfig { + std::size_t utDmNumBucket = 65536U; + std::size_t utDmInitialAllocationSize = 4096U; + std::size_t ctDmDmMultNumBucket = 16384U; + std::size_t ctDmAddNumBucket = 16384U; + std::size_t ctDmTraceNumBucket = 4096U; +}; + +/// Configuration for the (borrowed) matrix/vector `dd::Package` used by the +/// deterministic (density-matrix) noise-aware simulator. +constexpr auto DENSITY_MATRIX_SIMULATOR_DD_PACKAGE_CONFIG = []() { + dd::DDPackageConfig config{}; + config.utMatNumBucket = 16384U; + config.ctMatAddNumBucket = 4096U; + config.ctVecAddNumBucket = 4096U; + config.ctMatConjTransNumBucket = 4096U; + config.ctMatMatMultNumBucket = 1U; + config.ctMatVecMultNumBucket = 1U; + config.utVecNumBucket = 1U; + config.utVecInitialAllocationSize = 1U; + config.utMatInitialAllocationSize = 1U; + config.ctVecKronNumBucket = 1U; + config.ctMatKronNumBucket = 1U; + config.ctMatTraceNumBucket = 1U; + config.ctVecInnerProdNumBucket = 1U; + config.ctVecAddMagNumBucket = 1U; + config.ctMatAddMagNumBucket = 1U; + config.ctVecConjNumBucket = 1U; + return config; +}(); + +/// Configuration for the (borrowed) matrix/vector `dd::Package` used by the +/// stochastic noise-aware simulator. +constexpr auto STOCHASTIC_NOISE_SIMULATOR_DD_PACKAGE_CONFIG = []() { + dd::DDPackageConfig config{}; + config.ctVecAddMagNumBucket = 1U; + config.ctMatAddMagNumBucket = 1U; + config.ctVecConjNumBucket = 1U; + return config; +}(); + +/// Reinterpret a density-matrix edge as a matrix edge (inverse of +/// @ref densityFromMatrixEdge). Only safe for aligned (untagged) edges. +inline dd::mEdge matrixFromDensityEdge(const dEdge& e) { + return dd::mEdge{reinterpret_cast(e.p), e.w}; +} + +/// Engine that provides density-matrix DD operations on top of a `dd::Package`. +class DensityDDPackage { +public: + DensityDDPackage(dd::Package& package, std::size_t nqubits, + const DensityDDPackageConfig& config = {}) + : pkg(&package), dMemoryManager(dd::MemoryManager::create( + config.utDmInitialAllocationSize)), + dUniqueTable(dMemoryManager, + {.nVars = nqubits, .nBuckets = config.utDmNumBucket}), + densityAdd(config.ctDmAddNumBucket), + densityDensityMultiplication(config.ctDmDmMultNumBucket), + densityTrace(config.ctDmTraceNumBucket) {} + + /** + * @brief Construct the all-zero density operator \f$|0...0><0...0|\f$. + */ + dEdge makeZeroDensityOperator(std::size_t n); + + /** + * @brief Apply a matrix operation to a density matrix. + */ + dEdge applyOperationToDensity(dEdge& e, const dd::mEdge& operation); + + /** + * @brief Perform a collapsing measurement of a single qubit. + */ + char measureOneCollapsing(dEdge& e, dd::Qubit index, std::mt19937_64& mt); + + /** + * @brief Multiply two density-matrix DDs. + */ + dEdge multiply(const dEdge& x, const dEdge& y, + bool generateDensityMatrix = false); + + /** + * @brief Add two (cached) density-matrix DDs. + */ + dCachedEdge add2(const dCachedEdge& x, const dCachedEdge& y, dd::Qubit var); + + /** + * @brief Compute the trace of a density-matrix DD. + */ + dd::ComplexValue trace(const dEdge& a, std::size_t numQubits); + + /** + * @brief Create a normalized density-matrix node from a list of edges. + */ + dEdge makeDDNode(dd::Qubit var, const std::array& edges, + bool generateDensityMatrix = false); + dCachedEdge makeDDNode(dd::Qubit var, + const std::array& edges, + bool generateDensityMatrix = false); + + /// Increase the reference count of a density DD. + void incRef(const dEdge& e); + /// Decrease the reference count of a density DD. + void decRef(const dEdge& e); + + /// Trigger garbage collection of the density node space. + bool garbageCollect(bool force = false); + + /// Number of active density-matrix nodes. + [[nodiscard]] std::size_t computeActiveNodeCount() const; + + [[nodiscard]] dd::Package& package() const { return *pkg; } + +private: + dCachedEdge multiply2(const dEdge& x, const dEdge& y, dd::Qubit var, + bool generateDensityMatrix); + + dCachedEdge trace(const dEdge& a, const std::vector& eliminate, + std::size_t level, std::size_t alreadyEliminated = 0); + + dd::Package* pkg; + dd::MemoryManager dMemoryManager; + DensityUniqueTable dUniqueTable; + DensityComputeTable densityAdd; + DensityComputeTable densityDensityMultiplication; + dd::UnaryComputeTable densityTrace; + std::unordered_map dRoots; +}; + +} // namespace dd::ddsim diff --git a/include/DensityNode.hpp b/include/DensityNode.hpp new file mode 100644 index 000000000..70d7a8f00 --- /dev/null +++ b/include/DensityNode.hpp @@ -0,0 +1,330 @@ +/* + * Copyright (c) 2023 - 2026 Chair for Design Automation, TUM + * Copyright (c) 2025 - 2026 Munich Quantum Software Company GmbH + * All rights reserved. + * + * SPDX-License-Identifier: MIT + * + * Licensed under the MIT License + */ + +/** + * @file DensityNode.hpp + * @brief Density-matrix DD node, edge and cached-edge types. + * + * @details This is a self-contained copy of the density-matrix support that + * used to live in the MQT Core DD package (`dd::dNode`, `dd::dEdge`, + * `dd::dCachedEdge`) before it was removed in + * https://github.com/munich-quantum-toolkit/core/pull/1466. It lives in the + * separate namespace `dd::ddsim` so that it can coexist with a Core version + * that still provides the old types. + */ + +#pragma once + +#include "dd/Complex.hpp" +#include "dd/ComplexValue.hpp" +#include "dd/DDDefinitions.hpp" +#include "dd/Edge.hpp" +#include "dd/Node.hpp" + +#include +#include +#include +#include +#include + +namespace dd::ddsim { + +struct dNode; // NOLINT(readability-identifier-naming) + +/** + * @brief A weighted edge pointing to a density-matrix DD node. + */ +struct dEdge { // NOLINT(readability-identifier-naming) + dNode* p; + dd::Complex w; + + constexpr bool operator==(const dEdge& other) const { + return p == other.p && w.approximatelyEquals(other.w); + } + constexpr bool operator!=(const dEdge& other) const { + return !operator==(other); + } + + static constexpr dEdge zero() { return terminal(dd::Complex::zero()); } + static constexpr dEdge one() { return terminal(dd::Complex::one()); } + [[nodiscard]] static constexpr dEdge terminal(const dd::Complex& w); + + [[nodiscard]] static constexpr bool trackingRequired(const dEdge& e) { + return !e.isTerminal() || !dd::constants::isStaticNumber(e.w.r) || + !dd::constants::isStaticNumber(e.w.i); + } + + [[nodiscard]] bool isTerminal() const; + [[nodiscard]] bool isZeroTerminal() const { + return isTerminal() && w.exactlyZero(); + } + [[nodiscard]] bool isOneTerminal() const { + return isTerminal() && w.exactlyOne(); + } + + [[nodiscard]] bool isIdentity(bool upToGlobalPhase = true) const { + if (!isTerminal()) { + return false; + } + if (upToGlobalPhase) { + return !w.exactlyZero(); + } + return w.exactlyOne(); + } + + /** + * @brief Get the size of the DD (number of nodes including the terminal). + */ + [[nodiscard]] std::size_t size() const; + + /// @brief Mark the edge (and its sub-DD) as used. + void mark() const noexcept; + /// @brief Unmark the edge (and its sub-DD). + void unmark() const noexcept; + + /** + * @brief Get a normalized density-matrix DD from a fresh node and its edges. + */ + static auto normalize(dNode* p, const std::array& e, + dd::MemoryManager& mm, dd::ComplexNumbers& cn) -> dEdge; + + [[maybe_unused]] static void setDensityConjugateTrue(dEdge& e); + [[maybe_unused]] static void setFirstEdgeDensityPathTrue(dEdge& e); + static void setDensityMatrixTrue(dEdge& e); + static void alignDensityEdge(dEdge& e); + static void revertDmChangesToEdges(dEdge& x, dEdge& y); + static void revertDmChangesToEdge(dEdge& x); + static void applyDmChangesToEdges(dEdge& x, dEdge& y); + static void applyDmChangesToEdge(dEdge& x); + + /** + * @brief Get the sparse probability vector for the underlying density matrix. + */ + [[nodiscard]] dd::SparsePVec + getSparseProbabilityVector(std::size_t numQubits, + dd::fp threshold = 0.) const; + + /** + * @brief Get the sparse probability vector using strings as keys. + */ + [[nodiscard]] dd::SparsePVecStrKeys + getSparseProbabilityVectorStrKeys(std::size_t numQubits, + dd::fp threshold = 0.) const; + +private: + [[nodiscard]] std::size_t + size(std::unordered_set& visited) const; + + void traverseDiagonal(const dd::fp& prob, std::size_t i, + dd::ProbabilityFunc f, std::size_t level, + dd::fp threshold = 0.) const; +}; + +using DensityMatrixDD = dEdge; + +/** + * @brief A density-matrix DD node with a cached (non-canonical) edge weight. + */ +struct dCachedEdge { // NOLINT(readability-identifier-naming) + dNode* p{}; + dd::ComplexValue w; + + dCachedEdge() = default; + dCachedEdge(dNode* n, const dd::ComplexValue& v) : p(n), w(v) {} + dCachedEdge(dNode* n, const dd::Complex& c) + : p(n), w(static_cast(c)) {} + + bool operator==(const dCachedEdge& other) const { + return p == other.p && w.approximatelyEquals(other.w); + } + bool operator!=(const dCachedEdge& other) const { return !operator==(other); } + + [[nodiscard]] static dCachedEdge terminal(const dd::ComplexValue& w); + [[nodiscard]] static dCachedEdge terminal(const std::complex& w); + [[nodiscard]] static dCachedEdge terminal(const dd::Complex& w); + [[nodiscard]] static dCachedEdge zero() { + return terminal(dd::ComplexValue(0.)); + } + [[nodiscard]] static dCachedEdge one() { + return terminal(dd::ComplexValue(1.)); + } + + [[nodiscard]] bool isTerminal() const; + + [[nodiscard]] bool isIdentity(bool upToGlobalPhase = true) const { + if (!isTerminal()) { + return false; + } + if (upToGlobalPhase) { + return !w.exactlyZero(); + } + return w.exactlyOne(); + } + + /** + * @brief Get a normalized density-matrix DD from a fresh node and its edges. + */ + static auto normalize(dNode* p, const std::array& e, + dd::MemoryManager& mm, dd::ComplexNumbers& cn) + -> dCachedEdge; +}; + +/** + * @brief A density-matrix DD node. + * @details Data Layout (8)|(2|2|4)|(24|24|24|24) = 112B + */ +struct dNode final : dd::NodeBase { // NOLINT(readability-identifier-naming) + std::array e{}; // edges out of this node + + /// Getter for the next object + [[nodiscard]] dNode* next() const noexcept { + // NOLINTNEXTLINE(cppcoreguidelines-pro-type-static-cast-downcast) + return static_cast(next_); + } + /// Getter for the terminal object + static constexpr dNode* getTerminal() noexcept { return nullptr; } + + [[nodiscard]] [[maybe_unused]] static constexpr bool + tempDensityMatrixFlagsEqual(const std::uint8_t a, + const std::uint8_t b) noexcept { + return getDensityMatrixTempFlags(a) == getDensityMatrixTempFlags(b); + } + + [[nodiscard]] static constexpr bool + isConjugateTempFlagSet(const std::uintptr_t p) noexcept { + return (p & (1ULL << 0)) != 0U; + } + [[nodiscard]] static constexpr bool + isNonReduceTempFlagSet(const std::uintptr_t p) noexcept { + return (p & (1ULL << 1)) != 0U; + } + [[nodiscard]] static constexpr bool + isDensityMatrixTempFlagSet(const std::uintptr_t p) noexcept { + return (p & (1ULL << 2)) != 0U; + } + [[nodiscard]] static bool + isDensityMatrixNode(const std::uintptr_t p) noexcept { + return (p & (1ULL << 3)) != 0U; + } + + [[nodiscard]] static bool isConjugateTempFlagSet(const dNode* p) noexcept { + return isConjugateTempFlagSet(reinterpret_cast(p)); + } + [[nodiscard]] static bool isNonReduceTempFlagSet(const dNode* p) noexcept { + return isNonReduceTempFlagSet(reinterpret_cast(p)); + } + [[nodiscard]] static bool + isDensityMatrixTempFlagSet(const dNode* p) noexcept { + return isDensityMatrixTempFlagSet(reinterpret_cast(p)); + } + [[nodiscard]] static bool isDensityMatrixNode(const dNode* p) noexcept { + return isDensityMatrixNode(reinterpret_cast(p)); + } + + static void setConjugateTempFlagTrue(dNode*& p) noexcept { + p = reinterpret_cast(reinterpret_cast(p) | + (1ULL << 0)); + } + static void setNonReduceTempFlagTrue(dNode*& p) noexcept { + p = reinterpret_cast(reinterpret_cast(p) | + (1ULL << 1)); + } + static void setDensityMatTempFlagTrue(dNode*& p) noexcept { + p = reinterpret_cast(reinterpret_cast(p) | + (1ULL << 2)); + } + static void alignDensityNode(dNode*& p) noexcept { + p = reinterpret_cast(reinterpret_cast(p) & (~7ULL)); + } + + [[nodiscard]] static std::uintptr_t + getDensityMatrixTempFlags(dNode*& p) noexcept { + return getDensityMatrixTempFlags(reinterpret_cast(p)); + } + [[nodiscard]] static constexpr std::uintptr_t + getDensityMatrixTempFlags(const std::uintptr_t a) noexcept { + return a & (7ULL); + } + + constexpr void unsetTempDensityMatrixFlags() noexcept { + flags = flags & static_cast(~7U); + } + + void setDensityMatrixNodeFlag(bool densityMatrix) noexcept; + + static std::uint8_t alignDensityNodeNode(dNode*& p) noexcept; + + static void getAlignedNodeRevertModificationsOnSubEdges(dNode* p) noexcept; + + static void applyDmChangesToNode(dNode*& p) noexcept; + + static void revertDmChangesToNode(dNode*& p) noexcept; +}; + +/// Reinterpret a Core matrix edge as a density-matrix edge. +inline dEdge densityFromMatrixEdge(const dd::mEdge& e) { + return dEdge{reinterpret_cast(e.p), e.w}; +} + +///----------------------------------------------------------------------------- +/// Inline definitions that require the complete `dNode` type +///----------------------------------------------------------------------------- + +constexpr dEdge dEdge::terminal(const dd::Complex& w) { + return dEdge{dNode::getTerminal(), w}; +} + +inline bool dEdge::isTerminal() const { return dNode::isTerminal(p); } + +inline void dEdge::setDensityConjugateTrue(dEdge& e) { + dNode::setConjugateTempFlagTrue(e.p); +} +inline void dEdge::setFirstEdgeDensityPathTrue(dEdge& e) { + dNode::setNonReduceTempFlagTrue(e.p); +} +inline void dEdge::setDensityMatrixTrue(dEdge& e) { + dNode::setDensityMatTempFlagTrue(e.p); +} +inline void dEdge::alignDensityEdge(dEdge& e) { dNode::alignDensityNode(e.p); } +inline void dEdge::revertDmChangesToEdges(dEdge& x, dEdge& y) { + revertDmChangesToEdge(x); + revertDmChangesToEdge(y); +} +inline void dEdge::revertDmChangesToEdge(dEdge& x) { + dNode::revertDmChangesToNode(x.p); +} +inline void dEdge::applyDmChangesToEdges(dEdge& x, dEdge& y) { + applyDmChangesToEdge(x); + applyDmChangesToEdge(y); +} +inline void dEdge::applyDmChangesToEdge(dEdge& x) { + dNode::applyDmChangesToNode(x.p); +} + +inline dCachedEdge dCachedEdge::terminal(const dd::ComplexValue& w) { + return dCachedEdge{dNode::getTerminal(), w}; +} +inline dCachedEdge dCachedEdge::terminal(const std::complex& w) { + return dCachedEdge{dNode::getTerminal(), static_cast(w)}; +} +inline dCachedEdge dCachedEdge::terminal(const dd::Complex& w) { + return terminal(static_cast(w)); +} +inline bool dCachedEdge::isTerminal() const { return dNode::isTerminal(p); } + +} // namespace dd::ddsim + +template <> struct std::hash { + std::size_t operator()(const dd::ddsim::dEdge& e) const noexcept; +}; + +template <> struct std::hash { + std::size_t operator()(const dd::ddsim::dCachedEdge& e) const noexcept; +}; diff --git a/include/DensityNoiseTable.hpp b/include/DensityNoiseTable.hpp new file mode 100644 index 000000000..56a59b640 --- /dev/null +++ b/include/DensityNoiseTable.hpp @@ -0,0 +1,118 @@ +/* + * Copyright (c) 2023 - 2026 Chair for Design Automation, TUM + * Copyright (c) 2025 - 2026 Munich Quantum Software Company GmbH + * All rights reserved. + * + * SPDX-License-Identifier: MIT + * + * Licensed under the MIT License + */ + +/** + * @file DensityNoiseTable.hpp + * @brief Data structure for caching computed results of noise operations + * + * @details Self-contained copy of the MQT Core `dd::DensityNoiseTable` that was + * removed alongside the density-matrix support. + */ + +#pragma once + +#include "dd/DDDefinitions.hpp" +#include "dd/statistics/TableStatistics.hpp" + +#include +#include +#include +#include + +namespace dd::ddsim { + +/** + * @brief Data structure for caching computed results of noise operations + * @tparam OperandType type of the operation's operand + * @tparam ResultType type of the operation's result + */ +template +class DensityNoiseTable { // todo: Inherit from UnaryComputerTable +public: + /** + * Default constructor + * @param numBuckets Number of hash table buckets. Must be a power of two. + */ + explicit DensityNoiseTable(const std::size_t numBuckets = 32768U) { + // numBuckets must be a power of two + if ((numBuckets & (numBuckets - 1)) != 0) { + throw std::invalid_argument("Number of buckets must be a power of two."); + } + stats.entrySize = sizeof(Entry); + stats.numBuckets = numBuckets; + valid = std::vector(numBuckets, false); + table = std::vector(numBuckets); + } + + struct Entry { + OperandType operand; + ResultType result; + std::vector usedQubits; + }; + + /// Get a reference to the table + [[nodiscard]] const auto& getTable() const { return table; } + + /// Get a reference to the statistics + [[nodiscard]] const auto& getStats() const noexcept { return stats; } + + [[nodiscard]] std::size_t + hash(const OperandType& a, const std::vector& usedQubits) const { + std::size_t i = 0; + for (const auto qubit : usedQubits) { + i = (i << 3U) + (i * static_cast(qubit)) + + static_cast(qubit); + } + const std::size_t mask = stats.numBuckets - 1; + return (std::hash{}(a) + i) & mask; + } + + void insert(const OperandType& operand, const ResultType& result, + const std::vector& usedQubits) { + const auto key = hash(operand, usedQubits); + if (valid[key]) { + ++stats.collisions; + } else { + stats.trackInsert(); + valid[key] = true; + } + table[key] = {operand, result, usedQubits}; + } + + ResultType lookup(const OperandType& operand, + const std::vector& usedQubits) { + ResultType result{}; + ++stats.lookups; + const auto key = hash(operand, usedQubits); + + if (!valid[key]) { + return result; + } + + auto& entry = table[key]; + if (entry.operand != operand) { + return result; + } + if (entry.usedQubits != usedQubits) { + return result; + } + ++stats.hits; + return entry.result; + } + + void clear() { valid = std::vector(stats.numBuckets, false); } + +private: + std::vector table; + std::vector valid; + dd::TableStatistics stats{}; +}; + +} // namespace dd::ddsim diff --git a/include/DensityUniqueTable.hpp b/include/DensityUniqueTable.hpp new file mode 100644 index 000000000..6e7627353 --- /dev/null +++ b/include/DensityUniqueTable.hpp @@ -0,0 +1,135 @@ +/* + * Copyright (c) 2023 - 2026 Chair for Design Automation, TUM + * Copyright (c) 2025 - 2026 Munich Quantum Software Company GmbH + * All rights reserved. + * + * SPDX-License-Identifier: MIT + * + * Licensed under the MIT License + */ + +/** + * @file DensityUniqueTable.hpp + * @brief Unique table for density-matrix DD nodes. + * + * @details Self-contained copy of the relevant parts of the MQT Core + * `dd::UniqueTable`, specialized for the density-matrix node type + * `dd::ddsim::dNode`. Unlike the (post-removal) Core unique table, the node + * equality check accounts for the persistent density-matrix flag, which is + * required for the correctness of the reduced density-matrix representation. + */ + +#pragma once + +#include "DensityNode.hpp" +#include "dd/MemoryManager.hpp" +#include "dd/Node.hpp" +#include "dd/statistics/UniqueTableStatistics.hpp" +#include "ir/Definitions.hpp" + +#include +#include +#include +#include + +namespace dd::ddsim { + +/// Data structure for uniquely storing density-matrix DD nodes. +class DensityUniqueTable { +public: + static constexpr std::size_t INITIAL_GC_LIMIT = 131072U; + + struct UniqueTableConfig { + std::size_t nVars = 0U; + std::size_t nBuckets = 32768; + std::size_t initialGCLimit = INITIAL_GC_LIMIT; + }; + + DensityUniqueTable(dd::MemoryManager& manager, + const UniqueTableConfig& config); + + void resize(std::size_t nVars); + + [[nodiscard]] std::size_t hash(const dNode& p) const { + const std::size_t mask = cfg.nBuckets - 1; + std::size_t key = 0U; + for (const auto& succ : p.e) { + qc::hashCombine(key, std::hash{}(succ)); + } + key &= mask; + return key; + } + + [[nodiscard]] static bool nodesAreEqual(const dNode& p, const dNode& q) { + return (p.e == q.e && (p.flags == q.flags)); + } + + // Lookup a node in the unique table and insert it if it has not been found. + // Only normalized nodes shall be stored. + [[nodiscard]] dNode* lookup(dNode* p) { + // there are unique terminal nodes + if (dd::NodeBase::isTerminal(p)) { + return p; + } + + const auto key = hash(*p); + const auto v = p->v; + ++stats[v].lookups; + + if (auto* hashedNode = searchTable(*p, key); + !dNode::isTerminal(hashedNode)) { + return hashedNode; + } + + p->setNext(tables[v][key]); + tables[v][key] = p; + stats[v].trackInsert(); + + return p; + } + + [[nodiscard]] const auto& getTables() const { return tables; } + [[nodiscard]] const auto& getStats() const noexcept { return stats; } + [[nodiscard]] const dd::UniqueTableStatistics& + getStats(std::size_t idx) const noexcept; + [[nodiscard]] nlohmann::basic_json<> + getStatsJson(bool includeIndividualTables = false) const; + [[nodiscard]] std::size_t getNumEntries() const noexcept; + [[nodiscard]] std::size_t countMarkedEntries() const noexcept; + [[nodiscard]] bool possiblyNeedsCollection() const; + std::size_t garbageCollect(bool force = false); + void clear(); + +private: + using Bucket = dd::NodeBase*; + using Table = std::vector; + + UniqueTableConfig cfg; + std::size_t gcLimit; + dd::MemoryManager* memoryManager; + std::vector tables; + std::vector stats; + + [[nodiscard]] dNode* searchTable(dNode& p, const std::size_t& key) { + const auto v = p.v; + auto* bucket = static_cast(tables[v][key]); + while (bucket != nullptr) { + if (nodesAreEqual(p, *bucket)) { + // Match found + if (&p != bucket) { + // put node pointed to by p on available chain + memoryManager->returnEntry(p); + } + ++stats[v].hits; + return bucket; + } + ++stats[v].collisions; + bucket = bucket->next(); + } + + // Node not found in bucket + return dNode::getTerminal(); + } +}; + +} // namespace dd::ddsim diff --git a/include/DeterministicNoiseSimulator.hpp b/include/DeterministicNoiseSimulator.hpp index b1db3b593..ba7dfa585 100644 --- a/include/DeterministicNoiseSimulator.hpp +++ b/include/DeterministicNoiseSimulator.hpp @@ -11,11 +11,11 @@ #pragma once #include "CircuitSimulator.hpp" +#include "DensityDDPackage.hpp" +#include "DensityNode.hpp" +#include "NoiseFunctionality.hpp" #include "Simulator.hpp" #include "dd/DDDefinitions.hpp" -#include "dd/DDpackageConfig.hpp" -#include "dd/Node.hpp" -#include "dd/NoiseFunctionality.hpp" #include "dd/Package.hpp" #include "ir/QuantumComputation.hpp" #include "ir/operations/NonUnitaryOperation.hpp" @@ -37,7 +37,7 @@ class DeterministicNoiseSimulator : public CircuitSimulator { std::optional ampDampingProbability_ = std::nullopt, double multiQubitGateFactor_ = 2) : CircuitSimulator(std::move(qc_), approximationInfo_, - dd::DENSITY_MATRIX_SIMULATOR_DD_PACKAGE_CONFIG), + dd::ddsim::DENSITY_MATRIX_SIMULATOR_DD_PACKAGE_CONFIG), noiseEffects(std::move(noiseEffects_)), noiseProbSingleQubit(noiseProbability_), ampDampingProbSingleQubit(ampDampingProbability_ @@ -46,11 +46,12 @@ class DeterministicNoiseSimulator : public CircuitSimulator { noiseProbMultiQubit(noiseProbability_ * multiQubitGateFactor_), ampDampingProbMultiQubit(ampDampingProbSingleQubit * multiQubitGateFactor_), + densityDD(*dd, CircuitSimulator::getNumberOfQubits()), deterministicNoiseFunctionality( - *dd, CircuitSimulator::getNumberOfQubits(), noiseProbSingleQubit, - noiseProbMultiQubit, ampDampingProbSingleQubit, - ampDampingProbMultiQubit, noiseEffects) { - dd::sanityCheckOfNoiseProbabilities( + densityDD, CircuitSimulator::getNumberOfQubits(), + noiseProbSingleQubit, noiseProbMultiQubit, + ampDampingProbSingleQubit, ampDampingProbMultiQubit, noiseEffects) { + dd::ddsim::sanityCheckOfNoiseProbabilities( noiseProbability_, ampDampingProbSingleQubit, multiQubitGateFactor_); } @@ -71,7 +72,7 @@ class DeterministicNoiseSimulator : public CircuitSimulator { std::optional ampDampingProbability_ = std::nullopt, double multiQubitGateFactor_ = 2) : CircuitSimulator(std::move(qc_), approximationInfo_, seed_, - dd::DENSITY_MATRIX_SIMULATOR_DD_PACKAGE_CONFIG), + dd::ddsim::DENSITY_MATRIX_SIMULATOR_DD_PACKAGE_CONFIG), noiseEffects(std::move(noiseEffects_)), noiseProbSingleQubit(noiseProbability_), ampDampingProbSingleQubit(ampDampingProbability_ @@ -80,11 +81,12 @@ class DeterministicNoiseSimulator : public CircuitSimulator { noiseProbMultiQubit(noiseProbability_ * multiQubitGateFactor_), ampDampingProbMultiQubit(ampDampingProbSingleQubit * multiQubitGateFactor_), + densityDD(*dd, CircuitSimulator::getNumberOfQubits()), deterministicNoiseFunctionality( - *dd, CircuitSimulator::getNumberOfQubits(), noiseProbSingleQubit, - noiseProbMultiQubit, ampDampingProbSingleQubit, - ampDampingProbMultiQubit, noiseEffects) { - dd::sanityCheckOfNoiseProbabilities( + densityDD, CircuitSimulator::getNumberOfQubits(), + noiseProbSingleQubit, noiseProbMultiQubit, + ampDampingProbSingleQubit, ampDampingProbMultiQubit, noiseEffects) { + dd::ddsim::sanityCheckOfNoiseProbabilities( noiseProbability_, ampDampingProbSingleQubit, multiQubitGateFactor_); } @@ -106,19 +108,17 @@ class DeterministicNoiseSimulator : public CircuitSimulator { std::size_t shots); [[nodiscard]] std::size_t getActiveNodeCount() const override { - const auto [vectorNodes, matrixNodes, densityNodes, realNumbers] = - dd->computeActiveCounts(); - return densityNodes; + return densityDD.computeActiveNodeCount(); } [[nodiscard]] std::size_t countNodesFromRoot() override { - dd::DensityMatrixDD::alignDensityEdge(rootEdge); + dd::ddsim::DensityMatrixDD::alignDensityEdge(rootEdge); const std::size_t tmp = rootEdge.size(); - dd::DensityMatrixDD::setDensityMatrixTrue(rootEdge); + dd::ddsim::DensityMatrixDD::setDensityMatrixTrue(rootEdge); return tmp; } - dd::DensityMatrixDD rootEdge{}; + dd::ddsim::DensityMatrixDD rootEdge{}; private: std::string noiseEffects; @@ -129,5 +129,6 @@ class DeterministicNoiseSimulator : public CircuitSimulator { double ampDampingProbMultiQubit{}; double measurementThreshold = 0.01; - dd::DeterministicNoiseFunctionality deterministicNoiseFunctionality; + dd::ddsim::DensityDDPackage densityDD; + dd::ddsim::DeterministicNoiseFunctionality deterministicNoiseFunctionality; }; diff --git a/include/NoiseFunctionality.hpp b/include/NoiseFunctionality.hpp new file mode 100644 index 000000000..3f89b6079 --- /dev/null +++ b/include/NoiseFunctionality.hpp @@ -0,0 +1,174 @@ +/* + * Copyright (c) 2023 - 2026 Chair for Design Automation, TUM + * Copyright (c) 2025 - 2026 Munich Quantum Software Company GmbH + * All rights reserved. + * + * SPDX-License-Identifier: MIT + * + * Licensed under the MIT License + */ + +/** + * @file NoiseFunctionality.hpp + * @brief Stochastic and deterministic noise functionality. + * + * @details Self-contained copy of the MQT Core `dd::NoiseFunctionality` that + * was removed alongside the density-matrix support. The stochastic noise + * operation cache is now owned by @ref StochasticNoiseFunctionality (instead of + * the DD package), and the removed noise `qc::OpType`s + * (`ATrue`/`AFalse`/`MultiATrue`/`MultiAFalse`) are replaced by a local + * enumeration used only for cache keys and noise-operation selection. + */ + +#pragma once + +#include "DensityDDPackage.hpp" +#include "DensityNode.hpp" +#include "StochasticNoiseOperationTable.hpp" +#include "dd/DDDefinitions.hpp" +#include "dd/Package.hpp" +#include "ir/Definitions.hpp" +#include "ir/operations/Operation.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace dd::ddsim { + +using NrEdges = std::tuple_size; +using ArrayOfEdges = std::array; + +// noise operations available for deterministic noise aware quantum circuit +// simulation +enum NoiseOperations : std::uint8_t { + AmplitudeDamping, + PhaseFlip, + Depolarization, + Identity +}; + +void sanityCheckOfNoiseProbabilities(double noiseProbability, + double amplitudeDampingProb, + double multiQubitGateFactor); + +class StochasticNoiseFunctionality { +public: + StochasticNoiseFunctionality(dd::Package& dd, std::size_t nq, + double gateNoiseProbability, + double amplitudeDampingProb, + double multiQubitGateFactor, + const std::string& cNoiseEffects); + + ~StochasticNoiseFunctionality() { package->decRef(identityDD); } + +protected: + /// Local identifiers for cached stochastic noise operations, replacing the + /// removed noise `qc::OpType`s. + enum StochasticNoiseKind : std::uint8_t { + StochX, + StochY, + StochZ, + StochATrue, + StochAFalse, + StochMultiATrue, + StochMultiAFalse, + StochIdentity, + StochasticNoiseKindEnd + }; + + dd::Package* package; + std::size_t nQubits; + std::uniform_real_distribution dist; + + double noiseProbability; + double noiseProbabilityMulti; + dd::fp sqrtAmplitudeDampingProbability; + dd::fp oneMinusSqrtAmplitudeDampingProbability; + dd::fp sqrtAmplitudeDampingProbabilityMulti; + dd::fp oneMinusSqrtAmplitudeDampingProbabilityMulti; + dd::GateMatrix ampDampingTrue{}; + dd::GateMatrix ampDampingTrueMulti{}; + dd::GateMatrix ampDampingFalse{}; + dd::GateMatrix ampDampingFalseMulti{}; + std::vector noiseEffects; + dd::mEdge identityDD; + StochasticNoiseOperationTable stochasticNoiseOperationCache; + + [[nodiscard]] std::size_t getNumberOfQubits() const { return nQubits; } + [[nodiscard]] double getNoiseProbability(bool multiQubitNoiseFlag) const; + + [[nodiscard]] static StochasticNoiseKind + getAmplitudeDampingOperationType(bool multiQubitNoiseFlag, + bool amplitudeDampingFlag); + + [[nodiscard]] dd::GateMatrix + getAmplitudeDampingOperationMatrix(bool multiQubitNoiseFlag, + bool amplitudeDampingFlag) const; + +public: + void applyNoiseOperation(const std::set& targets, + dd::mEdge operation, dd::vEdge& state, + std::mt19937_64& generator); + +protected: + [[nodiscard]] dd::mEdge stackOperation(const dd::mEdge& operation, + qc::Qubit target, + StochasticNoiseKind noiseOperation, + const dd::GateMatrix& matrix); + + dd::mEdge generateNoiseOperation(dd::mEdge operation, qc::Qubit target, + std::mt19937_64& generator, + bool amplitudeDamping, + bool multiQubitOperation); + + [[nodiscard]] StochasticNoiseKind + returnNoiseOperation(NoiseOperations noiseOperation, double prob, + bool multiQubitNoiseFlag) const; +}; + +class DeterministicNoiseFunctionality { +public: + DeterministicNoiseFunctionality(DensityDDPackage& dd, std::size_t nq, + double noiseProbabilitySingleQubit, + double noiseProbabilityMultiQubit, + double ampDampProbSingleQubit, + double ampDampProbMultiQubit, + const std::string& cNoiseEffects); + +protected: + DensityDDPackage* package; + std::size_t nQubits; + + double noiseProbSingleQubit; + double noiseProbMultiQubit; + double ampDampingProbSingleQubit; + double ampDampingProbMultiQubit; + + std::vector noiseEffects; + + [[nodiscard]] std::size_t getNumberOfQubits() const { return nQubits; } + +public: + void applyNoiseEffects(dEdge& originalEdge, + const std::unique_ptr& qcOperation); + +private: + dCachedEdge applyNoiseEffects(dEdge& originalEdge, + const std::set& usedQubits, + bool firstPathEdge, dd::Qubit level); + + static void applyPhaseFlipToEdges(ArrayOfEdges& e, double probability); + + void applyAmplitudeDampingToEdges(ArrayOfEdges& e, double probability) const; + + void applyDepolarisationToEdges(ArrayOfEdges& e, double probability) const; +}; + +} // namespace dd::ddsim diff --git a/include/StochasticNoiseOperationTable.hpp b/include/StochasticNoiseOperationTable.hpp new file mode 100644 index 000000000..7ff41025e --- /dev/null +++ b/include/StochasticNoiseOperationTable.hpp @@ -0,0 +1,90 @@ +/* + * Copyright (c) 2023 - 2026 Chair for Design Automation, TUM + * Copyright (c) 2025 - 2026 Munich Quantum Software Company GmbH + * All rights reserved. + * + * SPDX-License-Identifier: MIT + * + * Licensed under the MIT License + */ + +/** + * @file StochasticNoiseOperationTable.hpp + * @brief Data structure for caching computed results of stochastic operations + * + * @details Self-contained copy of the MQT Core + * `dd::StochasticNoiseOperationTable` that was removed alongside the + * density-matrix support. The number of cached operations is provided + * explicitly instead of being derived from the (removed) noise `OpType`s. + */ + +#pragma once + +#include "dd/statistics/TableStatistics.hpp" +#include "ir/Definitions.hpp" + +#include +#include +#include +#include +#include + +namespace dd::ddsim { + +template class StochasticNoiseOperationTable { +public: + StochasticNoiseOperationTable(const std::size_t nv, + const std::size_t numberOfStochasticOperations) + : nvars(nv), numberOfOperations(numberOfStochasticOperations), + table(nv, std::vector(numberOfStochasticOperations)) { + stats.entrySize = sizeof(Edge); + stats.numBuckets = nv * numberOfStochasticOperations; + } + + /// Get a reference to the table + [[nodiscard]] const auto& getTable() const { return table; } + + /// Get a reference to the statistics + [[nodiscard]] const auto& getStats() const noexcept { return stats; } + + void resize(const std::size_t nq) { + nvars = nq; + table.resize(nvars, std::vector(numberOfOperations)); + } + + void insert(std::uint8_t kind, qc::Qubit target, const Edge& r) { + // Increase numberOfOperations if this assertion is hit for a valid kind. + assert(kind < numberOfOperations); + table.at(target).at(kind) = r; + stats.trackInsert(); + } + + Edge* lookup(std::uint8_t kind, qc::Qubit target) { + // Increase numberOfOperations if this assertion is hit for a valid kind. + assert(kind < numberOfOperations); + ++stats.lookups; + auto& entry = table.at(target).at(kind); + if (entry.w.r == nullptr) { + return nullptr; + } + ++stats.hits; + return &entry; + } + + void clear() { + if (stats.numEntries > 0) { + for (auto& t : table) { + std::fill(t.begin(), t.end(), Edge{}); + } + stats.numEntries = 0; + } + } + +private: + std::size_t nvars; + std::size_t numberOfOperations; + std::vector> table; + dd::TableStatistics stats{}; +}; + +} // namespace dd::ddsim diff --git a/include/StochasticNoiseSimulator.hpp b/include/StochasticNoiseSimulator.hpp index 58847ceb3..1b72659f6 100644 --- a/include/StochasticNoiseSimulator.hpp +++ b/include/StochasticNoiseSimulator.hpp @@ -11,8 +11,8 @@ #pragma once #include "CircuitSimulator.hpp" -#include "dd/DDpackageConfig.hpp" -#include "dd/NoiseFunctionality.hpp" +#include "DensityDDPackage.hpp" +#include "NoiseFunctionality.hpp" #include "ir/Definitions.hpp" #include "ir/QuantumComputation.hpp" @@ -34,8 +34,9 @@ class StochasticNoiseSimulator final : public CircuitSimulator { std::string noiseEffects_ = "APD", double noiseProbability_ = 0.001, std::optional ampDampingProbability_ = std::nullopt, double multiQubitGateFactor_ = 2) - : CircuitSimulator(std::move(qc_), approximationInfo_, - dd::STOCHASTIC_NOISE_SIMULATOR_DD_PACKAGE_CONFIG), + : CircuitSimulator( + std::move(qc_), approximationInfo_, + dd::ddsim::STOCHASTIC_NOISE_SIMULATOR_DD_PACKAGE_CONFIG), noiseProbability(noiseProbability_), amplitudeDampingProb((ampDampingProbability_) ? ampDampingProbability_.value() @@ -45,8 +46,8 @@ class StochasticNoiseSimulator final : public CircuitSimulator { ? std::thread::hardware_concurrency() - 4 : 1), noiseEffects(std::move(noiseEffects_)) { - dd::sanityCheckOfNoiseProbabilities(noiseProbability, amplitudeDampingProb, - multiQubitGateFactor); + dd::ddsim::sanityCheckOfNoiseProbabilities( + noiseProbability, amplitudeDampingProb, multiQubitGateFactor); } explicit StochasticNoiseSimulator( @@ -64,8 +65,9 @@ class StochasticNoiseSimulator final : public CircuitSimulator { std::string noiseEffects_ = "APD", double noiseProbability_ = 0.001, std::optional ampDampingProbability_ = std::nullopt, double multiQubitGateFactor_ = 2) - : CircuitSimulator(std::move(qc_), approximationInfo_, seed_, - dd::STOCHASTIC_NOISE_SIMULATOR_DD_PACKAGE_CONFIG), + : CircuitSimulator( + std::move(qc_), approximationInfo_, seed_, + dd::ddsim::STOCHASTIC_NOISE_SIMULATOR_DD_PACKAGE_CONFIG), noiseProbability(noiseProbability_), amplitudeDampingProb((ampDampingProbability_) ? ampDampingProbability_.value() @@ -75,8 +77,8 @@ class StochasticNoiseSimulator final : public CircuitSimulator { ? std::thread::hardware_concurrency() - 4 : 1), noiseEffects(std::move(noiseEffects_)) { - dd::sanityCheckOfNoiseProbabilities(noiseProbability, amplitudeDampingProb, - multiQubitGateFactor); + dd::ddsim::sanityCheckOfNoiseProbabilities( + noiseProbability, amplitudeDampingProb, multiQubitGateFactor); } std::vector> classicalMeasurementsMaps; diff --git a/src/DensityDDPackage.cpp b/src/DensityDDPackage.cpp new file mode 100644 index 000000000..b4c9bc21a --- /dev/null +++ b/src/DensityDDPackage.cpp @@ -0,0 +1,488 @@ +/* + * Copyright (c) 2023 - 2026 Chair for Design Automation, TUM + * Copyright (c) 2025 - 2026 Munich Quantum Software Company GmbH + * All rights reserved. + * + * SPDX-License-Identifier: MIT + * + * Licensed under the MIT License + */ + +#include "DensityDDPackage.hpp" + +#include "DensityNode.hpp" +#include "dd/Complex.hpp" +#include "dd/ComplexNumbers.hpp" +#include "dd/ComplexValue.hpp" +#include "dd/DDDefinitions.hpp" +#include "dd/GateMatrixDefinitions.hpp" +#include "dd/Node.hpp" +#include "dd/Package.hpp" + +#include +#include +#include +#include +#include + +namespace dd::ddsim { + +///----------------------------------------------------------------------------- +/// \n Node creation \n +///----------------------------------------------------------------------------- + +dEdge DensityDDPackage::makeDDNode(const dd::Qubit var, + const std::array& edges, + const bool generateDensityMatrix) { + auto& mm = dMemoryManager; + auto* p = mm.get(); + p->v = var; + p->flags = 0; + p->setDensityMatrixNodeFlag(generateDensityMatrix); + + auto e = dEdge::normalize(p, edges, mm, pkg->cn); + if (!e.isTerminal()) { + const auto& es = e.p->e; + // Check if node resembles the identity. If so, skip it. + if ((es[0].p == es[3].p) && + (es[0].w.exactlyOne() && es[1].w.exactlyZero() && + es[2].w.exactlyZero() && es[3].w.exactlyOne())) { + auto* ptr = es[0].p; + mm.returnEntry(*e.p); + return dEdge{ptr, e.w}; + } + } + + auto* l = dUniqueTable.lookup(e.p); + return dEdge{l, e.w}; +} + +dCachedEdge +DensityDDPackage::makeDDNode(const dd::Qubit var, + const std::array& edges, + const bool generateDensityMatrix) { + auto& mm = dMemoryManager; + auto* p = mm.get(); + p->v = var; + p->flags = 0; + p->setDensityMatrixNodeFlag(generateDensityMatrix); + + auto e = dCachedEdge::normalize(p, edges, mm, pkg->cn); + if (!e.isTerminal()) { + const auto& es = e.p->e; + if ((es[0].p == es[3].p) && + (es[0].w.exactlyOne() && es[1].w.exactlyZero() && + es[2].w.exactlyZero() && es[3].w.exactlyOne())) { + auto* ptr = es[0].p; + mm.returnEntry(*e.p); + return dCachedEdge{ptr, e.w}; + } + } + + auto* l = dUniqueTable.lookup(e.p); + return dCachedEdge{l, e.w}; +} + +///----------------------------------------------------------------------------- +/// \n Addition \n +///----------------------------------------------------------------------------- + +dCachedEdge DensityDDPackage::add2(const dCachedEdge& x, const dCachedEdge& y, + const dd::Qubit var) { + if (x.w.exactlyZero()) { + if (y.w.exactlyZero()) { + return dCachedEdge::zero(); + } + return y; + } + if (y.w.exactlyZero()) { + return x; + } + if (x.p == y.p) { + const auto rWeight = x.w + y.w; + return {x.p, rWeight}; + } + + if (const auto* r = densityAdd.lookup(x, y); r != nullptr) { + return *r; + } + + constexpr std::size_t n = dd::NEDGE; + std::array edge{}; + for (std::size_t i = 0U; i < n; i++) { + dCachedEdge e1{}; + if (x.isIdentity() || x.p->v < var) { + // [ 0 | 1 ] [ x | 0 ] + // --------- = --------- + // [ 2 | 3 ] [ 0 | x ] + if (i == 0 || i == 3) { + e1 = x; + } + } else { + auto& xSuccessor = x.p->e[i]; + e1 = {xSuccessor.p, 0}; + if (!xSuccessor.w.exactlyZero()) { + e1.w = x.w * xSuccessor.w; + } + } + dCachedEdge e2{}; + if (y.isIdentity() || y.p->v < var) { + // [ 0 | 1 ] [ y | 0 ] + // --------- = --------- + // [ 2 | 3 ] [ 0 | y ] + if (i == 0 || i == 3) { + e2 = y; + } + } else { + auto& ySuccessor = y.p->e[i]; + e2 = {ySuccessor.p, 0}; + if (!ySuccessor.w.exactlyZero()) { + e2.w = y.w * ySuccessor.w; + } + } + + dNode::applyDmChangesToNode(e1.p); + dNode::applyDmChangesToNode(e2.p); + edge[i] = add2(e1, e2, var - 1); + dNode::revertDmChangesToNode(e2.p); + dNode::revertDmChangesToNode(e1.p); + } + auto r = makeDDNode(var, edge); + densityAdd.insert(x, y, r); + return r; +} + +///----------------------------------------------------------------------------- +/// \n Multiplication \n +///----------------------------------------------------------------------------- + +dEdge DensityDDPackage::multiply(const dEdge& x, const dEdge& y, + const bool generateDensityMatrix) { + dd::Qubit var{}; + auto xCopy = x; + auto yCopy = y; + dEdge::applyDmChangesToEdges(xCopy, yCopy); + + if (!xCopy.isTerminal()) { + var = xCopy.p->v; + } + if (!y.isTerminal() && yCopy.p->v > var) { + var = yCopy.p->v; + } + + const auto e = multiply2(xCopy, yCopy, var, generateDensityMatrix); + dEdge::revertDmChangesToEdges(xCopy, yCopy); + return dEdge{e.p, pkg->cn.lookup(e.w)}; +} + +dCachedEdge DensityDDPackage::multiply2(const dEdge& x, const dEdge& y, + const dd::Qubit var, + const bool generateDensityMatrix) { + using ResultEdge = dCachedEdge; + + if (x.w.exactlyZero() || y.w.exactlyZero()) { + return ResultEdge::zero(); + } + + const auto xWeight = static_cast(x.w); + const auto yWeight = static_cast(y.w); + const auto rWeight = xWeight * yWeight; + if (x.isIdentity()) { + if (y.isIdentity() || + (dNode::isDensityMatrixTempFlagSet(y.p->flags) && + generateDensityMatrix) || + (!dNode::isDensityMatrixTempFlagSet(y.p->flags) && + !generateDensityMatrix)) { + return {y.p, rWeight}; + } + } + + if (y.isIdentity()) { + if (x.isIdentity() || + (dNode::isDensityMatrixTempFlagSet(x.p->flags) && + generateDensityMatrix) || + (!dNode::isDensityMatrixTempFlagSet(x.p->flags) && + !generateDensityMatrix)) { + return {x.p, rWeight}; + } + } + + if (const auto* r = + densityDensityMultiplication.lookup(x.p, y.p, generateDensityMatrix); + r != nullptr) { + return {r->p, r->w * rWeight}; + } + + constexpr std::size_t n = dd::NEDGE; + constexpr std::size_t rows = dd::RADIX; + constexpr std::size_t cols = dd::RADIX; + + std::array edge{}; + for (auto i = 0U; i < rows; i++) { + for (auto j = 0U; j < cols; j++) { + auto idx = (cols * i) + j; + edge[idx] = ResultEdge::zero(); + for (auto k = 0U; k < rows; k++) { + const auto xIdx = (rows * i) + k; + dEdge e1{}; + if (x.p != nullptr && x.p->v == var) { + e1 = x.p->e[xIdx]; + } else { + if (xIdx == 0 || xIdx == 3) { + e1 = dEdge{x.p, dd::Complex::one()}; + } else { + e1 = dEdge::zero(); + } + } + + const auto yIdx = j + (cols * k); + dEdge e2{}; + if (y.p != nullptr && y.p->v == var) { + e2 = y.p->e[yIdx]; + } else { + if (yIdx == 0 || yIdx == 3) { + e2 = dEdge{y.p, dd::Complex::one()}; + } else { + e2 = dEdge::zero(); + } + } + + const auto v = static_cast(var - 1); + dCachedEdge m; + dEdge::applyDmChangesToEdges(e1, e2); + if (!generateDensityMatrix || idx == 1) { + // When generateDensityMatrix is false or I have the first edge I + // don't optimize anything and set generateDensityMatrix to false + // for all child edges + m = multiply2(e1, e2, v, false); + } else if (idx == 2) { + // When I have the second edge and generateDensityMatrix == false, + // then edge[2] == edge[1] + if (k == 0) { + if (edge[1].w.approximatelyZero()) { + edge[2] = ResultEdge::zero(); + } else { + edge[2] = edge[1]; + } + } + continue; + } else { + m = multiply2(e1, e2, v, generateDensityMatrix); + } + + if (k == 0 || edge[idx].w.exactlyZero()) { + edge[idx] = m; + } else if (!m.w.exactlyZero()) { + dNode::applyDmChangesToNode(edge[idx].p); + dNode::applyDmChangesToNode(m.p); + edge[idx] = add2(edge[idx], m, v); + dNode::revertDmChangesToNode(m.p); + dNode::revertDmChangesToNode(edge[idx].p); + } + // Undo modifications on density matrices + dEdge::revertDmChangesToEdges(e1, e2); + } + } + } + + auto e = makeDDNode(var, edge, generateDensityMatrix); + densityDensityMultiplication.insert(x.p, y.p, e); + + e.w = e.w * rWeight; + return e; +} + +///----------------------------------------------------------------------------- +/// \n Trace \n +///----------------------------------------------------------------------------- + +dd::ComplexValue DensityDDPackage::trace(const dEdge& a, + const std::size_t numQubits) { + if (a.isIdentity()) { + return static_cast(a.w); + } + const auto eliminate = std::vector(numQubits, true); + return trace(a, eliminate, numQubits).w; +} + +dCachedEdge DensityDDPackage::trace(const dEdge& a, + const std::vector& eliminate, + std::size_t level, + std::size_t alreadyEliminated) { + const auto aWeight = static_cast(a.w); + if (aWeight.approximatelyZero()) { + return dCachedEdge::zero(); + } + + // If `a` is the identity matrix or there is nothing left to eliminate, + // then simply return `a` + if (a.isIdentity() || + std::none_of(eliminate.begin(), + eliminate.begin() + + static_cast::difference_type>(level), + [](bool v) { return v; })) { + return dCachedEdge{a.p, aWeight}; + } + + const auto v = a.p->v; + if (eliminate[v]) { + const auto eliminateAll = + std::all_of(eliminate.begin(), + eliminate.begin() + + static_cast::difference_type>(level), + [](bool e) { return e; }); + if (eliminateAll) { + if (const auto* r = densityTrace.lookup(a.p); r != nullptr) { + return {r->p, r->w * aWeight}; + } + } + + const auto elims = alreadyEliminated + 1; + auto r = add2(trace(a.p->e[0], eliminate, level - 1, elims), + trace(a.p->e[3], eliminate, level - 1, elims), v - 1); + + // Unlike for matrices, no normalization is applied for density matrices as + // their trace is always 1 by definition. + + if (eliminateAll) { + densityTrace.insert(a.p, r); + } + r.w = r.w * aWeight; + return r; + } + + std::array edge{}; + std::transform(a.p->e.cbegin(), a.p->e.cend(), edge.begin(), + [this, &eliminate, &alreadyEliminated, + &level](const dEdge& e) -> dCachedEdge { + return trace(e, eliminate, level - 1, alreadyEliminated); + }); + const auto adjustedV = + static_cast(static_cast(a.p->v) - + (static_cast(std::count( + eliminate.begin(), eliminate.end(), true)) - + alreadyEliminated)); + auto r = makeDDNode(adjustedV, edge); + r.w = r.w * aWeight; + return r; +} + +///----------------------------------------------------------------------------- +/// \n High-level operations \n +///----------------------------------------------------------------------------- + +dEdge DensityDDPackage::makeZeroDensityOperator(const std::size_t n) { + auto f = dEdge::one(); + for (std::size_t p = 0; p < n; p++) { + f = makeDDNode(static_cast(p), + std::array{f, dEdge::zero(), dEdge::zero(), dEdge::zero()}); + } + incRef(f); + return f; +} + +dEdge DensityDDPackage::applyOperationToDensity(dEdge& e, + const dd::mEdge& operation) { + const auto tmp0 = pkg->conjugateTranspose(operation); + const auto tmp1 = multiply(e, dd::ddsim::densityFromMatrixEdge(tmp0), false); + const auto tmp2 = + multiply(dd::ddsim::densityFromMatrixEdge(operation), tmp1, true); + incRef(tmp2); + dEdge::alignDensityEdge(e); + decRef(e); + e = tmp2; + dEdge::setDensityMatrixTrue(e); + return e; +} + +char DensityDDPackage::measureOneCollapsing(dEdge& e, const dd::Qubit index, + std::mt19937_64& mt) { + char measuredResult = '0'; + dEdge::alignDensityEdge(e); + const auto nrQubits = e.p->v + 1U; + dEdge::setDensityMatrixTrue(e); + + auto const measZeroDd = pkg->makeGateDD(dd::MEAS_ZERO_MAT, index); + + auto tmp0 = pkg->conjugateTranspose(measZeroDd); + auto tmp1 = multiply(e, dd::ddsim::densityFromMatrixEdge(tmp0), false); + auto tmp2 = + multiply(dd::ddsim::densityFromMatrixEdge(measZeroDd), tmp1, true); + auto densityMatrixTrace = trace(tmp2, nrQubits); + + std::uniform_real_distribution dist(0., 1.); + if (const auto threshold = dist(mt); threshold > densityMatrixTrace.r) { + auto const measOneDd = pkg->makeGateDD(dd::MEAS_ONE_MAT, index); + tmp0 = pkg->conjugateTranspose(measOneDd); + tmp1 = multiply(e, dd::ddsim::densityFromMatrixEdge(tmp0), false); + tmp2 = multiply(dd::ddsim::densityFromMatrixEdge(measOneDd), tmp1, true); + measuredResult = '1'; + densityMatrixTrace = trace(tmp2, nrQubits); + } + + dEdge::alignDensityEdge(e); + tmp2.w = pkg->cn.lookup(e.w / densityMatrixTrace); // Normalize density matrix + incRef(tmp2); + decRef(e); + e = tmp2; + dEdge::setDensityMatrixTrue(e); + + return measuredResult; +} + +///----------------------------------------------------------------------------- +/// \n Reference counting and GC \n +///----------------------------------------------------------------------------- + +void DensityDDPackage::incRef(const dEdge& e) { + if (dEdge::trackingRequired(e)) { + ++dRoots[e]; + } + // Keep shared matrix nodes and complex numbers alive during the borrowed + // package's own garbage collection by registering the (aligned) density DD + // as a matrix root. + pkg->incRef(matrixFromDensityEdge(e)); +} + +void DensityDDPackage::decRef(const dEdge& e) { + if (dEdge::trackingRequired(e)) { + if (auto it = dRoots.find(e); it != dRoots.end()) { + if (--it->second == 0U) { + dRoots.erase(it); + } + } + } + pkg->decRef(matrixFromDensityEdge(e)); +} + +bool DensityDDPackage::garbageCollect(const bool force) { + if (!force && !dUniqueTable.possiblyNeedsCollection()) { + return false; + } + for (const auto& edge : dRoots) { + edge.first.mark(); + } + const bool collected = dUniqueTable.garbageCollect(force) > 0; + for (const auto& edge : dRoots) { + edge.first.unmark(); + } + if (collected) { + densityAdd.clear(); + densityDensityMultiplication.clear(); + densityTrace.clear(); + } + return collected; +} + +std::size_t DensityDDPackage::computeActiveNodeCount() const { + for (const auto& edge : dRoots) { + edge.first.mark(); + } + const auto count = dUniqueTable.countMarkedEntries(); + for (const auto& edge : dRoots) { + edge.first.unmark(); + } + return count; +} + +} // namespace dd::ddsim diff --git a/src/DensityNode.cpp b/src/DensityNode.cpp new file mode 100644 index 000000000..6d774975b --- /dev/null +++ b/src/DensityNode.cpp @@ -0,0 +1,400 @@ +/* + * Copyright (c) 2023 - 2026 Chair for Design Automation, TUM + * Copyright (c) 2025 - 2026 Munich Quantum Software Company GmbH + * All rights reserved. + * + * SPDX-License-Identifier: MIT + * + * Licensed under the MIT License + */ + +#include "DensityNode.hpp" + +#include "dd/Complex.hpp" +#include "dd/ComplexNumbers.hpp" +#include "dd/ComplexValue.hpp" +#include "dd/DDDefinitions.hpp" +#include "dd/MemoryManager.hpp" +#include "dd/RealNumber.hpp" +#include "ir/Definitions.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace dd::ddsim { + +///----------------------------------------------------------------------------- +/// \n dNode methods \n +///----------------------------------------------------------------------------- + +void dNode::setDensityMatrixNodeFlag(const bool densityMatrix) noexcept { + if (densityMatrix) { + flags = (flags | static_cast(8U)); + } else { + flags = (flags & static_cast(~8U)); + } +} + +std::uint8_t dNode::alignDensityNodeNode(dNode*& p) noexcept { + const auto flags = static_cast(getDensityMatrixTempFlags(p)); + // Get an aligned node + alignDensityNode(p); + + if (dNode::isTerminal(p)) { + return 0U; + } + + if (isNonReduceTempFlagSet(flags) && !isConjugateTempFlagSet(flags)) { + // nothing more to do for first edge path (inherited by all child paths) + return flags; + } + + if (!isConjugateTempFlagSet(flags)) { + p->e[2].w = dd::ComplexNumbers::conj(p->e[2].w); + setConjugateTempFlagTrue(p->e[2].p); + // Mark the first edge + setNonReduceTempFlagTrue(p->e[1].p); + + for (auto& edge : p->e) { + setDensityMatTempFlagTrue(edge.p); + } + + } else { + std::swap(p->e[2], p->e[1]); + for (auto& edge : p->e) { + edge.w = dd::ComplexNumbers::conj(edge.w); + setConjugateTempFlagTrue(edge.p); + setDensityMatTempFlagTrue(edge.p); + } + } + return flags; +} + +void dNode::getAlignedNodeRevertModificationsOnSubEdges(dNode* p) noexcept { + // Get an aligned node and revert the modifications on the sub edges + alignDensityNode(p); + + for (auto& edge : p->e) { + // remove the set properties from the node pointers of edge.p->e + alignDensityNode(edge.p); + } + + if (isNonReduceTempFlagSet(p->flags) && !isConjugateTempFlagSet(p->flags)) { + // nothing more to do for a first edge path + return; + } + + if (!isConjugateTempFlagSet(p->flags)) { + p->e[2].w = dd::ComplexNumbers::conj(p->e[2].w); + return; + } + for (auto& edge : p->e) { + edge.w = dd::ComplexNumbers::conj(edge.w); + } + std::swap(p->e[2], p->e[1]); +} + +void dNode::applyDmChangesToNode(dNode*& p) noexcept { + if (isDensityMatrixTempFlagSet(p)) { + const auto tmp = alignDensityNodeNode(p); + if (p == nullptr) { + return; + } + assert(getDensityMatrixTempFlags(p->flags) == 0); + p->flags = p->flags | tmp; + } +} + +void dNode::revertDmChangesToNode(dNode*& p) noexcept { + if (!dNode::isTerminal(p) && isDensityMatrixTempFlagSet(p->flags)) { + getAlignedNodeRevertModificationsOnSubEdges(p); + p->unsetTempDensityMatrixFlags(); + } +} + +///----------------------------------------------------------------------------- +/// \n General purpose dEdge methods \n +///----------------------------------------------------------------------------- + +auto dEdge::size() const -> std::size_t { + static constexpr std::size_t NODECOUNT_BUCKETS = 200000U; + static std::unordered_set visited{NODECOUNT_BUCKETS}; + visited.max_load_factor(10); + visited.clear(); + return size(visited); +} + +auto dEdge::size(std::unordered_set& visited) const + -> std::size_t { + visited.emplace(p); + std::size_t sum = 1U; + if (!isTerminal()) { + for (const auto& e : p->e) { + if (!visited.contains(e.p)) { + sum += e.size(visited); + } + } + } + return sum; +} + +void dEdge::mark() const noexcept { + w.mark(); + if (isTerminal() || p->isMarked()) { + return; + } + p->mark(); + for (const dEdge& e : p->e) { + e.mark(); + } +} + +void dEdge::unmark() const noexcept { + w.unmark(); + if (isTerminal() || !p->isMarked()) { + return; + } + p->unmark(); + for (const dEdge& e : p->e) { + e.unmark(); + } +} + +///----------------------------------------------------------------------------- +/// \n Normalization (matrix variant) \n +///----------------------------------------------------------------------------- + +auto dEdge::normalize(dNode* p, const std::array& e, + dd::MemoryManager& mm, dd::ComplexNumbers& cn) -> dEdge { + assert(p != nullptr && "Node pointer passed to normalize is null."); + const auto zero = std::array{e[0].w.exactlyZero(), e[1].w.exactlyZero(), + e[2].w.exactlyZero(), e[3].w.exactlyZero()}; + + if (std::all_of(zero.begin(), zero.end(), [](auto b) { return b; })) { + mm.returnEntry(*p); + return dEdge::zero(); + } + + const auto weights = std::array{static_cast(e[0].w), + static_cast(e[1].w), + static_cast(e[2].w), + static_cast(e[3].w)}; + + std::optional argMax = std::nullopt; + dd::fp maxMag2 = 0.; + auto maxVal = dd::Complex::one(); + // determine max amplitude + for (auto i = 0U; i < dd::NEDGE; ++i) { + if (zero[i]) { + p->e[i] = dEdge::zero(); + continue; + } + const auto& w = weights[i]; + if (!argMax.has_value()) { + argMax = i; + maxMag2 = w.mag2(); + maxVal = e[i].w; + } else { + if (const auto mag2 = w.mag2(); mag2 - maxMag2 > dd::RealNumber::eps) { + argMax = i; + maxMag2 = mag2; + maxVal = e[i].w; + } + } + } + assert(argMax.has_value() && "argMax should have been set by now"); + + const auto argMaxValue = *argMax; + const auto argMaxWeight = weights[argMaxValue]; + for (auto i = 0U; i < dd::NEDGE; ++i) { + if (zero[i]) { + continue; + } + if (i == argMaxValue) { + p->e[i] = {e[i].p, dd::Complex::one()}; + continue; + } + p->e[i] = {e[i].p, cn.lookup(weights[i] / argMaxWeight)}; + if (p->e[i].w.exactlyZero()) { + p->e[i].p = dNode::getTerminal(); + } + } + return dEdge{p, maxVal}; +} + +///----------------------------------------------------------------------------- +/// \n Methods for density matrix DDs \n +///----------------------------------------------------------------------------- + +auto dEdge::getSparseProbabilityVector(const std::size_t numQubits, + const dd::fp threshold) const + -> dd::SparsePVec { + if (numQubits == 0U) { + return {{0, static_cast>(w).real()}}; + } + + auto e = *this; + dEdge::alignDensityEdge(e); + + auto probabilities = dd::SparsePVec{}; + e.traverseDiagonal( + 1, 0, + [&probabilities](const std::size_t i, const dd::fp& prob) { + probabilities[i] = prob; + }, + numQubits, threshold); + return probabilities; +} + +auto dEdge::getSparseProbabilityVectorStrKeys(const std::size_t numQubits, + const dd::fp threshold) const + -> dd::SparsePVecStrKeys { + if (numQubits == 0U) { + return {{"0", static_cast>(w).real()}}; + } + + auto e = *this; + dEdge::alignDensityEdge(e); + const auto nqubits = static_cast(e.p->v) + 1U; + + auto probabilities = dd::SparsePVecStrKeys{}; + e.traverseDiagonal( + 1, 0, + [&probabilities, &nqubits](const std::size_t i, const dd::fp& prob) { + probabilities[dd::intToBinaryString(i, nqubits)] = prob; + }, + numQubits, threshold); + return probabilities; +} + +void dEdge::traverseDiagonal(const dd::fp& prob, const std::size_t i, + dd::ProbabilityFunc f, const std::size_t level, + const dd::fp threshold) const { + // calculate new accumulated probability + const auto c = static_cast>(w); + const auto val = prob * c.real(); + + if (val < threshold) { + return; + } + + if (level == 0) { + assert(isTerminal()); + f(i, val); + return; + } + + const auto nextLevel = static_cast(level - 1U); + if (isTerminal() || p->v < nextLevel) { + traverseDiagonal(prob, i, f, nextLevel, threshold); + traverseDiagonal(prob, i | (1ULL << nextLevel), f, nextLevel, threshold); + return; + } + + if (auto& e = p->e[0]; !e.w.exactlyZero()) { + e.traverseDiagonal(val, i, f, nextLevel, threshold); + } + if (auto& e = p->e[3]; !e.w.exactlyZero()) { + e.traverseDiagonal(val, i | (1ULL << nextLevel), f, nextLevel, threshold); + } +} + +///----------------------------------------------------------------------------- +/// \n dCachedEdge normalization \n +///----------------------------------------------------------------------------- + +auto dCachedEdge::normalize(dNode* p, + const std::array& e, + dd::MemoryManager& mm, dd::ComplexNumbers& cn) + -> dCachedEdge { + assert(p != nullptr && "Node pointer passed to normalize is null."); + const auto zero = + std::array{e[0].w.approximatelyZero(), e[1].w.approximatelyZero(), + e[2].w.approximatelyZero(), e[3].w.approximatelyZero()}; + + if (std::all_of(zero.begin(), zero.end(), [](auto b) { return b; })) { + mm.returnEntry(*p); + return dCachedEdge::zero(); + } + + std::optional argMax = std::nullopt; + dd::fp maxMag2 = 0.; + dd::ComplexValue maxVal = 1.; + // determine max amplitude + for (auto i = 0U; i < dd::NEDGE; ++i) { + if (zero[i]) { + continue; + } + const auto& w = e[i].w; + if (!argMax.has_value()) { + argMax = i; + maxMag2 = w.mag2(); + maxVal = w; + } else { + if (const auto mag2 = w.mag2(); mag2 - maxMag2 > dd::RealNumber::eps) { + argMax = i; + maxMag2 = mag2; + maxVal = w; + } + } + } + assert(argMax.has_value() && "argMax should have been set by now"); + + const auto argMaxValue = *argMax; + for (auto i = 0U; i < dd::NEDGE; ++i) { + // The approximation below is really important for numerical stability. + // An exactly zero check will lead to numerical instabilities. + if (zero[i]) { + p->e[i] = dEdge::zero(); + continue; + } + if (i == argMaxValue) { + p->e[i] = {e[i].p, dd::Complex::one()}; + continue; + } + p->e[i] = {e[i].p, cn.lookup(e[i].w / maxVal)}; + if (p->e[i].w.exactlyZero()) { + p->e[i].p = dNode::getTerminal(); + } + } + return dCachedEdge{p, maxVal}; +} + +} // namespace dd::ddsim + +///----------------------------------------------------------------------------- +/// \n Hash related code \n +///----------------------------------------------------------------------------- + +namespace std { + +std::size_t +hash::operator()(const dd::ddsim::dEdge& e) const noexcept { + const auto h1 = dd::murmur64(reinterpret_cast(e.p)); + const auto h2 = std::hash{}(e.w); + auto h3 = qc::combineHash(h1, h2); + if (e.isTerminal()) { + return h3; + } + assert(dd::ddsim::dNode::isDensityMatrixTempFlagSet(e.p) == false); + const auto h4 = dd::ddsim::dNode::getDensityMatrixTempFlags(e.p->flags); + h3 = qc::combineHash(h3, h4); + return h3; +} + +std::size_t hash::operator()( + const dd::ddsim::dCachedEdge& e) const noexcept { + const auto h1 = dd::murmur64(reinterpret_cast(e.p)); + const auto h2 = std::hash{}(e.w); + return qc::combineHash(h1, h2); +} + +} // namespace std diff --git a/src/DensityUniqueTable.cpp b/src/DensityUniqueTable.cpp new file mode 100644 index 000000000..ed06d184c --- /dev/null +++ b/src/DensityUniqueTable.cpp @@ -0,0 +1,163 @@ +/* + * Copyright (c) 2023 - 2026 Chair for Design Automation, TUM + * Copyright (c) 2025 - 2026 Munich Quantum Software Company GmbH + * All rights reserved. + * + * SPDX-License-Identifier: MIT + * + * Licensed under the MIT License + */ + +#include "DensityUniqueTable.hpp" + +#include "dd/MemoryManager.hpp" +#include "dd/Node.hpp" +#include "dd/statistics/UniqueTableStatistics.hpp" + +#include +#include +#include +#include +#include + +namespace dd::ddsim { + +DensityUniqueTable::DensityUniqueTable(dd::MemoryManager& manager, + const UniqueTableConfig& config) + : cfg(config), gcLimit(config.initialGCLimit), memoryManager(&manager), + tables(config.nVars, Table(config.nBuckets)), stats(config.nVars) { + for (auto& stat : stats) { + stat.entrySize = sizeof(Bucket); + stat.numBuckets = cfg.nBuckets; + } +} + +void DensityUniqueTable::resize(const std::size_t nVars) { + cfg.nVars = nVars; + tables.resize(nVars, Table(cfg.nBuckets)); + stats.resize(nVars); + for (auto& stat : stats) { + stat.entrySize = sizeof(Bucket); + stat.numBuckets = cfg.nBuckets; + } +} + +bool DensityUniqueTable::possiblyNeedsCollection() const { + return getNumEntries() >= gcLimit; +} + +std::size_t DensityUniqueTable::garbageCollect(const bool force) { + const std::size_t numEntriesBefore = getNumEntries(); + if ((!force && numEntriesBefore < gcLimit) || numEntriesBefore == 0U) { + return 0U; + } + + std::size_t v = 0U; + for (auto& table : tables) { + auto& stat = stats[v]; + ++stat.gcRuns; + for (auto& bucket : table) { + dd::NodeBase* p = bucket; + dd::NodeBase* lastp = nullptr; + while (p != nullptr) { + if (!p->isMarked()) { + dd::NodeBase* next = p->next(); + if (lastp == nullptr) { + bucket = next; + } else { + lastp->setNext(next); + } + memoryManager->returnEntry(*p); + p = next; + --stat.numEntries; + } else { + lastp = p; + p = p->next(); + } + } + } + ++v; + } + + const auto numEntries = getNumEntries(); + if (numEntries > gcLimit / 10 * 9) { + gcLimit = numEntries + cfg.initialGCLimit; + } + return numEntriesBefore - numEntries; +} + +void DensityUniqueTable::clear() { + for (auto& table : tables) { + for (auto& bucket : table) { + bucket = nullptr; + } + } + gcLimit = cfg.initialGCLimit; + for (auto& stat : stats) { + stat.reset(); + } +} + +const dd::UniqueTableStatistics& +DensityUniqueTable::getStats(const std::size_t idx) const noexcept { + return stats.at(idx); +} + +nlohmann::basic_json<> +DensityUniqueTable::getStatsJson(const bool includeIndividualTables) const { + if (std::ranges::all_of(stats, [](const dd::UniqueTableStatistics& stat) { + return stat.peakNumEntries == 0U; + })) { + return "unused"; + } + + dd::UniqueTableStatistics totalStats; + for (const auto& stat : stats) { + totalStats.entrySize = std::max(totalStats.entrySize, stat.entrySize); + totalStats.numBuckets += stat.numBuckets; + totalStats.numEntries += stat.numEntries; + totalStats.peakNumEntries += stat.peakNumEntries; + totalStats.collisions += stat.collisions; + totalStats.hits += stat.hits; + totalStats.lookups += stat.lookups; + totalStats.inserts += stat.inserts; + totalStats.gcRuns = std::max(totalStats.gcRuns, stat.gcRuns); + } + + nlohmann::basic_json<> j; + j["total"] = totalStats.json(); + if (includeIndividualTables) { + std::size_t v = 0U; + for (const auto& stat : stats) { + j[std::to_string(v)] = stat.json(); + ++v; + } + } + return j; +} + +std::size_t DensityUniqueTable::getNumEntries() const noexcept { + return std::accumulate( + stats.begin(), stats.end(), std::size_t{0}, + [](const std::size_t& sum, const dd::UniqueTableStatistics& stat) { + return sum + stat.numEntries; + }); +} + +std::size_t DensityUniqueTable::countMarkedEntries() const noexcept { + std::size_t count = 0U; + for (const auto& table : tables) { + for (auto* bucket : table) { + auto* p = bucket; + while (p != nullptr) { + if (p->isMarked()) { + ++count; + } + p = p->next(); + } + } + } + return count; +} + +} // namespace dd::ddsim diff --git a/src/DeterministicNoiseSimulator.cpp b/src/DeterministicNoiseSimulator.cpp index 9a29a6853..85e183e73 100644 --- a/src/DeterministicNoiseSimulator.cpp +++ b/src/DeterministicNoiseSimulator.cpp @@ -31,30 +31,32 @@ using CN = dd::ComplexNumbers; void DeterministicNoiseSimulator::initializeSimulation( const std::size_t nQubits) { - rootEdge = dd->makeZeroDensityOperator(static_cast(nQubits)); + rootEdge = densityDD.makeZeroDensityOperator(static_cast(nQubits)); } void DeterministicNoiseSimulator::applyOperationToState( std::unique_ptr& op) { auto operation = dd::getDD(*op, *Simulator::dd); - dd->applyOperationToDensity(DeterministicNoiseSimulator::rootEdge, operation); + densityDD.applyOperationToDensity(DeterministicNoiseSimulator::rootEdge, + operation); deterministicNoiseFunctionality.applyNoiseEffects( DeterministicNoiseSimulator::rootEdge, op); + densityDD.garbageCollect(); } char DeterministicNoiseSimulator::measure(const dd::Qubit i) { - return Simulator::dd->measureOneCollapsing( - rootEdge, static_cast(i), Simulator::mt); + return densityDD.measureOneCollapsing(rootEdge, static_cast(i), + Simulator::mt); } void DeterministicNoiseSimulator::reset(qc::NonUnitaryOperation* nonUnitaryOp) { for (const auto& qubit : nonUnitaryOp->getTargets()) { - auto const result = - dd->measureOneCollapsing(rootEdge, static_cast(qubit), mt); + auto const result = densityDD.measureOneCollapsing( + rootEdge, static_cast(qubit), mt); if (result == '1') { const auto x = qc::StandardOperation(qubit, qc::X); const auto operation = dd::getDD(x, *dd); - rootEdge = dd->applyOperationToDensity(rootEdge, operation); + rootEdge = densityDD.applyOperationToDensity(rootEdge, operation); } } } diff --git a/src/NoiseFunctionality.cpp b/src/NoiseFunctionality.cpp new file mode 100644 index 000000000..23d4b84d2 --- /dev/null +++ b/src/NoiseFunctionality.cpp @@ -0,0 +1,485 @@ +/* + * Copyright (c) 2023 - 2026 Chair for Design Automation, TUM + * Copyright (c) 2025 - 2026 Munich Quantum Software Company GmbH + * All rights reserved. + * + * SPDX-License-Identifier: MIT + * + * Licensed under the MIT License + */ + +#include "NoiseFunctionality.hpp" + +#include "DensityDDPackage.hpp" +#include "DensityNode.hpp" +#include "dd/Complex.hpp" +#include "dd/ComplexNumbers.hpp" +#include "dd/ComplexValue.hpp" +#include "dd/DDDefinitions.hpp" +#include "dd/GateMatrixDefinitions.hpp" +#include "dd/Package.hpp" +#include "ir/Definitions.hpp" +#include "ir/operations/OpType.hpp" +#include "ir/operations/Operation.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace { + +std::vector +initializeNoiseEffects(const std::string& cNoiseEffects) { + std::vector noiseOperationVector{}; + noiseOperationVector.reserve(cNoiseEffects.size()); + for (const auto noise : cNoiseEffects) { + switch (noise) { + case 'A': + noiseOperationVector.emplace_back(dd::ddsim::AmplitudeDamping); + break; + case 'P': + noiseOperationVector.emplace_back(dd::ddsim::PhaseFlip); + break; + case 'D': + noiseOperationVector.emplace_back(dd::ddsim::Depolarization); + break; + case 'I': + noiseOperationVector.emplace_back(dd::ddsim::Identity); + break; + default: + throw std::runtime_error("Unknown noise operation '" + cNoiseEffects + + "'\n"); + } + } + return noiseOperationVector; +} +} // namespace + +namespace dd::ddsim { +StochasticNoiseFunctionality::StochasticNoiseFunctionality( + dd::Package& dd, const std::size_t nq, const double gateNoiseProbability, + const double amplitudeDampingProb, const double multiQubitGateFactor, + const std::string& cNoiseEffects) + : package(&dd), nQubits(nq), dist(0.0, 1.0L), + noiseProbability(gateNoiseProbability), + noiseProbabilityMulti(gateNoiseProbability * multiQubitGateFactor), + sqrtAmplitudeDampingProbability(std::sqrt(amplitudeDampingProb)), + oneMinusSqrtAmplitudeDampingProbability( + std::sqrt(1 - amplitudeDampingProb)), + sqrtAmplitudeDampingProbabilityMulti(std::sqrt(gateNoiseProbability) * + multiQubitGateFactor), + oneMinusSqrtAmplitudeDampingProbabilityMulti( + std::sqrt(1 - (multiQubitGateFactor * amplitudeDampingProb))), + ampDampingTrue({0, sqrtAmplitudeDampingProbability, 0, 0}), + ampDampingTrueMulti({0, sqrtAmplitudeDampingProbabilityMulti, 0, 0}), + ampDampingFalse({1, 0, 0, oneMinusSqrtAmplitudeDampingProbability}), + ampDampingFalseMulti( + {1, 0, 0, oneMinusSqrtAmplitudeDampingProbabilityMulti}), + noiseEffects(initializeNoiseEffects(cNoiseEffects)), + identityDD(dd::Package::makeIdent()), + stochasticNoiseOperationCache(nq, StochasticNoiseKindEnd) { + sanityCheckOfNoiseProbabilities(gateNoiseProbability, amplitudeDampingProb, + multiQubitGateFactor); + package->incRef(identityDD); +} + +double StochasticNoiseFunctionality::getNoiseProbability( + const bool multiQubitNoiseFlag) const { + return multiQubitNoiseFlag ? noiseProbabilityMulti : noiseProbability; +} + +StochasticNoiseFunctionality::StochasticNoiseKind +StochasticNoiseFunctionality::getAmplitudeDampingOperationType( + const bool multiQubitNoiseFlag, const bool amplitudeDampingFlag) { + if (amplitudeDampingFlag) { + return multiQubitNoiseFlag ? StochMultiATrue : StochATrue; + } + return multiQubitNoiseFlag ? StochMultiAFalse : StochAFalse; +} + +dd::GateMatrix StochasticNoiseFunctionality::getAmplitudeDampingOperationMatrix( + const bool multiQubitNoiseFlag, const bool amplitudeDampingFlag) const { + if (amplitudeDampingFlag) { + return multiQubitNoiseFlag ? ampDampingTrueMulti : ampDampingTrue; + } + return multiQubitNoiseFlag ? ampDampingFalseMulti : ampDampingFalse; +} + +void StochasticNoiseFunctionality::applyNoiseOperation( + const std::set& targets, dd::mEdge operation, dd::vEdge& state, + std::mt19937_64& generator) { + const bool multiQubitOperation = targets.size() > 1; + + for (const auto& target : targets) { + auto stackedOperation = generateNoiseOperation(operation, target, generator, + false, multiQubitOperation); + auto tmp = package->multiply(stackedOperation, state); + + if (dd::ComplexNumbers::mag2(tmp.w) < dist(generator)) { + // The probability of amplitude damping does not only depend on the + // noise probability, but also the quantum state. Due to the + // normalization constraint of decision diagrams the probability for + // applying amplitude damping stands in the root edge weight, of the dd + // after the noise has been applied + stackedOperation = generateNoiseOperation(operation, target, generator, + true, multiQubitOperation); + tmp = package->multiply(stackedOperation, state); + } + tmp.w = dd::Complex::one(); + + package->incRef(tmp); + package->decRef(state); + state = tmp; + + // I only need to apply the operations once + operation = identityDD; + } +} + +dd::mEdge StochasticNoiseFunctionality::stackOperation( + const dd::mEdge& operation, const qc::Qubit target, + const StochasticNoiseKind noiseOperation, const dd::GateMatrix& matrix) { + const auto kind = static_cast(noiseOperation); + if (const auto* op = stochasticNoiseOperationCache.lookup(kind, target); + op != nullptr) { + return package->multiply(*op, operation); + } + const auto gateDD = package->makeGateDD(matrix, target); + stochasticNoiseOperationCache.insert(kind, target, gateDD); + return package->multiply(gateDD, operation); +} + +dd::mEdge StochasticNoiseFunctionality::generateNoiseOperation( + dd::mEdge operation, const qc::Qubit target, std::mt19937_64& generator, + const bool amplitudeDamping, const bool multiQubitOperation) { + for (const auto& noiseType : noiseEffects) { + const auto effect = noiseType == AmplitudeDamping + ? getAmplitudeDampingOperationType( + multiQubitOperation, amplitudeDamping) + : returnNoiseOperation(noiseType, dist(generator), + multiQubitOperation); + switch (effect) { + case StochIdentity: { + continue; + } + case StochMultiATrue: + case StochATrue: { + const dd::GateMatrix amplitudeDampingMatrix = + getAmplitudeDampingOperationMatrix(multiQubitOperation, true); + operation = + stackOperation(operation, target, effect, amplitudeDampingMatrix); + break; + } + case StochMultiAFalse: + case StochAFalse: { + const dd::GateMatrix amplitudeDampingMatrix = + getAmplitudeDampingOperationMatrix(multiQubitOperation, false); + operation = + stackOperation(operation, target, effect, amplitudeDampingMatrix); + break; + } + case StochX: { + operation = stackOperation(operation, target, effect, + dd::opToSingleQubitGateMatrix(qc::X)); + break; + } + case StochY: { + operation = stackOperation(operation, target, effect, + dd::opToSingleQubitGateMatrix(qc::Y)); + break; + } + case StochZ: { + operation = stackOperation(operation, target, effect, + dd::opToSingleQubitGateMatrix(qc::Z)); + break; + } + default: { + throw std::runtime_error("Unknown noise operation '" + + std::to_string(effect) + "'\n"); + } + } + } + return operation; +} + +StochasticNoiseFunctionality::StochasticNoiseKind +StochasticNoiseFunctionality::returnNoiseOperation( + const NoiseOperations noiseOperation, const double prob, + const bool multiQubitNoiseFlag) const { + switch (noiseOperation) { + case Depolarization: { + if (prob >= (getNoiseProbability(multiQubitNoiseFlag) * 0.75)) { + // prob > prob apply StochIdentity, also 25 % of the time when + // depolarization is applied nothing happens + return StochIdentity; + } + if (prob < (getNoiseProbability(multiQubitNoiseFlag) * 0.25)) { + // if 0 < prob < 0.25 (25 % of the time when applying depolarization) + // apply StochX + return StochX; + } + if (prob < (getNoiseProbability(multiQubitNoiseFlag) * 0.5)) { + // if 0.25 < prob < 0.5 (25 % of the time when applying depolarization) + // apply StochY + return StochY; + } + // if 0.5 < prob < 0.75 (25 % of the time when applying depolarization) + // apply StochZ + return StochZ; + } + case PhaseFlip: { + if (prob > getNoiseProbability(multiQubitNoiseFlag)) { + return StochIdentity; + } + return StochZ; + } + case Identity: { + return StochIdentity; + } + default: + throw std::runtime_error(std::string{"Unknown noise effect '"} + + std::to_string(noiseOperation) + "'"); + } +} + +DeterministicNoiseFunctionality::DeterministicNoiseFunctionality( + DensityDDPackage& dd, const std::size_t nq, + const double noiseProbabilitySingleQubit, + const double noiseProbabilityMultiQubit, + const double ampDampProbSingleQubit, const double ampDampProbMultiQubit, + const std::string& cNoiseEffects) + : package(&dd), nQubits(nq), + noiseProbSingleQubit(noiseProbabilitySingleQubit), + noiseProbMultiQubit(noiseProbabilityMultiQubit), + ampDampingProbSingleQubit(ampDampProbSingleQubit), + ampDampingProbMultiQubit(ampDampProbMultiQubit), + noiseEffects(initializeNoiseEffects(cNoiseEffects)) { + sanityCheckOfNoiseProbabilities(noiseProbabilitySingleQubit, + ampDampProbSingleQubit, 1); + sanityCheckOfNoiseProbabilities(noiseProbabilityMultiQubit, + ampDampProbMultiQubit, 1); +} + +void DeterministicNoiseFunctionality::applyNoiseEffects( + dEdge& originalEdge, const std::unique_ptr& qcOperation) { + const auto usedQubits = qcOperation->getUsedQubits(); + dCachedEdge nodeAfterNoise = {}; + dEdge::applyDmChangesToEdge(originalEdge); + nodeAfterNoise = applyNoiseEffects(originalEdge, usedQubits, false, + static_cast(nQubits)); + dEdge::revertDmChangesToEdge(originalEdge); + const auto r = + dEdge{nodeAfterNoise.p, package->package().cn.lookup(nodeAfterNoise.w)}; + package->incRef(r); + dEdge::alignDensityEdge(originalEdge); + package->decRef(originalEdge); + originalEdge = r; + dEdge::setDensityMatrixTrue(originalEdge); +} + +dCachedEdge DeterministicNoiseFunctionality::applyNoiseEffects( + dEdge& originalEdge, const std::set& usedQubits, + const bool firstPathEdge, const dd::Qubit level) { + + const auto originalWeight = static_cast(originalEdge.w); + if (originalEdge.isZeroTerminal() || level <= *usedQubits.begin()) { + return {originalEdge.p, originalWeight}; + } + + auto originalCopy = dEdge{originalEdge.p, dd::Complex::one()}; + ArrayOfEdges newEdges{}; + const auto nextLevel = static_cast(level - 1U); + if (originalEdge.isIdentity()) { + newEdges[0] = + applyNoiseEffects(originalCopy, usedQubits, firstPathEdge, nextLevel); + newEdges[3] = + applyNoiseEffects(originalCopy, usedQubits, firstPathEdge, nextLevel); + } else { + for (std::size_t i = 0; i < newEdges.size(); i++) { + auto& successor = originalCopy.p->e[i]; + if (firstPathEdge || i == 1) { + // If I am to the firstPathEdge I cannot minimize the necessary + // operations anymore + dEdge::applyDmChangesToEdge(successor); + newEdges[i] = applyNoiseEffects(successor, usedQubits, true, nextLevel); + dEdge::revertDmChangesToEdge(successor); + } else if (i == 2) { + // Since e[1] == e[2] (due to density matrix representation), I can skip + // calculating e[2] + newEdges[2] = newEdges[1]; + } else { + dEdge::applyDmChangesToEdge(successor); + newEdges[i] = + applyNoiseEffects(successor, usedQubits, false, nextLevel); + dEdge::revertDmChangesToEdge(successor); + } + } + } + if (std::ranges::any_of(usedQubits, [&nextLevel](const qc::Qubit qubit) { + return nextLevel == qubit; + })) { + for (auto const& type : noiseEffects) { + switch (type) { + case AmplitudeDamping: + applyAmplitudeDampingToEdges(newEdges, (usedQubits.size() == 1) + ? ampDampingProbSingleQubit + : ampDampingProbMultiQubit); + break; + case PhaseFlip: + applyPhaseFlipToEdges(newEdges, (usedQubits.size() == 1) + ? noiseProbSingleQubit + : noiseProbMultiQubit); + break; + case Depolarization: + applyDepolarisationToEdges(newEdges, (usedQubits.size() == 1) + ? noiseProbSingleQubit + : noiseProbMultiQubit); + break; + case Identity: + continue; + } + } + } + + auto e = package->makeDDNode(nextLevel, newEdges, firstPathEdge); + if (e.w.exactlyZero()) { + return e; + } + e.w = e.w * originalWeight; + return e; +} + +void DeterministicNoiseFunctionality::applyPhaseFlipToEdges( + ArrayOfEdges& e, const double probability) { + const auto complexProb = 1. - (2. * probability); + + // e[0] = e[0] + // e[1] = (1-2p)*e[1] + if (!e[1].w.exactlyZero()) { + e[1].w *= complexProb; + } + // e[2] = (1-2p)*e[2] + if (!e[2].w.exactlyZero()) { + e[2].w *= complexProb; + } + // e[3] = e[3] +} + +void DeterministicNoiseFunctionality::applyAmplitudeDampingToEdges( + ArrayOfEdges& e, const double probability) const { + // e[0] = e[0] + p*e[3] + if (!e[3].w.exactlyZero()) { + if (!e[0].w.exactlyZero()) { + const auto var = static_cast(std::max( + {e[0].p != nullptr ? e[0].p->v : 0, e[1].p != nullptr ? e[1].p->v : 0, + e[2].p != nullptr ? e[2].p->v : 0, + e[3].p != nullptr ? e[3].p->v : 0})); + e[0] = package->add2(e[0], {e[3].p, e[3].w * probability}, var); + } else { + e[0] = {e[3].p, e[3].w * probability}; + } + } + + // e[1] = sqrt(1-p)*e[1] + if (!e[1].w.exactlyZero()) { + e[1].w *= std::sqrt(1 - probability); + } + + // e[2] = sqrt(1-p)*e[2] + if (!e[2].w.exactlyZero()) { + e[2].w *= std::sqrt(1 - probability); + } + + // e[3] = (1-p)*e[3] + if (!e[3].w.exactlyZero()) { + e[3].w *= (1 - probability); + } +} + +void DeterministicNoiseFunctionality::applyDepolarisationToEdges( + ArrayOfEdges& e, const double probability) const { + std::array helperEdge{}; + + const auto var = static_cast(std::max( + {e[0].p != nullptr ? e[0].p->v : 0, e[1].p != nullptr ? e[1].p->v : 0, + e[2].p != nullptr ? e[2].p->v : 0, e[3].p != nullptr ? e[3].p->v : 0})); + + const auto oldE0Edge = e[0]; + + // e[0] = 0.5*((2-p)*e[0] + p*e[3]) + { + // helperEdge[0] = 0.5*((2-p)*e[0] + helperEdge[0].p = e[0].p; + if (!e[0].w.exactlyZero()) { + helperEdge[0].w = e[0].w * (2 - probability) * 0.5; + } else { + helperEdge[0].w = 0; + } + + // helperEdge[1] = 0.5*p*e[3] + helperEdge[1].p = e[3].p; + if (!e[3].w.exactlyZero()) { + helperEdge[1].w = e[3].w * probability * 0.5; + } else { + helperEdge[1].w = 0; + } + + // e[0] = helperEdge[0] + helperEdge[1] + e[0] = package->add2(helperEdge[0], helperEdge[1], var); + } + + // e[1]=(1-p)*e[1] + if (!e[1].w.exactlyZero()) { + e[1].w *= (1 - probability); + } + // e[2]=(1-p)*e[2] + if (!e[2].w.exactlyZero()) { + e[2].w *= (1 - probability); + } + + // e[3] = 0.5*((2-p)*e[3]) + 0.5*(p*e[0]) + { + // helperEdge[0] = 0.5*((2-p)*e[3]) + helperEdge[0].p = e[3].p; + if (!e[3].w.exactlyZero()) { + helperEdge[0].w = e[3].w * (2 - probability) * 0.5; + } else { + helperEdge[0].w = 0; + } + + // helperEdge[1] = 0.5*p*e[0] + helperEdge[1].p = oldE0Edge.p; + if (!oldE0Edge.w.exactlyZero()) { + helperEdge[1].w = oldE0Edge.w * probability * 0.5; + } else { + helperEdge[1].w = 0; + } + e[3] = package->add2(helperEdge[0], helperEdge[1], var); + } +} + +void sanityCheckOfNoiseProbabilities(const double noiseProbability, + const double amplitudeDampingProb, + const double multiQubitGateFactor) { + if (noiseProbability < 0 || amplitudeDampingProb < 0 || + noiseProbability * multiQubitGateFactor > 1 || + amplitudeDampingProb * multiQubitGateFactor > 1) { + throw std::runtime_error( + "Error probabilities are faulty!" + "\n single qubit error probability: " + + std::to_string(noiseProbability) + " multi qubit error probability: " + + std::to_string(noiseProbability * multiQubitGateFactor) + + "\n single qubit amplitude damping probability: " + + std::to_string(amplitudeDampingProb) + + " multi qubit amplitude damping probability: " + + std::to_string(amplitudeDampingProb * multiQubitGateFactor)); + } +} +} // namespace dd::ddsim diff --git a/src/StochasticNoiseSimulator.cpp b/src/StochasticNoiseSimulator.cpp index debb748e8..052537ba2 100644 --- a/src/StochasticNoiseSimulator.cpp +++ b/src/StochasticNoiseSimulator.cpp @@ -11,7 +11,6 @@ #include "StochasticNoiseSimulator.hpp" #include "dd/DDDefinitions.hpp" -#include "dd/DDpackageConfig.hpp" #include "dd/Node.hpp" #include "dd/Operations.hpp" #include "dd/Package.hpp" @@ -78,10 +77,13 @@ void StochasticNoiseSimulator::runStochSimulationForId( for (std::size_t currentRun = 0U; currentRun < numberOfRuns; currentRun++) { auto localDD = std::make_unique( - getNumberOfQubits(), dd::STOCHASTIC_NOISE_SIMULATOR_DD_PACKAGE_CONFIG); - auto stochasticNoiseFunctionality = dd::StochasticNoiseFunctionality( + getNumberOfQubits(), + dd::ddsim::STOCHASTIC_NOISE_SIMULATOR_DD_PACKAGE_CONFIG); + auto stochasticNoiseFunctionality = dd::ddsim::StochasticNoiseFunctionality( *localDD, static_cast(nQubits), noiseProbability, amplitudeDampingProb, multiQubitGateFactor, noiseEffects); + dd::ddsim::StochasticNoiseOperationTable gateCache( + getNumberOfQubits(), qc::OpType::OpTypeEnd); std::vector classicValues(qc->getNcbits(), false); @@ -118,13 +120,12 @@ void StochasticNoiseSimulator::runStochSimulationForId( const auto& controls = op->getControls(); if (targets.size() == 1 && controls.empty()) { - const auto* oper = localDD->stochasticNoiseOperationCache.lookup( + const auto* oper = gateCache.lookup( op->getType(), static_cast(targets.front())); if (oper == nullptr) { operation = getDD(*op, *localDD); - localDD->stochasticNoiseOperationCache.insert( - op->getType(), static_cast(targets.front()), - operation); + gateCache.insert(op->getType(), + static_cast(targets.front()), operation); } else { operation = *oper; } From a0fcf88363f952533290386ea01e6e9164b491c6 Mon Sep 17 00:00:00 2001 From: Daniel Haag <121057143+denialhaag@users.noreply.github.com> Date: Sun, 2 Aug 2026 16:26:12 +0200 Subject: [PATCH 2/5] Fix Windows build and linter errors Assisted-by: Claude Opus 4.8 via Claude Code --- include/DensityNode.hpp | 8 ++++---- src/DensityDDPackage.cpp | 21 +++++++++++---------- src/DensityNode.cpp | 18 +++++++++--------- 3 files changed, 24 insertions(+), 23 deletions(-) diff --git a/include/DensityNode.hpp b/include/DensityNode.hpp index 70d7a8f00..850e82aa4 100644 --- a/include/DensityNode.hpp +++ b/include/DensityNode.hpp @@ -56,7 +56,7 @@ struct dEdge { // NOLINT(readability-identifier-naming) static constexpr dEdge one() { return terminal(dd::Complex::one()); } [[nodiscard]] static constexpr dEdge terminal(const dd::Complex& w); - [[nodiscard]] static constexpr bool trackingRequired(const dEdge& e) { + [[nodiscard]] static bool trackingRequired(const dEdge& e) { return !e.isTerminal() || !dd::constants::isStaticNumber(e.w.r) || !dd::constants::isStaticNumber(e.w.i); } @@ -123,7 +123,7 @@ struct dEdge { // NOLINT(readability-identifier-naming) size(std::unordered_set& visited) const; void traverseDiagonal(const dd::fp& prob, std::size_t i, - dd::ProbabilityFunc f, std::size_t level, + const dd::ProbabilityFunc& f, std::size_t level, dd::fp threshold = 0.) const; }; @@ -270,7 +270,7 @@ struct dNode final : dd::NodeBase { // NOLINT(readability-identifier-naming) /// Reinterpret a Core matrix edge as a density-matrix edge. inline dEdge densityFromMatrixEdge(const dd::mEdge& e) { - return dEdge{reinterpret_cast(e.p), e.w}; + return dEdge{.p = reinterpret_cast(e.p), .w = e.w}; } ///----------------------------------------------------------------------------- @@ -278,7 +278,7 @@ inline dEdge densityFromMatrixEdge(const dd::mEdge& e) { ///----------------------------------------------------------------------------- constexpr dEdge dEdge::terminal(const dd::Complex& w) { - return dEdge{dNode::getTerminal(), w}; + return dEdge{.p = dNode::getTerminal(), .w = w}; } inline bool dEdge::isTerminal() const { return dNode::isTerminal(p); } diff --git a/src/DensityDDPackage.cpp b/src/DensityDDPackage.cpp index b4c9bc21a..9b443e545 100644 --- a/src/DensityDDPackage.cpp +++ b/src/DensityDDPackage.cpp @@ -49,12 +49,12 @@ dEdge DensityDDPackage::makeDDNode(const dd::Qubit var, es[2].w.exactlyZero() && es[3].w.exactlyOne())) { auto* ptr = es[0].p; mm.returnEntry(*e.p); - return dEdge{ptr, e.w}; + return dEdge{.p = ptr, .w = e.w}; } } auto* l = dUniqueTable.lookup(e.p); - return dEdge{l, e.w}; + return dEdge{.p = l, .w = e.w}; } dCachedEdge @@ -172,7 +172,7 @@ dEdge DensityDDPackage::multiply(const dEdge& x, const dEdge& y, const auto e = multiply2(xCopy, yCopy, var, generateDensityMatrix); dEdge::revertDmChangesToEdges(xCopy, yCopy); - return dEdge{e.p, pkg->cn.lookup(e.w)}; + return dEdge{.p = e.p, .w = pkg->cn.lookup(e.w)}; } dCachedEdge DensityDDPackage::multiply2(const dEdge& x, const dEdge& y, @@ -229,7 +229,7 @@ dCachedEdge DensityDDPackage::multiply2(const dEdge& x, const dEdge& y, e1 = x.p->e[xIdx]; } else { if (xIdx == 0 || xIdx == 3) { - e1 = dEdge{x.p, dd::Complex::one()}; + e1 = dEdge{.p = x.p, .w = dd::Complex::one()}; } else { e1 = dEdge::zero(); } @@ -241,7 +241,7 @@ dCachedEdge DensityDDPackage::multiply2(const dEdge& x, const dEdge& y, e2 = y.p->e[yIdx]; } else { if (yIdx == 0 || yIdx == 3) { - e2 = dEdge{y.p, dd::Complex::one()}; + e2 = dEdge{.p = y.p, .w = dd::Complex::one()}; } else { e2 = dEdge::zero(); } @@ -352,11 +352,12 @@ dCachedEdge DensityDDPackage::trace(const dEdge& a, } std::array edge{}; - std::transform(a.p->e.cbegin(), a.p->e.cend(), edge.begin(), - [this, &eliminate, &alreadyEliminated, - &level](const dEdge& e) -> dCachedEdge { - return trace(e, eliminate, level - 1, alreadyEliminated); - }); + std::ranges::transform(a.p->e, edge.begin(), + [this, &eliminate, &alreadyEliminated, + &level](const dEdge& e) -> dCachedEdge { + return trace(e, eliminate, level - 1, + alreadyEliminated); + }); const auto adjustedV = static_cast(static_cast(a.p->v) - (static_cast(std::count( diff --git a/src/DensityNode.cpp b/src/DensityNode.cpp index 6d774975b..5e05f006a 100644 --- a/src/DensityNode.cpp +++ b/src/DensityNode.cpp @@ -25,7 +25,6 @@ #include #include #include -#include #include #include @@ -178,7 +177,7 @@ auto dEdge::normalize(dNode* p, const std::array& e, const auto zero = std::array{e[0].w.exactlyZero(), e[1].w.exactlyZero(), e[2].w.exactlyZero(), e[3].w.exactlyZero()}; - if (std::all_of(zero.begin(), zero.end(), [](auto b) { return b; })) { + if (std::ranges::all_of(zero, [](auto b) { return b; })) { mm.returnEntry(*p); return dEdge::zero(); } @@ -219,15 +218,15 @@ auto dEdge::normalize(dNode* p, const std::array& e, continue; } if (i == argMaxValue) { - p->e[i] = {e[i].p, dd::Complex::one()}; + p->e[i] = {.p = e[i].p, .w = dd::Complex::one()}; continue; } - p->e[i] = {e[i].p, cn.lookup(weights[i] / argMaxWeight)}; + p->e[i] = {.p = e[i].p, .w = cn.lookup(weights[i] / argMaxWeight)}; if (p->e[i].w.exactlyZero()) { p->e[i].p = dNode::getTerminal(); } } - return dEdge{p, maxVal}; + return dEdge{.p = p, .w = maxVal}; } ///----------------------------------------------------------------------------- @@ -276,7 +275,8 @@ auto dEdge::getSparseProbabilityVectorStrKeys(const std::size_t numQubits, } void dEdge::traverseDiagonal(const dd::fp& prob, const std::size_t i, - dd::ProbabilityFunc f, const std::size_t level, + const dd::ProbabilityFunc& f, + const std::size_t level, const dd::fp threshold) const { // calculate new accumulated probability const auto c = static_cast>(w); @@ -320,7 +320,7 @@ auto dCachedEdge::normalize(dNode* p, std::array{e[0].w.approximatelyZero(), e[1].w.approximatelyZero(), e[2].w.approximatelyZero(), e[3].w.approximatelyZero()}; - if (std::all_of(zero.begin(), zero.end(), [](auto b) { return b; })) { + if (std::ranges::all_of(zero, [](auto b) { return b; })) { mm.returnEntry(*p); return dCachedEdge::zero(); } @@ -357,10 +357,10 @@ auto dCachedEdge::normalize(dNode* p, continue; } if (i == argMaxValue) { - p->e[i] = {e[i].p, dd::Complex::one()}; + p->e[i] = {.p = e[i].p, .w = dd::Complex::one()}; continue; } - p->e[i] = {e[i].p, cn.lookup(e[i].w / maxVal)}; + p->e[i] = {.p = e[i].p, .w = cn.lookup(e[i].w / maxVal)}; if (p->e[i].w.exactlyZero()) { p->e[i].p = dNode::getTerminal(); } From d26abc7b5fb1b956855b581f77d2480d2b237b3c Mon Sep 17 00:00:00 2001 From: Daniel Haag <121057143+denialhaag@users.noreply.github.com> Date: Sun, 2 Aug 2026 17:07:08 +0200 Subject: [PATCH 3/5] Fix linter errors Assisted-by: Claude Opus 4.8 via Claude Code --- src/DensityNode.cpp | 1 + src/NoiseFunctionality.cpp | 8 +++++--- src/StochasticNoiseSimulator.cpp | 2 ++ 3 files changed, 8 insertions(+), 3 deletions(-) diff --git a/src/DensityNode.cpp b/src/DensityNode.cpp index 5e05f006a..91706b0ae 100644 --- a/src/DensityNode.cpp +++ b/src/DensityNode.cpp @@ -24,6 +24,7 @@ #include #include #include +#include #include #include #include diff --git a/src/NoiseFunctionality.cpp b/src/NoiseFunctionality.cpp index 23d4b84d2..b8a7107c1 100644 --- a/src/NoiseFunctionality.cpp +++ b/src/NoiseFunctionality.cpp @@ -17,6 +17,7 @@ #include "dd/ComplexValue.hpp" #include "dd/DDDefinitions.hpp" #include "dd/GateMatrixDefinitions.hpp" +#include "dd/Node.hpp" #include "dd/Package.hpp" #include "ir/Definitions.hpp" #include "ir/operations/OpType.hpp" @@ -26,6 +27,7 @@ #include #include #include +#include #include #include #include @@ -275,8 +277,8 @@ void DeterministicNoiseFunctionality::applyNoiseEffects( nodeAfterNoise = applyNoiseEffects(originalEdge, usedQubits, false, static_cast(nQubits)); dEdge::revertDmChangesToEdge(originalEdge); - const auto r = - dEdge{nodeAfterNoise.p, package->package().cn.lookup(nodeAfterNoise.w)}; + const auto r = dEdge{.p = nodeAfterNoise.p, + .w = package->package().cn.lookup(nodeAfterNoise.w)}; package->incRef(r); dEdge::alignDensityEdge(originalEdge); package->decRef(originalEdge); @@ -293,7 +295,7 @@ dCachedEdge DeterministicNoiseFunctionality::applyNoiseEffects( return {originalEdge.p, originalWeight}; } - auto originalCopy = dEdge{originalEdge.p, dd::Complex::one()}; + auto originalCopy = dEdge{.p = originalEdge.p, .w = dd::Complex::one()}; ArrayOfEdges newEdges{}; const auto nextLevel = static_cast(level - 1U); if (originalEdge.isIdentity()) { diff --git a/src/StochasticNoiseSimulator.cpp b/src/StochasticNoiseSimulator.cpp index 052537ba2..fdb39cc52 100644 --- a/src/StochasticNoiseSimulator.cpp +++ b/src/StochasticNoiseSimulator.cpp @@ -10,6 +10,8 @@ #include "StochasticNoiseSimulator.hpp" +#include "DensityDDPackage.hpp" +#include "StochasticNoiseOperationTable.hpp" #include "dd/DDDefinitions.hpp" #include "dd/Node.hpp" #include "dd/Operations.hpp" From fb98e5c1119808720c009c1bfed4d4bc6a6fc695 Mon Sep 17 00:00:00 2001 From: Daniel Haag <121057143+denialhaag@users.noreply.github.com> Date: Tue, 4 Aug 2026 21:23:29 +0200 Subject: [PATCH 4/5] Improve test coverage Assisted-by: Claude Opus 4.8 via Claude Code --- test/CMakeLists.txt | 1 + test/test_noise_functionality.cpp | 322 ++++++++++++++++++++++++++++++ 2 files changed, 323 insertions(+) create mode 100644 test/test_noise_functionality.cpp diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 5ab6e4cd1..9d36628a6 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -16,6 +16,7 @@ package_add_test( test_hybridsim.cpp test_stoch_noise_sim.cpp test_det_noise_sim.cpp + test_noise_functionality.cpp test_unitary_sim.cpp test_path_sim.cpp) diff --git a/test/test_noise_functionality.cpp b/test/test_noise_functionality.cpp new file mode 100644 index 000000000..18311638e --- /dev/null +++ b/test/test_noise_functionality.cpp @@ -0,0 +1,322 @@ +/* + * Copyright (c) 2023 - 2026 Chair for Design Automation, TUM + * Copyright (c) 2025 - 2026 Munich Quantum Software Company GmbH + * All rights reserved. + * + * SPDX-License-Identifier: MIT + * + * Licensed under the MIT License + */ + +#include "DensityDDPackage.hpp" +#include "DensityNode.hpp" +#include "NoiseFunctionality.hpp" +#include "dd/DDDefinitions.hpp" +#include "dd/Operations.hpp" +#include "dd/Package.hpp" +#include "dd/StateGeneration.hpp" +#include "ir/QuantumComputation.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +class DDNoiseFunctionalityTest : public ::testing::Test { +protected: + void SetUp() override { + // circuit taken from https://github.com/pnnl/qasmbench + qc.addQubitRegister(4U); + qc.x(0); + qc.x(1); + qc.h(3); + qc.cx(2, 3); + qc.t(0); + qc.t(1); + qc.t(2); + qc.tdg(3); + qc.cx(0, 1); + qc.cx(2, 3); + qc.cx(3, 0); + qc.cx(1, 2); + qc.cx(0, 1); + qc.cx(2, 3); + qc.tdg(0); + qc.tdg(1); + qc.tdg(2); + qc.t(3); + qc.cx(0, 1); + qc.cx(2, 3); + qc.s(3); + qc.cx(3, 0); + qc.h(3); + } + + qc::QuantumComputation qc; + std::size_t stochRuns = 1000U; +}; + +TEST_F(DDNoiseFunctionalityTest, DetSimulateAdder4TrackAPD) { + const dd::SparsePVecStrKeys reference = { + {"0000", 0.0969332192741}, {"1000", 0.0907888041538}, + {"0100", 0.0141409660985}, {"1100", 0.0092413539333}, + {"0010", 0.0238203475524}, {"1010", 0.0235097990017}, + {"0110", 0.0244576087400}, {"1110", 0.0116282811276}, + {"0001", 0.1731941264570}, {"1001", 0.4145855071998}, + {"0101", 0.0138062113213}, {"1101", 0.0184033482066}, + {"0011", 0.0242454336917}, {"1011", 0.0262779844799}, + {"0111", 0.0239296920989}, {"1111", 0.0110373166627}}; + + auto dd = std::make_unique( + qc.getNqubits(), dd::ddsim::DENSITY_MATRIX_SIMULATOR_DD_PACKAGE_CONFIG); + dd::ddsim::DensityDDPackage densityDD(*dd, qc.getNqubits()); + + auto rootEdge = densityDD.makeZeroDensityOperator(qc.getNqubits()); + + const auto* const noiseEffects = "APDI"; + + auto deterministicNoiseFunctionality = + dd::ddsim::DeterministicNoiseFunctionality( + densityDD, qc.getNqubits(), 0.01, 0.02, 0.02, 0.04, noiseEffects); + + for (auto const& op : qc) { + densityDD.applyOperationToDensity(rootEdge, dd::getDD(*op, *dd)); + deterministicNoiseFunctionality.applyNoiseEffects(rootEdge, op); + } + + // Expect that all results are the same + const auto m = + rootEdge.getSparseProbabilityVectorStrKeys(qc.getNqubits(), 0.001); + static constexpr dd::fp TOLERANCE = 1e-10; + for (const auto& [key, value] : m) { + EXPECT_NEAR(value, reference.at(key), TOLERANCE); + } +} + +TEST_F(DDNoiseFunctionalityTest, DetSimulateAdder4TrackD) { + const dd::SparsePVecStrKeys reference = { + {"0000", 0.0332328704931}, {"0001", 0.0683938280189}, + {"0011", 0.0117061689898}, {"0100", 0.0129643065735}, + {"0101", 0.0107812802908}, {"0111", 0.0160082331009}, + {"1000", 0.0328434857577}, {"1001", 0.7370101351171}, + {"1011", 0.0186346925411}, {"1101", 0.0275086747656}}; + + auto dd = std::make_unique( + qc.getNqubits(), dd::ddsim::DENSITY_MATRIX_SIMULATOR_DD_PACKAGE_CONFIG); + dd::ddsim::DensityDDPackage densityDD(*dd, qc.getNqubits()); + + auto rootEdge = densityDD.makeZeroDensityOperator(qc.getNqubits()); + + const auto* const noiseEffects = "D"; + + auto deterministicNoiseFunctionality = + dd::ddsim::DeterministicNoiseFunctionality( + densityDD, qc.getNqubits(), 0.01, 0.02, 0.02, 0.04, noiseEffects); + + for (auto const& op : qc) { + densityDD.applyOperationToDensity(rootEdge, dd::getDD(*op, *dd)); + deterministicNoiseFunctionality.applyNoiseEffects(rootEdge, op); + } + + // Expect that all results are the same + const auto m = + rootEdge.getSparseProbabilityVectorStrKeys(qc.getNqubits(), 0.01); + static constexpr dd::fp TOLERANCE = 1e-10; + for (const auto& [key, value] : m) { + EXPECT_NEAR(value, reference.at(key), TOLERANCE); + } +} + +TEST_F(DDNoiseFunctionalityTest, testingMeasure) { + constexpr double tolerance = 1e-10; + + qc::QuantumComputation qcOp{}; + qcOp.addQubitRegister(3U); + qcOp.h(0); + qcOp.h(1); + qcOp.h(2); + + auto dd = std::make_unique( + qcOp.getNqubits(), dd::ddsim::DENSITY_MATRIX_SIMULATOR_DD_PACKAGE_CONFIG); + dd::ddsim::DensityDDPackage densityDD(*dd, qcOp.getNqubits()); + + auto rootEdge = densityDD.makeZeroDensityOperator(qcOp.getNqubits()); + + auto deterministicNoiseFunctionality = + dd::ddsim::DeterministicNoiseFunctionality(densityDD, qcOp.getNqubits(), + 0.01, 0.02, 0.02, 0.04, {}); + + for (auto const& op : qcOp) { + densityDD.applyOperationToDensity(rootEdge, dd::getDD(*op, *dd)); + deterministicNoiseFunctionality.applyNoiseEffects(rootEdge, op); + } + + auto tmp = rootEdge.getSparseProbabilityVectorStrKeys(qc.getNqubits()); + auto prob = 0.125; + EXPECT_NEAR(tmp["000"], prob, tolerance); + EXPECT_NEAR(tmp["001"], prob, tolerance); + EXPECT_NEAR(tmp["010"], prob, tolerance); + EXPECT_NEAR(tmp["011"], prob, tolerance); + EXPECT_NEAR(tmp["100"], prob, tolerance); + EXPECT_NEAR(tmp["101"], prob, tolerance); + EXPECT_NEAR(tmp["110"], prob, tolerance); + EXPECT_NEAR(tmp["111"], prob, tolerance); + + densityDD.measureOneCollapsing(rootEdge, 0, qc.getGenerator()); + + auto tmp0 = rootEdge.getSparseProbabilityVectorStrKeys(qc.getNqubits()); + prob = 0.25; + + EXPECT_TRUE(std::fabs(tmp0["000"] + tmp0["001"] - prob) < tolerance); + EXPECT_TRUE(std::fabs(tmp0["010"] + tmp0["011"] - prob) < tolerance); + EXPECT_TRUE(std::fabs(tmp0["100"] + tmp0["101"] - prob) < tolerance); + EXPECT_TRUE(std::fabs(tmp0["110"] + tmp0["111"] - prob) < tolerance); + + densityDD.measureOneCollapsing(rootEdge, 1, qc.getGenerator()); + + auto tmp1 = rootEdge.getSparseProbabilityVectorStrKeys(qc.getNqubits()); + prob = 0.5; + EXPECT_TRUE(std::fabs(tmp0["000"] + tmp0["001"] + tmp0["010"] + tmp0["011"] - + prob) < tolerance); + EXPECT_TRUE(std::fabs(tmp0["100"] + tmp0["101"] + tmp0["110"] + tmp0["111"] - + prob) < tolerance); + + densityDD.measureOneCollapsing(rootEdge, 2, qc.getGenerator()); + + auto tmp2 = rootEdge.getSparseProbabilityVectorStrKeys(qc.getNqubits()); + EXPECT_TRUE(std::fabs(tmp2["000"] - 1) < tolerance || + std::fabs(tmp2["001"] - 1) < tolerance || + std::fabs(tmp2["010"] - 1) < tolerance || + std::fabs(tmp2["011"] - 1) < tolerance || + std::fabs(tmp2["100"] - 1) < tolerance || + std::fabs(tmp2["101"] - 1) < tolerance || + std::fabs(tmp2["111"] - 1) < tolerance); +} + +TEST_F(DDNoiseFunctionalityTest, StochSimulateAdder4TrackAPD) { + auto dd = std::make_unique( + qc.getNqubits(), dd::ddsim::STOCHASTIC_NOISE_SIMULATOR_DD_PACKAGE_CONFIG); + + std::map> measSummary = { + {"0000", 0.}, {"0001", 0.}, {"0010", 0.}, {"0011", 0.}, {"0100", 0.}, + {"0101", 0.}, {"0110", 0.}, {"0111", 0.}, {"1000", 0.}, {"1001", 0.}, + {"1010", 0.}, {"1011", 0.}, {"1100", 0.}, {"1101", 0.}}; + + const auto* const noiseEffects = "APDI"; + + auto stochasticNoiseFunctionality = dd::ddsim::StochasticNoiseFunctionality( + *dd, qc.getNqubits(), 0.01, 0.02, 2., noiseEffects); + + for (std::size_t i = 0U; i < stochRuns; i++) { + auto rootEdge = dd::makeZeroState(qc.getNqubits(), *dd); + dd->incRef(rootEdge); + + for (auto const& op : qc) { + auto operation = dd::getDD(*op, *dd); + auto usedQubits = op->getUsedQubits(); + stochasticNoiseFunctionality.applyNoiseOperation( + usedQubits, operation, rootEdge, qc.getGenerator()); + } + + const auto amplitudes = rootEdge.getVector(); + for (std::size_t m = 0U; m < amplitudes.size(); m++) { + auto state = std::bitset<4U>(m).to_string(); + std::ranges::reverse(state); + const auto amplitude = amplitudes[m]; + const auto prob = std::norm(amplitude); + measSummary[state] += prob / static_cast(stochRuns); + } + } + + const double tolerance = 0.1; + EXPECT_NEAR(measSummary["0000"], 0.09693321927412533, tolerance); + EXPECT_NEAR(measSummary["0001"], 0.09078880415385877, tolerance); + EXPECT_NEAR(measSummary["0010"], 0.01414096609854787, tolerance); + EXPECT_NEAR(measSummary["0100"], 0.02382034755245074, tolerance); + EXPECT_NEAR(measSummary["0101"], 0.023509799001774703, tolerance); + EXPECT_NEAR(measSummary["0110"], 0.02445760874001203, tolerance); + EXPECT_NEAR(measSummary["0111"], 0.011628281127642115, tolerance); + EXPECT_NEAR(measSummary["1000"], 0.1731941264570172, tolerance); + EXPECT_NEAR(measSummary["1001"], 0.41458550719988047, tolerance); + EXPECT_NEAR(measSummary["1010"], 0.013806211321349706, tolerance); + EXPECT_NEAR(measSummary["1011"], 0.01840334820660922, tolerance); + EXPECT_NEAR(measSummary["1100"], 0.024245433691737584, tolerance); + EXPECT_NEAR(measSummary["1101"], 0.026277984479993615, tolerance); + EXPECT_NEAR(measSummary["1110"], 0.023929692098939092, tolerance); + EXPECT_NEAR(measSummary["1111"], 0.011037316662706232, tolerance); +} + +TEST_F(DDNoiseFunctionalityTest, StochSimulateAdder4IdentityError) { + auto dd = std::make_unique( + qc.getNqubits(), dd::ddsim::STOCHASTIC_NOISE_SIMULATOR_DD_PACKAGE_CONFIG); + + std::map> measSummary = { + {"0000", 0.}, {"0001", 0.}, {"0010", 0.}, {"0011", 0.}, {"0100", 0.}, + {"0101", 0.}, {"0110", 0.}, {"0111", 0.}, {"1000", 0.}, {"1001", 0.}, + {"1010", 0.}, {"1011", 0.}, {"1100", 0.}, {"1101", 0.}}; + + const auto* const noiseEffects = "I"; + + auto stochasticNoiseFunctionality = dd::ddsim::StochasticNoiseFunctionality( + *dd, qc.getNqubits(), 0.01, 0.02, 2., noiseEffects); + + for (std::size_t i = 0U; i < stochRuns; i++) { + auto rootEdge = dd::makeZeroState(qc.getNqubits(), *dd); + dd->incRef(rootEdge); + + for (auto const& op : qc) { + auto operation = dd::getDD(*op, *dd); + stochasticNoiseFunctionality.applyNoiseOperation( + op->getUsedQubits(), operation, rootEdge, qc.getGenerator()); + } + + const auto amplitudes = rootEdge.getVector(); + for (std::size_t m = 0U; m < amplitudes.size(); m++) { + auto state = std::bitset<4U>(m).to_string(); + std::ranges::reverse(state); + const auto amplitude = amplitudes[m]; + const auto prob = std::norm(amplitude); + measSummary[state] += prob / static_cast(stochRuns); + } + } + + const double tolerance = 0.1; + EXPECT_NEAR(measSummary["0000"], 0., tolerance); + EXPECT_NEAR(measSummary["0001"], 0., tolerance); + EXPECT_NEAR(measSummary["0010"], 0., tolerance); + EXPECT_NEAR(measSummary["0100"], 0., tolerance); + EXPECT_NEAR(measSummary["0101"], 0., tolerance); + EXPECT_NEAR(measSummary["0110"], 0., tolerance); + EXPECT_NEAR(measSummary["0111"], 0., tolerance); + EXPECT_NEAR(measSummary["1000"], 0., tolerance); + EXPECT_NEAR(measSummary["1001"], 1., tolerance); + EXPECT_NEAR(measSummary["1010"], 0., tolerance); + EXPECT_NEAR(measSummary["1011"], 0., tolerance); + EXPECT_NEAR(measSummary["1100"], 0., tolerance); + EXPECT_NEAR(measSummary["1101"], 0., tolerance); + EXPECT_NEAR(measSummary["1110"], 0., tolerance); + EXPECT_NEAR(measSummary["1111"], 0., tolerance); +} + +TEST_F(DDNoiseFunctionalityTest, invalidNoiseEffect) { + auto dd = std::make_unique( + qc.getNqubits(), dd::ddsim::STOCHASTIC_NOISE_SIMULATOR_DD_PACKAGE_CONFIG); + EXPECT_THROW(dd::ddsim::StochasticNoiseFunctionality(*dd, qc.getNqubits(), + 0.01, 0.02, 2., "APK"), + std::runtime_error); +} + +TEST_F(DDNoiseFunctionalityTest, invalidNoiseProbabilities) { + auto dd = std::make_unique( + qc.getNqubits(), dd::ddsim::STOCHASTIC_NOISE_SIMULATOR_DD_PACKAGE_CONFIG); + EXPECT_THROW(dd::ddsim::StochasticNoiseFunctionality(*dd, qc.getNqubits(), + 0.3, 0.6, 2, "APD"), + std::runtime_error); +} From 6eb4c010434017886e5dfa4d52ba0ec5f75e357a Mon Sep 17 00:00:00 2001 From: Daniel Haag <121057143+denialhaag@users.noreply.github.com> Date: Tue, 4 Aug 2026 21:26:11 +0200 Subject: [PATCH 5/5] Update upgrade guide Assisted-by: Claude Opus 4.8 via Claude Code --- UPGRADING.md | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/UPGRADING.md b/UPGRADING.md index df0c34d9e..9a93c3dbd 100644 --- a/UPGRADING.md +++ b/UPGRADING.md @@ -8,6 +8,16 @@ of changes including minor and patch releases, please refer to the This release updates the minimum required `mqt-core` version to 3.8.0. +### Vendored density-matrix decision-diagram support + +Support for density-matrix decision diagrams will be removed from `mqt-core` in +version 4.0.0. Because the deterministic and stochastic noise-aware simulators +depend on it, MQT DDSIM now vendors this functionality under the new `dd::ddsim` +namespace. This is an internal change: the `DeterministicNoiseSimulator` and +`StochasticNoiseSimulator` APIs are unaffected. If your C++ code previously used +the density-matrix DD types (e.g. `dd::dEdge` or `dd::DensityMatrixDD`) through +MQT DDSIM, switch to their `dd::ddsim::` counterparts. + ## [2.4.0] This release updates the minimum required `mqt-core` version to 3.7.0 as well as