From d6a976228c64fc74689e5ebe44aa8a90af133f34 Mon Sep 17 00:00:00 2001 From: burgholzer Date: Sat, 17 Jan 2026 00:20:11 +0100 Subject: [PATCH 1/2] =?UTF-8?q?=F0=9F=94=A5=20Remove=20density=20matrix=20?= =?UTF-8?q?support=20from=20DD=20package?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: burgholzer --- .clang-tidy | 1 + CHANGELOG.md | 2 + UPGRADING.md | 7 + bindings/dd/register_matrix_dds.cpp | 10 +- include/mqt-core/dd/Approximation.hpp | 3 + include/mqt-core/dd/CachedEdge.hpp | 33 +- include/mqt-core/dd/Complex.hpp | 4 +- include/mqt-core/dd/ComputeTable.hpp | 32 +- include/mqt-core/dd/DDDefinitions.hpp | 16 +- include/mqt-core/dd/DDpackageConfig.hpp | 52 -- include/mqt-core/dd/DensityNoiseTable.hpp | 114 ----- include/mqt-core/dd/Edge.hpp | 172 ++----- include/mqt-core/dd/MemoryManager.hpp | 8 +- include/mqt-core/dd/Node.hpp | 115 +---- include/mqt-core/dd/NoiseFunctionality.hpp | 144 ------ include/mqt-core/dd/Package.hpp | 383 +++----------- include/mqt-core/dd/Package_fwd.hpp | 1 - include/mqt-core/dd/RealNumber.hpp | 2 +- include/mqt-core/dd/StateGeneration.hpp | 3 +- .../dd/StochasticNoiseOperationTable.hpp | 88 ---- include/mqt-core/dd/UniqueTable.hpp | 13 +- .../dd/statistics/PackageStatistics.hpp | 4 +- .../dd/statistics/TableStatistics.hpp | 2 +- .../dd/statistics/UniqueTableStatistics.hpp | 2 +- include/mqt-core/ir/operations/OpType.inc | 6 - src/dd/CachedEdge.cpp | 27 +- src/dd/Edge.cpp | 269 ++-------- src/dd/Node.cpp | 106 ---- src/dd/NoiseFunctionality.cpp | 469 ------------------ src/dd/Package.cpp | 163 +++--- src/dd/StateGeneration.cpp | 7 +- src/dd/statistics/PackageStatistics.cpp | 64 +-- test/dd/test_approximations.cpp | 18 +- test/dd/test_dd_functionality.cpp | 63 ++- test/dd/test_dd_noise_functionality.cpp | 344 ------------- test/dd/test_edge_functionality.cpp | 115 ----- test/dd/test_package.cpp | 235 +-------- test/dd/test_state_generation.cpp | 35 +- 38 files changed, 368 insertions(+), 2764 deletions(-) delete mode 100644 include/mqt-core/dd/DensityNoiseTable.hpp delete mode 100644 include/mqt-core/dd/NoiseFunctionality.hpp delete mode 100644 include/mqt-core/dd/StochasticNoiseOperationTable.hpp delete mode 100644 src/dd/Node.cpp delete mode 100644 src/dd/NoiseFunctionality.cpp delete mode 100644 test/dd/test_dd_noise_functionality.cpp diff --git a/.clang-tidy b/.clang-tidy index 8683b18a06..3d9e3b1bfd 100644 --- a/.clang-tidy +++ b/.clang-tidy @@ -14,6 +14,7 @@ Checks: | -cppcoreguidelines-macro-usage, -cppcoreguidelines-pro-type-reinterpret-cast, -cppcoreguidelines-pro-bounds-constant-array-index, + -cppcoreguidelines-pro-bounds-avoid-unchecked-container-access, -cppcoreguidelines-avoid-do-while, google-*, -google-readability-todo, diff --git a/CHANGELOG.md b/CHANGELOG.md index 5f7a707a48..938f4ada36 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -25,6 +25,7 @@ This project adheres to [Semantic Versioning], with the exception that minor rel ### Removed +- 🔥 Remove the density matrix support from the MQT Core DD package ([#1466]) ([**@burgholzer**]) - 🔥 Remove `datastructures` (`ds`) (sub)library from MQT Core ([#1458]) ([**@burgholzer**]) - 🔥 No longer actively type check Python code with `mypy` and solely rely on `ty` ([#1437]) ([**@burgholzer**]) @@ -314,6 +315,7 @@ _📚 Refer to the [GitHub Release Notes](https://github.com/munich-quantum-tool +[#1466]: https://github.com/munich-quantum-toolkit/core/pull/1466 [#1465]: https://github.com/munich-quantum-toolkit/core/pull/1465 [#1458]: https://github.com/munich-quantum-toolkit/core/pull/1458 [#1453]: https://github.com/munich-quantum-toolkit/core/pull/1453 diff --git a/UPGRADING.md b/UPGRADING.md index 4d8b72e206..5589af48b9 100644 --- a/UPGRADING.md +++ b/UPGRADING.md @@ -24,6 +24,13 @@ Known limitations: - AppleClang 17+ is required to build MQT Core with MLIR enabled due to some C++20 features being used that are not yet properly supported by older versions. - Our pre-built distributions are compiled in Release mode. On Windows, this leads to ABI incompatibilities with debug builds. Either build in Release mode or build LLVM from source in Debug mode to resolve this. +### Removal of the density matrix support from the DD package + +The density matrix support within the DD package has been removed. +This change was made to reduce the maintenance burden of the package. +Any libraries that depend on the density matrix functionality, such as [MQT DDSIM], need to implement it on their own or use an alternative solution. +In a related fashion, this PR also removes the noise operations from the MQT Core IR as they no longer serve a purpose. + ### Removal of the `datastructures` (sub)library The `datastructures` (sub)library has been removed from the MQT Core repository. diff --git a/bindings/dd/register_matrix_dds.cpp b/bindings/dd/register_matrix_dds.cpp index 259fad608d..81dd70187f 100644 --- a/bindings/dd/register_matrix_dds.cpp +++ b/bindings/dd/register_matrix_dds.cpp @@ -74,8 +74,7 @@ void registerMatrixDDs(const nb::module_& m) { mat.def("is_zero_terminal", &dd::mEdge::isZeroTerminal, "Check if the DD is a zero terminal node."); - mat.def("is_identity", &dd::mEdge::isIdentity<>, - "up_to_global_phase"_a = true, + mat.def("is_identity", &dd::mEdge::isIdentity, "up_to_global_phase"_a = true, R"pb(Check if the DD represents the identity matrix. Args: @@ -87,8 +86,11 @@ void registerMatrixDDs(const nb::module_& m) { mat.def("size", nb::overload_cast<>(&dd::mEdge::size, nb::const_), "Get the size of the DD by traversing it once."); - mat.def("get_entry", &dd::mEdge::getValueByIndex<>, "num_qubits"_a, "row"_a, - "col"_a, "Get the entry of the matrix by row and column index."); + mat.def("get_entry", + nb::overload_cast(&dd::mEdge::getValueByIndex, + nb::const_), + "num_qubits"_a, "row"_a, "col"_a, + "Get the entry of the matrix by row and column index."); mat.def("get_entry_by_path", &dd::mEdge::getValueByPath, "num_qubits"_a, "decisions"_a, R"pb(Get the entry of the matrix by decisions. diff --git a/include/mqt-core/dd/Approximation.hpp b/include/mqt-core/dd/Approximation.hpp index 9ca955f264..c76823f1e8 100644 --- a/include/mqt-core/dd/Approximation.hpp +++ b/include/mqt-core/dd/Approximation.hpp @@ -10,9 +10,12 @@ #pragma once +#include "dd/DDDefinitions.hpp" #include "dd/Node.hpp" #include "dd/Package.hpp" +#include + namespace dd { /** diff --git a/include/mqt-core/dd/CachedEdge.hpp b/include/mqt-core/dd/CachedEdge.hpp index 5527718f8a..95db0c1cae 100644 --- a/include/mqt-core/dd/CachedEdge.hpp +++ b/include/mqt-core/dd/CachedEdge.hpp @@ -18,25 +18,12 @@ #include #include #include -#include namespace dd { -struct vNode; // NOLINT(readability-identifier-naming) -struct mNode; // NOLINT(readability-identifier-naming) -struct dNode; // NOLINT(readability-identifier-naming) class ComplexNumbers; class MemoryManager; -template -using isVector = std::enable_if_t, bool>; -template -using isMatrix = std::enable_if_t, bool>; -template -using isMatrixVariant = - std::enable_if_t || std::is_same_v, - bool>; - /** * @brief A DD node with a cached edge weight * @details Some DD operations create intermediate results that are not part of @@ -116,7 +103,6 @@ template struct CachedEdge { /** * @brief Get a normalized vector DD from a fresh node and a list of edges. - * @tparam T template parameter to enable this method only for vNode * @param p the fresh node * @param e the list of edges that form the successor nodes * @param mm a reference to the memory manager (for returning unused nodes) @@ -124,32 +110,31 @@ template struct CachedEdge { * complex numbers) * @return the normalized vector DD */ - template = true> static auto normalize(Node* p, const std::array& e, - MemoryManager& mm, ComplexNumbers& cn) -> CachedEdge; + MemoryManager& mm, ComplexNumbers& cn) -> CachedEdge + requires IsVector; /** - * @brief Get a normalized (density) matrix) DD from a fresh node and a list + * @brief Get a normalized matrix DD from a fresh node and a list * of edges. - * @tparam T template parameter to enable this method only for matrix nodes * @param p the fresh node * @param e the list of edges that form the successor nodes * @param mm a reference to the memory manager (for returning unused nodes) * @param cn a reference to the complex number manager (for adding new * complex numbers) - * @return the normalized (density) matrix DD + * @return the normalized matrix DD */ - template = true> static auto normalize(Node* p, const std::array& e, - MemoryManager& mm, ComplexNumbers& cn) -> CachedEdge; + MemoryManager& mm, ComplexNumbers& cn) -> CachedEdge + requires IsMatrix; /** * @brief Check whether the matrix represented by the DD is the identity. - * @tparam T template parameter to enable this function only for matrix nodes * @return whether the matrix is the identity */ - template = true> - [[nodiscard]] bool isIdentity(const bool upToGlobalPhase = true) const { + [[nodiscard]] bool isIdentity(const bool upToGlobalPhase = true) const + requires IsMatrix + { if (!isTerminal()) { return false; } diff --git a/include/mqt-core/dd/Complex.hpp b/include/mqt-core/dd/Complex.hpp index 2764529230..d1ec6f108c 100644 --- a/include/mqt-core/dd/Complex.hpp +++ b/include/mqt-core/dd/Complex.hpp @@ -37,7 +37,7 @@ struct Complex { * @return A complex number with real and imaginary part equal to zero. */ static constexpr Complex zero() noexcept { - return {&constants::zero, &constants::zero}; + return {.r = &constants::zero, .i = &constants::zero}; } /** @@ -46,7 +46,7 @@ struct Complex { * equal to zero. */ static constexpr Complex one() noexcept { - return {&constants::one, &constants::zero}; + return {.r = &constants::one, .i = &constants::zero}; } /** diff --git a/include/mqt-core/dd/ComputeTable.hpp b/include/mqt-core/dd/ComputeTable.hpp index 243afce601..4f8a46c459 100644 --- a/include/mqt-core/dd/ComputeTable.hpp +++ b/include/mqt-core/dd/ComputeTable.hpp @@ -15,7 +15,6 @@ #pragma once -#include "dd/Node.hpp" #include "dd/statistics/TableStatistics.hpp" #include "ir/Definitions.hpp" @@ -73,20 +72,8 @@ class ComputeTable { */ [[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, dd::dNode::getDensityMatrixTempFlags(leftOperand->flags)); - } - } - auto h2 = std::hash{}(rightOperand); - if constexpr (std::is_same_v) { - if (!dNode::isTerminal(rightOperand)) { - h2 = qc::combineHash( - h2, dd::dNode::getDensityMatrixTempFlags(rightOperand->flags)); - } - } + const auto h1 = std::hash{}(leftOperand); + const auto h2 = std::hash{}(rightOperand); const auto hash = qc::combineHash(h1, h2); const auto mask = stats.numBuckets - 1; return hash & mask; @@ -121,12 +108,10 @@ class ComputeTable { * @brief Look up a result in the compute table * @param leftOperand The left operand * @param rightOperand The right operand - * @param useDensityMatrix Whether a density matrix is expected * @return A pointer to the result if it is found, otherwise nullptr. */ ResultType* lookup(const LeftOperandType& leftOperand, - const RightOperandType& rightOperand, - [[maybe_unused]] const bool useDensityMatrix = false) { + const RightOperandType& rightOperand) { ResultType* result = nullptr; ++stats.lookups; const auto key = hash(leftOperand, rightOperand); @@ -142,17 +127,6 @@ class ComputeTable { 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; } diff --git a/include/mqt-core/dd/DDDefinitions.hpp b/include/mqt-core/dd/DDDefinitions.hpp index 0b334189b5..da936fa2d6 100644 --- a/include/mqt-core/dd/DDDefinitions.hpp +++ b/include/mqt-core/dd/DDDefinitions.hpp @@ -17,6 +17,7 @@ #include #include #include +#include #include #include #include @@ -57,8 +58,7 @@ enum class BasisStates : std::uint8_t { static constexpr fp SQRT2_2 = static_cast( 0.707106781186547524400844362104849039284835937688474036588L); -static constexpr fp PI = static_cast( - 3.141592653589793238462643383279502884197169399375105820974L); +static constexpr fp PI = std::numbers::pi; static constexpr fp PI_2 = static_cast( 1.570796326794896619231321691639751442098584699687552910487L); static constexpr fp PI_4 = static_cast( @@ -91,7 +91,7 @@ using TwoQubitGateMatrix = * @param nbits The number of bits to use for the binary representation * @return The binary representation of the decimal number */ -[[nodiscard]] static inline std::string +[[nodiscard, maybe_unused]] static std::string intToBinaryString(const std::size_t value, const std::size_t nbits) { std::string binary(nbits, '0'); for (std::size_t j = 0; j < nbits; ++j) { @@ -104,7 +104,7 @@ intToBinaryString(const std::size_t value, const std::size_t nbits) { // calculates the Units in Last Place (ULP) distance of two floating point // numbers -[[maybe_unused]] static std::size_t ulpDistance(fp a, fp b) { +[[nodiscard, maybe_unused]] static std::size_t ulpDistance(fp a, fp b) { // NOLINTNEXTLINE(clang-diagnostic-float-equal) if (a == b) { return 0; @@ -136,4 +136,12 @@ intToBinaryString(const std::size_t value, const std::size_t nbits) { return k; } +struct vNode; +struct mNode; + +template +concept IsVector = std::is_same_v; +template +concept IsMatrix = std::is_same_v; + } // namespace dd diff --git a/include/mqt-core/dd/DDpackageConfig.hpp b/include/mqt-core/dd/DDpackageConfig.hpp index fdf5e9b99c..f29cc0579e 100644 --- a/include/mqt-core/dd/DDpackageConfig.hpp +++ b/include/mqt-core/dd/DDpackageConfig.hpp @@ -10,8 +10,6 @@ #pragma once -#include "ir/operations/OpType.hpp" - #include namespace dd { @@ -30,58 +28,10 @@ struct DDPackageConfig { std::size_t ctMatMatMultNumBucket = 16384U; std::size_t ctVecKronNumBucket = 4096U; std::size_t ctMatKronNumBucket = 4096U; - std::size_t ctDmTraceNumBucket = 1U; std::size_t ctMatTraceNumBucket = 4096U; std::size_t ctVecInnerProdNumBucket = 4096U; - std::size_t ctDmNoiseNumBucket = 1U; - std::size_t utDmNumBucket = 1U; - std::size_t utDmInitialAllocationSize = 1U; - std::size_t ctDmDmMultNumBucket = 1U; - std::size_t ctDmAddNumBucket = 1U; - - // The number of different quantum operations. i.e., the number of operations - // defined in OpType.hpp. This parameter is required to initialize the - // StochasticNoiseOperationTable.hpp - std::size_t stochasticCacheOps = 1; }; -constexpr auto STOCHASTIC_NOISE_SIMULATOR_DD_PACKAGE_CONFIG = []() { - DDPackageConfig config{}; - config.stochasticCacheOps = qc::OpType::OpTypeEnd; - config.ctVecAddMagNumBucket = 1U; - config.ctMatAddMagNumBucket = 1U; - config.ctVecConjNumBucket = 1U; - return config; -}(); - -constexpr auto DENSITY_MATRIX_SIMULATOR_DD_PACKAGE_CONFIG = []() { - DDPackageConfig config{}; - config.utDmNumBucket = 65536U; - config.utDmInitialAllocationSize = 4096U; - config.ctDmDmMultNumBucket = 16384U; - config.ctDmAddNumBucket = 16384U; - config.ctDmNoiseNumBucket = 4096U; - 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.ctDmTraceNumBucket = 4096U; - config.ctMatTraceNumBucket = 1U; - config.ctVecInnerProdNumBucket = 1U; - config.stochasticCacheOps = 1U; - config.ctVecAddMagNumBucket = 1U; - config.ctMatAddMagNumBucket = 1U; - config.ctVecConjNumBucket = 1U; - return config; -}(); - constexpr auto UNITARY_SIMULATOR_DD_PACKAGE_CONFIG = []() { DDPackageConfig config{}; config.utMatNumBucket = 65'536U; @@ -90,8 +40,6 @@ constexpr auto UNITARY_SIMULATOR_DD_PACKAGE_CONFIG = []() { config.utVecNumBucket = 1U; config.utVecInitialAllocationSize = 1U; config.ctVecAddNumBucket = 1U; - config.ctVecAddMagNumBucket = 1U; - config.ctMatAddMagNumBucket = 1U; config.ctVecConjNumBucket = 1U; config.ctMatConjTransNumBucket = 1U; config.ctMatVecMultNumBucket = 1U; diff --git a/include/mqt-core/dd/DensityNoiseTable.hpp b/include/mqt-core/dd/DensityNoiseTable.hpp deleted file mode 100644 index 506827e8cd..0000000000 --- a/include/mqt-core/dd/DensityNoiseTable.hpp +++ /dev/null @@ -1,114 +0,0 @@ -/* - * 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 - */ - -#pragma once - -#include "dd/DDDefinitions.hpp" -#include "dd/statistics/TableStatistics.hpp" - -#include -#include -#include -#include - -namespace dd { - -/** - * @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 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 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; - TableStatistics stats{}; -}; -} // namespace dd diff --git a/include/mqt-core/dd/Edge.hpp b/include/mqt-core/dd/Edge.hpp index 543f9cbe9b..7a44f96a35 100644 --- a/include/mqt-core/dd/Edge.hpp +++ b/include/mqt-core/dd/Edge.hpp @@ -12,34 +12,20 @@ #include "dd/Complex.hpp" #include "dd/DDDefinitions.hpp" +#include "dd/RealNumber.hpp" #include #include #include #include #include -#include #include namespace dd { -struct vNode; -struct mNode; -struct dNode; class ComplexNumbers; class MemoryManager; -template -using isVector = std::enable_if_t, bool>; -template -using isMatrix = std::enable_if_t, bool>; -template -using isDensityMatrix = std::enable_if_t, bool>; -template -using isMatrixVariant = - std::enable_if_t || std::is_same_v, - bool>; - using AmplitudeFunc = std::function&)>; using ProbabilityFunc = std::function; using MatrixEntryFunc = @@ -159,7 +145,6 @@ template struct Edge { public: /** * @brief Get a normalized vector DD from a fresh node and a list of edges - * @tparam T template parameter to enable this function only for vNode * @param p the fresh node * @param e the list of edges that form the successor nodes * @param mm a reference to the memory manager (for returning unused nodes) @@ -167,60 +152,55 @@ template struct Edge { * complex numbers) * @return the normalized vector DD */ - template = true> static auto normalize(Node* p, const std::array& e, - MemoryManager& mm, ComplexNumbers& cn) -> Edge; + MemoryManager& mm, ComplexNumbers& cn) -> Edge + requires IsVector; /** * @brief Get a single element of the vector represented by the DD - * @tparam T template parameter to enable this function only for vNode * @param i index of the element * @return the complex value of the amplitude */ - template = true> - [[nodiscard]] std::complex getValueByIndex(std::size_t i) const; + [[nodiscard]] std::complex getValueByIndex(std::size_t i) const + requires IsVector; /** * @brief Get the vector represented by the DD - * @tparam T template parameter to enable this function only for vNode * @param threshold amplitudes with a magnitude below this threshold will be * ignored * @return the vector */ - template = true> - [[nodiscard]] CVec getVector(fp threshold = 0.) const; + [[nodiscard]] CVec getVector(fp threshold = 0.) const + requires IsVector; /** * @brief Get the sparse vector represented by the DD - * @tparam T template parameter to enable this function only for vNode * @param threshold amplitudes with a magnitude below this threshold will be * ignored * @return the sparse vector */ - template = true> - [[nodiscard]] SparseCVec getSparseVector(fp threshold = 0.) const; + [[nodiscard]] SparseCVec getSparseVector(fp threshold = 0.) const + requires IsVector; /** * @brief Print the vector represented by the DD - * @tparam T template parameter to enable this function only for vNode * @note This function scales exponentially with the number of qubits. */ - template = true> void printVector() const; + void printVector() const + requires IsVector; /** * @brief Add the amplitudes of a vector DD to a vector - * @tparam T template parameter to enable this function only for vNode * @param amplitudes the vector to add to */ - template = true> - void addToVector(CVec& amplitudes) const; + void addToVector(CVec& amplitudes) const + requires IsVector; private: /** * @brief Recursively traverse the DD and call a function for each non-zero * amplitude. * @details Scales with the number of non-zero amplitudes. - * @tparam T template parameter to enable this function only for vNode * @param amp the accumulated amplitude from previous traversals * @param i the current index in the vector * @param f This function is called for each non-zero amplitude with the @@ -228,33 +208,32 @@ template struct Edge { * @param threshold amplitude with a magnitude below this threshold will be * ignored */ - template = true> void traverseVector(const std::complex& amp, std::size_t i, - AmplitudeFunc f, fp threshold = 0.) const; + AmplitudeFunc f, fp threshold = 0.) const + requires IsVector; public: /** - * @brief Get a normalized (density) matrix DD from a fresh node and a list + * @brief Get a normalized matrix DD from a fresh node and a list * of edges - * @tparam T template parameter to enable this function only for matrix nodes * @param p the fresh node * @param e the list of edges that form the successor nodes * @param mm a reference to the memory manager (for returning unused nodes) * @param cn a reference to the complex number manager (for adding new * complex numbers) - * @return the normalized (density) matrix DD + * @return the normalized matrix DD */ - template = true> static auto normalize(Node* p, const std::array& e, - MemoryManager& mm, ComplexNumbers& cn) -> Edge; + MemoryManager& mm, ComplexNumbers& cn) -> Edge + requires IsMatrix; /** * @brief Check whether the matrix represented by the DD is the identity - * @tparam T template parameter to enable this function only for matrix nodes * @return whether the matrix is the identity */ - template = true> - [[nodiscard]] bool isIdentity(const bool upToGlobalPhase = true) const { + [[nodiscard]] bool isIdentity(const bool upToGlobalPhase = true) const + requires IsMatrix + { if (!isTerminal()) { return false; } @@ -266,52 +245,47 @@ template struct Edge { /** * @brief Get a single element of the matrix represented by the DD - * @tparam T template parameter to enable this function only for matrix nodes * @param numQubits number of qubits in the considered DD * @param i row index of the element * @param j column index of the element * @return the complex value of the entry */ - template = true> [[nodiscard]] std::complex - getValueByIndex(std::size_t numQubits, std::size_t i, std::size_t j) const; + getValueByIndex(std::size_t numQubits, std::size_t i, std::size_t j) const + requires IsMatrix; /** * @brief Get the matrix represented by the DD - * @tparam T template parameter to enable this function only for matrix nodes * @param numQubits number of qubits in the considered DD * @param threshold entries with a magnitude below this threshold will be * ignored * @return the matrix */ - template = true> - [[nodiscard]] CMat getMatrix(std::size_t numQubits, fp threshold = 0.) const; + [[nodiscard]] CMat getMatrix(std::size_t numQubits, fp threshold = 0.) const + requires IsMatrix; /** * @brief Get the sparse matrix represented by the DD - * @tparam T template parameter to enable this function only for matrix nodes * @param numQubits number of qubits in the considered DD * @param threshold entries with a magnitude below this threshold will be * ignored * @return the sparse matrix */ - template = true> [[nodiscard]] SparseCMat getSparseMatrix(std::size_t numQubits, - fp threshold = 0.) const; + fp threshold = 0.) const + requires IsMatrix; /** * @brief Print the matrix represented by the DD - * @tparam T template parameter to enable this function only for matrix nodes * @param numQubits number of qubits in the considered DD * @note This function scales exponentially with the number of qubits. */ - template = true> - void printMatrix(std::size_t numQubits) const; + void printMatrix(std::size_t numQubits) const + requires IsMatrix; /** * @brief Recursively traverse the DD and call a function for each non-zero * matrix entry. - * @tparam T template parameter to enable this function only for matrix nodes * @param amp the accumulated amplitude from previous traversals * @param i the current row index in the matrix * @param j the current column index in the matrix @@ -322,92 +296,10 @@ template struct Edge { * @param threshold entries with a magnitude below this threshold will be * ignored */ - template = true> void traverseMatrix(const std::complex& amp, std::size_t i, std::size_t j, MatrixEntryFunc f, std::size_t level, - fp threshold = 0.) const; - - template = true> - [[maybe_unused]] static void setDensityConjugateTrue(Edge& e) { - Node::setConjugateTempFlagTrue(e.p); - } - - template = true> - [[maybe_unused]] static void setFirstEdgeDensityPathTrue(Edge& e) { - Node::setNonReduceTempFlagTrue(e.p); - } - - template = true> - static void setDensityMatrixTrue(Edge& e) { - Node::setDensityMatTempFlagTrue(e.p); - } - - template = true> - static void alignDensityEdge(Edge& e) { - Node::alignDensityNode(e.p); - } - - template = true> - static void revertDmChangesToEdges(Edge& x, Edge& y) { - revertDmChangesToEdge(x); - revertDmChangesToEdge(y); - } - - template = true> - static void revertDmChangesToEdge(Edge& x) { - Node::revertDmChangesToNode(x.p); - } - - template = true> - static void applyDmChangesToEdges(Edge& x, Edge& y) { - applyDmChangesToEdge(x); - applyDmChangesToEdge(y); - } - - template = true> - static void applyDmChangesToEdge(Edge& x) { - Node::applyDmChangesToNode(x.p); - } - - /** - * @brief Get the sparse probability vector for the underlying density matrix - * @tparam T template parameter to enable this function only for dNode - * @param numQubits number of qubits in the considered DD - * @param threshold probabilities below this threshold will be ignored - * @return the sparse probability vector - */ - template = true> - [[nodiscard]] SparsePVec getSparseProbabilityVector(std::size_t numQubits, - fp threshold = 0.) const; - - /** - * @brief Get the sparse probability vector for the underlying density matrix - * @tparam T template parameter to enable this function only for dNode - * @param numQubits number of qubits in the considered DD - * @param threshold probabilities below this threshold will be ignored - * @return the sparse probability vector (using strings as keys) - */ - template = true> - [[nodiscard]] SparsePVecStrKeys - getSparseProbabilityVectorStrKeys(std::size_t numQubits, - fp threshold = 0.) const; - -private: - /** - * @brief Recursively traverse diagonal of the DD and call a function for each - * non-zero entry. - * @tparam T template parameter to enable this function only for dNode - * @param prob the accumulated probability from previous traversals - * @param i the current diagonal index in the matrix - * @param f This function is called for each non-zero entry with the - * diagonal index and the probability as arguments. - * @param level the current level in the DD (ranges from 1 to n for regular - * nodes and is 0 for the terminal node) - * @param threshold probabilities below this threshold will be ignored - */ - template = true> - void traverseDiagonal(const fp& prob, std::size_t i, ProbabilityFunc f, - std::size_t level, fp threshold = 0.) const; + fp threshold = 0.) const + requires IsMatrix; }; } // namespace dd diff --git a/include/mqt-core/dd/MemoryManager.hpp b/include/mqt-core/dd/MemoryManager.hpp index da2b59015f..68fca18b31 100644 --- a/include/mqt-core/dd/MemoryManager.hpp +++ b/include/mqt-core/dd/MemoryManager.hpp @@ -13,8 +13,8 @@ #include "dd/statistics/MemoryManagerStatistics.hpp" #include +#include #include -#include #include namespace dd { @@ -86,9 +86,9 @@ class MemoryManager { * @tparam T The type of the entry. * @return A pointer to an entry. */ - template [[nodiscard]] T* get() { - static_assert(std::is_base_of_v, - "T must be derived from LLBase"); + template + requires std::derived_from + [[nodiscard]] T* get() { assert(sizeof(T) == entrySize_ && "Cannot get entry of different size"); return static_cast(get()); diff --git a/include/mqt-core/dd/Node.hpp b/include/mqt-core/dd/Node.hpp index 859a11d302..10d6b94288 100644 --- a/include/mqt-core/dd/Node.hpp +++ b/include/mqt-core/dd/Node.hpp @@ -37,16 +37,12 @@ struct NodeBase : LLBase { * @details Not required for all node types, but padding is required either * way. * - * 0b10000 = mark flag used for mark-and-sweep garbage collection, - * 0b1000 = marks a reduced dm node, - * 0b100 = marks a dm (tmp flag), - * 0b10 = mark first path edge (tmp flag), - * 0b1 = mark path is conjugated (tmp flag) + * 0b1 = mark flag used for mark-and-sweep garbage collection */ std::uint16_t flags = 0; /// Mark flag used for mark-and-sweep garbage collection - static constexpr std::uint16_t MARK_FLAG = 0b10000U; + static constexpr std::uint16_t MARK_FLAG = 0b1U; /// @brief Whether a node is marked as used. [[nodiscard]] bool isMarked() const noexcept { @@ -67,15 +63,11 @@ struct NodeBase : LLBase { /** * @brief Check if a node is terminal - * @details Generally, a node is terminal if it is nullptr. - * Some nodes (dNode) encode additional information in the least significant - * bits of the pointer. These bits are masked out before checking for - * terminal nodes. * @param p The node to check * @return true if the node is terminal, false otherwise. */ - [[nodiscard]] static bool isTerminal(const NodeBase* p) noexcept { - return (reinterpret_cast(p) & (~7ULL)) == 0ULL; + [[nodiscard]] static constexpr bool isTerminal(const NodeBase* p) noexcept { + return p == nullptr; } static constexpr NodeBase* getTerminal() noexcept { return nullptr; } }; @@ -121,103 +113,4 @@ using mEdge = Edge; using mCachedEdge = CachedEdge; using MatrixDD = mEdge; -/** - * @brief A density matrix DD node - * @details Data Layout (8)|(2|2|4)|(24|24|24|24) = 112B - */ -struct dNode final : NodeBase { // NOLINT(readability-identifier-naming) - std::array, NEDGE> 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; -}; -using dEdge = Edge; -using dCachedEdge = CachedEdge; -using DensityMatrixDD = dEdge; - -static inline dEdge densityFromMatrixEdge(const mEdge& e) { - return dEdge{reinterpret_cast(e.p), e.w}; -} - } // namespace dd diff --git a/include/mqt-core/dd/NoiseFunctionality.hpp b/include/mqt-core/dd/NoiseFunctionality.hpp deleted file mode 100644 index d06fc7eace..0000000000 --- a/include/mqt-core/dd/NoiseFunctionality.hpp +++ /dev/null @@ -1,144 +0,0 @@ -/* - * 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 - */ - -#pragma once - -#include "dd/DDDefinitions.hpp" -#include "dd/Node.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 dd { - -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(Package& dd, std::size_t nq, - double gateNoiseProbability, - double amplitudeDampingProb, - double multiQubitGateFactor, - const std::string& cNoiseEffects); - - ~StochasticNoiseFunctionality() { package->decRef(identityDD); } - -protected: - Package* package; - std::size_t nQubits; - std::uniform_real_distribution dist; - - double noiseProbability; - double noiseProbabilityMulti; - fp sqrtAmplitudeDampingProbability; - fp oneMinusSqrtAmplitudeDampingProbability; - fp sqrtAmplitudeDampingProbabilityMulti; - fp oneMinusSqrtAmplitudeDampingProbabilityMulti; - GateMatrix ampDampingTrue{}; - GateMatrix ampDampingTrueMulti{}; - GateMatrix ampDampingFalse{}; - GateMatrix ampDampingFalseMulti{}; - std::vector noiseEffects; - mEdge identityDD; - - [[nodiscard]] std::size_t getNumberOfQubits() const { return nQubits; } - [[nodiscard]] double getNoiseProbability(bool multiQubitNoiseFlag) const; - - [[nodiscard]] static qc::OpType - getAmplitudeDampingOperationType(bool multiQubitNoiseFlag, - bool amplitudeDampingFlag); - - [[nodiscard]] GateMatrix - getAmplitudeDampingOperationMatrix(bool multiQubitNoiseFlag, - bool amplitudeDampingFlag) const; - -public: - void applyNoiseOperation(const std::set& targets, mEdge operation, - vEdge& state, std::mt19937_64& generator); - -protected: - [[nodiscard]] mEdge stackOperation(const mEdge& operation, qc::Qubit target, - qc::OpType noiseOperation, - const GateMatrix& matrix) const; - - mEdge generateNoiseOperation(mEdge operation, qc::Qubit target, - std::mt19937_64& generator, - bool amplitudeDamping, bool multiQubitOperation); - - [[nodiscard]] qc::OpType returnNoiseOperation(NoiseOperations noiseOperation, - double prob, - bool multiQubitNoiseFlag) const; -}; - -class DeterministicNoiseFunctionality { -public: - DeterministicNoiseFunctionality(Package& dd, std::size_t nq, - double noiseProbabilitySingleQubit, - double noiseProbabilityMultiQubit, - double ampDampProbSingleQubit, - double ampDampProbMultiQubit, - const std::string& cNoiseEffects); - -protected: - Package* 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, 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 diff --git a/include/mqt-core/dd/Package.hpp b/include/mqt-core/dd/Package.hpp index 1ad91f8b71..1caf7936ea 100644 --- a/include/mqt-core/dd/Package.hpp +++ b/include/mqt-core/dd/Package.hpp @@ -17,23 +17,19 @@ #include "dd/ComputeTable.hpp" #include "dd/DDDefinitions.hpp" #include "dd/DDpackageConfig.hpp" -#include "dd/DensityNoiseTable.hpp" #include "dd/Edge.hpp" #include "dd/MemoryManager.hpp" #include "dd/Node.hpp" #include "dd/Package_fwd.hpp" // IWYU pragma: export #include "dd/RealNumber.hpp" #include "dd/RealNumberUniqueTable.hpp" -#include "dd/StochasticNoiseOperationTable.hpp" #include "dd/UnaryComputeTable.hpp" #include "dd/UniqueTable.hpp" #include "ir/Definitions.hpp" #include "ir/Permutation.hpp" #include "ir/operations/Control.hpp" -#include #include -#include #include #include #include @@ -41,6 +37,7 @@ #include #include #include +#include #include #include #include @@ -119,9 +116,6 @@ class Package { /// The memory manager for matrix nodes MemoryManager mMemoryManager{ MemoryManager::create(config_.utMatInitialAllocationSize)}; - /// The memory manager for density matrix nodes - MemoryManager dMemoryManager{ - MemoryManager::create(config_.utDmInitialAllocationSize)}; /** * @brief The memory manager for complex numbers * @note The real and imaginary part of complex numbers are treated @@ -139,8 +133,6 @@ class Package { return vMemoryManager; } else if constexpr (std::is_same_v) { return mMemoryManager; - } else if constexpr (std::is_same_v) { - return dMemoryManager; } else if constexpr (std::is_same_v) { return cMemoryManager; } @@ -155,11 +147,11 @@ class Package { void resetMemoryManagers(bool resizeToTotal = false); /// The unique table used for vector nodes - UniqueTable vUniqueTable{vMemoryManager, {0U, config_.utVecNumBucket}}; + UniqueTable vUniqueTable{vMemoryManager, + {.nVars = 0U, .nBuckets = config_.utVecNumBucket}}; /// The unique table used for matrix nodes - UniqueTable mUniqueTable{mMemoryManager, {0U, config_.utMatNumBucket}}; - /// The unique table used for density matrix nodes - UniqueTable dUniqueTable{dMemoryManager, {0U, config_.utDmNumBucket}}; + UniqueTable mUniqueTable{mMemoryManager, + {.nVars = 0U, .nBuckets = config_.utMatNumBucket}}; /** * @brief The unique table used for complex numbers * @note The table actually only stores real numbers in the interval [0, 1], @@ -179,8 +171,6 @@ class Package { return vUniqueTable; } else if constexpr (std::is_same_v) { return mUniqueTable; - } else if constexpr (std::is_same_v) { - return dUniqueTable; } else if constexpr (std::is_same_v) { return cUniqueTable; } @@ -224,13 +214,12 @@ class Package { private: struct RootSetManager { - public: template using RootSet = std::unordered_map, std::size_t>; /// @brief Add to respective root set. template void addToRoots(const Edge& e) noexcept { - getRoots()[e]++; + ++(getRoots()[e]); } /// @brief Remove from respective root set. @@ -257,20 +246,19 @@ class Package { void reset() { vRoots.clear(); mRoots.clear(); - dRoots.clear(); } private: /// @brief Mark edges contained in @p roots. template static void mark(const RootSet& roots) { - for (auto& [edge, _] : roots) { + for (auto& edge : roots | std::views::keys) { edge.mark(); } } /// @brief Unmark edges contained in @p roots. template static void unmark(const RootSet& roots) { - for (auto& [edge, _] : roots) { + for (auto& edge : roots | std::views::keys) { edge.unmark(); } } @@ -279,40 +267,32 @@ class Package { void mark() noexcept { RootSetManager::mark(vRoots); RootSetManager::mark(mRoots); - RootSetManager::mark(dRoots); } /// @brief Unmark edges contained in all root sets. void unmark() noexcept { RootSetManager::unmark(vRoots); RootSetManager::unmark(mRoots); - RootSetManager::unmark(dRoots); } /// @brief Return vector roots. - template , bool> = true> - auto& getRoots() noexcept { + template + auto& getRoots() noexcept + requires(IsVector) + { return vRoots; } /// @brief Return matrix roots. - template , bool> = true> - auto& getRoots() noexcept { + template + auto& getRoots() noexcept + requires(IsMatrix) + { return mRoots; } - /// @brief Return density roots. - template , bool> = true> - auto& getRoots() noexcept { - return dRoots; - } - RootSet vRoots; RootSet mRoots; - RootSet dRoots; template friend auto& Package::getRootSet() noexcept; }; @@ -331,15 +311,14 @@ class Package { * * @param force Force garbage collect, regardless of whether any * table reports that it may need collecting. - * @returns Whether at least one vector, matrix, density-matrix node, or - * any complex number was reclaimed. + * @returns Whether at least one vector, matrix, or any complex number was + * reclaimed. */ bool garbageCollect(bool force = false); struct ActiveCounts { std::size_t vector = 0U; std::size_t matrix = 0U; - std::size_t density = 0U; std::size_t reals = 0U; }; /** @@ -348,18 +327,6 @@ class Package { */ [[nodiscard]] ActiveCounts computeActiveCounts(); - /// - /// Vector nodes, edges and quantum states - /// - - /** - * @brief Construct the all-zero density operator - \f$|0...0\rangle\langle0...0|\f$ - * @param n The number of qubits - * @return A decision diagram for the all-zero density operator - */ - dEdge makeZeroDensityOperator(std::size_t n); - /// /// Matrix nodes, edges and quantum gates /// @@ -484,29 +451,23 @@ class Package { * @tparam EdgeType The type of the edge. * @param var The variable associated with the node. * @param edges The edges of the node. - * @param generateDensityMatrix Flag to indicate if a density matrix node - * should be generated. * @return An edge pointing to the normalized DD node. */ template class EdgeType> EdgeType makeDDNode(const Qubit var, const std::array, - std::tuple_size_v>& edges, - [[maybe_unused]] const bool generateDensityMatrix = false) { + std::tuple_size_v>& edges) { auto& memoryManager = getMemoryManager(); - auto p = memoryManager.template get(); + auto* p = memoryManager.template get(); p->v = var; - if constexpr (std::is_same_v || std::is_same_v) { + if constexpr (IsMatrix) { p->flags = 0; - if constexpr (std::is_same_v) { - p->setDensityMatrixNodeFlag(generateDensityMatrix); - } } auto e = EdgeType::normalize(p, edges, memoryManager, cn); - if constexpr (std::is_same_v || std::is_same_v) { + if constexpr (IsMatrix) { if (!e.isTerminal()) { const auto& es = e.p->e; // Check if node resembles the identity. If so, skip it. @@ -673,8 +634,6 @@ class Package { char measureOneCollapsing(vEdge& rootEdge, Qubit index, std::mt19937_64& mt, fp epsilon = 0.001); - char measureOneCollapsing(dEdge& e, Qubit index, std::mt19937_64& mt); - /** * @brief Performs a specific measurement on the given state vector decision * diagram. Collapses the state according to the measurement result. @@ -695,8 +654,6 @@ class Package { config_.ctVecAddNumBucket}; ComputeTable matrixAdd{ config_.ctMatAddNumBucket}; - ComputeTable densityAdd{ - config_.ctDmAddNumBucket}; /** * @brief Get the compute table for addition operations. @@ -706,12 +663,10 @@ class Package { * type. */ template [[nodiscard]] auto& getAddComputeTable() { - if constexpr (std::is_same_v) { + if constexpr (IsVector) { return vectorAdd; - } else if constexpr (std::is_same_v) { + } else if constexpr (IsMatrix) { return matrixAdd; - } else if constexpr (std::is_same_v) { - return densityAdd; } } @@ -728,9 +683,9 @@ class Package { * type. */ template [[nodiscard]] auto& getAddMagnitudesComputeTable() { - if constexpr (std::is_same_v) { + if constexpr (IsVector) { return vectorAddMagnitudes; - } else if constexpr (std::is_same_v) { + } else if constexpr (IsMatrix) { return matrixAddMagnitudes; } } @@ -803,8 +758,7 @@ class Package { std::array, n> edge{}; for (std::size_t i = 0U; i < n; i++) { CachedEdge e1{}; - if constexpr (std::is_same_v || - std::is_same_v) { + if constexpr (IsMatrix) { if (x.isIdentity() || x.p->v < var) { // [ 0 | 1 ] [ x | 0 ] // --------- = --------- @@ -827,8 +781,7 @@ class Package { } } CachedEdge e2{}; - if constexpr (std::is_same_v || - std::is_same_v) { + if constexpr (IsMatrix) { if (y.isIdentity() || y.p->v < var) { // [ 0 | 1 ] [ y | 0 ] // --------- = --------- @@ -850,16 +803,7 @@ class Package { e2.w = y.w * ySuccessor.w; } } - - if constexpr (std::is_same_v) { - dNode::applyDmChangesToNode(e1.p); - dNode::applyDmChangesToNode(e2.p); - edge[i] = add2(e1, e2, var - 1); - dNode::revertDmChangesToNode(e2.p); - dNode::revertDmChangesToNode(e1.p); - } else { - edge[i] = add2(e1, e2, var - 1); - } + edge[i] = add2(e1, e2, var - 1); } auto r = makeDDNode(var, edge); computeTable.insert(x, y, r); @@ -906,8 +850,7 @@ class Package { std::array, n> edge{}; for (std::size_t i = 0U; i < n; i++) { CachedEdge e1{}; - if constexpr (std::is_same_v || - std::is_same_v) { + if constexpr (IsMatrix) { if (x.isIdentity() || x.p->v < var) { if (i == 0 || i == 3) { e1 = x; @@ -927,8 +870,7 @@ class Package { } } CachedEdge e2{}; - if constexpr (std::is_same_v || - std::is_same_v) { + if constexpr (IsMatrix) { if (y.isIdentity() || y.p->v < var) { if (i == 0 || i == 3) { e2 = y; @@ -1003,8 +945,6 @@ class Package { config_.ctMatVecMultNumBucket}; ComputeTable matrixMatrixMultiplication{ config_.ctMatMatMultNumBucket}; - ComputeTable densityDensityMultiplication{ - config_.ctDmDmMultNumBucket}; /** * @brief Get the compute table for multiplication operations. @@ -1019,8 +959,6 @@ class Package { return matrixVectorMultiplication; } else if constexpr (std::is_same_v) { return matrixMatrixMultiplication; - } else if constexpr (std::is_same_v) { - return densityDensityMultiplication; } } @@ -1053,8 +991,6 @@ class Package { MatrixDD applyOperation(const MatrixDD& operation, const MatrixDD& e, bool applyFromLeft = true); - dEdge applyOperationToDensity(dEdge& e, const mEdge& operation); - /** * @brief Multiplies two decision diagrams. * @@ -1062,8 +998,6 @@ class Package { * @tparam RightOperandNode The type of the right operand node. * @param x The left operand decision diagram. * @param y The right operand decision diagram. - * @param generateDensityMatrix Flag to indicate if a density matrix node - * should be generated. * @return The resulting decision diagram after multiplication. * * @details This function performs the multiplication of two decision diagrams @@ -1076,44 +1010,20 @@ class Package { * canonical form. */ template - Edge - multiply(const Edge& x, const Edge& y, - [[maybe_unused]] const bool generateDensityMatrix = false) { - using LEdge = Edge; - using REdge = Edge; - static_assert(std::disjunction_v, - std::is_same>, - "Left operand must be a matrix or density matrix"); - static_assert(std::disjunction_v, - std::is_same, - std::is_same>, - "Right operand must be a vector, matrix or density matrix"); + requires IsMatrix && + (IsVector || IsMatrix) + Edge multiply(const Edge& x, + const Edge& y) { Qubit var{}; - if constexpr (std::is_same_v) { - 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 cn.lookup(e); - } else { - if (!x.isTerminal()) { - var = x.p->v; - } - if (!y.isTerminal() && y.p->v > var) { - var = y.p->v; - } - const auto e = multiply2(x, y, var); - return cn.lookup(e); + if (!x.isTerminal()) { + var = x.p->v; } + if (!y.isTerminal() && y.p->v > var) { + var = y.p->v; + } + const auto e = multiply2(x, y, var); + return cn.lookup(e); } private: @@ -1128,15 +1038,12 @@ class Package { * @param x The left operand decision diagram. * @param y The right operand decision diagram. * @param var The variable associated with the current level of recursion. - * @param generateDensityMatrix Flag to indicate if a density matrix node - * should be generated. * @return The resulting DD after multiplication. */ template - CachedEdge - multiply2(const Edge& x, const Edge& y, - const Qubit var, - [[maybe_unused]] const bool generateDensityMatrix = false) { + CachedEdge multiply2(const Edge& x, + const Edge& y, + const Qubit var) { using LEdge = Edge; using REdge = Edge; using ResultEdge = CachedEdge; @@ -1149,39 +1056,17 @@ class Package { const auto yWeight = static_cast(y.w); const auto rWeight = xWeight * yWeight; if (x.isIdentity()) { - if constexpr (!std::is_same_v) { - return {y.p, rWeight}; - } else { - if (y.isIdentity() || - (dNode::isDensityMatrixTempFlagSet(y.p->flags) && - generateDensityMatrix) || - (!dNode::isDensityMatrixTempFlagSet(y.p->flags) && - !generateDensityMatrix)) { - return {y.p, rWeight}; - } - } + return {y.p, rWeight}; } - if constexpr (std::is_same_v || - std::is_same_v) { + if constexpr (std::is_same_v) { if (y.isIdentity()) { - if constexpr (!std::is_same_v) { - return {x.p, rWeight}; - } else { - if (x.isIdentity() || - (dNode::isDensityMatrixTempFlagSet(x.p->flags) && - generateDensityMatrix) || - (!dNode::isDensityMatrixTempFlagSet(x.p->flags) && - !generateDensityMatrix)) { - return {x.p, rWeight}; - } - } + return {x.p, rWeight}; } } auto& computeTable = getMultiplicationComputeTable(); - if (const auto* r = computeTable.lookup(x.p, y.p, generateDensityMatrix); - r != nullptr) { + if (const auto* r = computeTable.lookup(x.p, y.p); r != nullptr) { return {r->p, r->w * rWeight}; } @@ -1221,54 +1106,18 @@ class Package { } const auto v = static_cast(var - 1); - if constexpr (std::is_same_v) { - 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); - } else { - auto m = multiply2(e1, e2, v); + auto m = multiply2(e1, e2, v); - if (k == 0 || edge[idx].w.exactlyZero()) { - edge[idx] = m; - } else if (!m.w.exactlyZero()) { - edge[idx] = add2(edge[idx], m, v); - } + if (k == 0 || edge[idx].w.exactlyZero()) { + edge[idx] = m; + } else if (!m.w.exactlyZero()) { + edge[idx] = add2(edge[idx], m, v); } } } } - auto e = makeDDNode(var, edge, generateDensityMatrix); + auto e = makeDDNode(var, edge); computeTable.insert(x.p, y.p, e); e.w = e.w * rWeight; @@ -1389,7 +1238,7 @@ class Package { * type. */ template [[nodiscard]] auto& getKroneckerComputeTable() { - if constexpr (std::is_same_v) { + if constexpr (IsVector) { return vectorKronecker; } else { return matrixKronecker; @@ -1407,17 +1256,10 @@ class Package { * decision diagram. * @return The resulting decision diagram after computing the Kronecker * product. - * @throws std::invalid_argument if the node type is `dNode` (density - * matrices). */ template Edge kronecker(const Edge& x, const Edge& y, const std::size_t yNumQubits, const bool incIdx = true) { - if constexpr (std::is_same_v) { - throw std::invalid_argument( - "Kronecker is currently not supported for density matrices"); - } - const auto e = kronecker2(x, y, yNumQubits, incIdx); return cn.lookup(e); } @@ -1462,7 +1304,7 @@ class Package { return {x.p, rWeight}; } - if constexpr (std::is_same_v || std::is_same_v) { + if constexpr (IsMatrix) { if (x.isIdentity()) { return {y.p, rWeight}; } @@ -1491,8 +1333,7 @@ class Package { Qubit idx = x.p->v; if (incIdx) { // use the given number of qubits if y is an identity - if constexpr (std::is_same_v || - std::is_same_v) { + if constexpr (IsMatrix) { if (y.isIdentity()) { idx += static_cast(yNumQubits); } else { @@ -1502,7 +1343,7 @@ class Package { idx += static_cast(y.p->v + 1U); } } - auto e = makeDDNode(idx, edge, true); + auto e = makeDDNode(idx, edge); computeTable.insert(x.p, y.p, {e.p, e.w}); return {e.p, rWeight}; } @@ -1511,8 +1352,6 @@ class Package { /// (Partial) trace /// public: - UnaryComputeTable densityTrace{ - config_.ctDmTraceNumBucket}; UnaryComputeTable matrixTrace{ config_.ctMatTraceNumBucket}; @@ -1523,13 +1362,7 @@ class Package { * @return A reference to the appropriate compute table for the given node * type. */ - template [[nodiscard]] auto& getTraceComputeTable() { - if constexpr (std::is_same_v) { - return matrixTrace; - } else { - return densityTrace; - } - } + [[nodiscard]] auto& getTraceComputeTable() { return matrixTrace; } /** * @brief Computes the partial trace of a matrix decision diagram. @@ -1541,21 +1374,13 @@ class Package { mEdge partialTrace(const mEdge& a, const std::vector& eliminate); /** - * @brief Computes the trace of a decision diagram. + * @brief Computes the trace of a matrix decision diagram. * - * @tparam Node The type of the node. * @param a The decision diagram. * @param numQubits The number of qubits in the decision diagram. * @return The trace of the decision diagram as a complex value. */ - template - ComplexValue trace(const Edge& 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; - } + ComplexValue trace(const mEdge& a, std::size_t numQubits); /** * @brief Checks if a given matrix is close to the identity matrix. @@ -1589,80 +1414,9 @@ class Package { * For matrices, normalization is continuously applied, dividing by two at * each level marked for elimination, thereby ensuring that the result is * mapped to the interval [0,1] (as opposed to the interval [0,2^N]). - * - * For density matrices, such normalization is not applied as the trace of - * density matrices is always 1 by definition. */ - template - CachedEdge trace(const Edge& a, - const std::vector& eliminate, std::size_t level, - std::size_t alreadyEliminated = 0) { - const auto aWeight = static_cast(a.w); - if (aWeight.approximatelyZero()) { - return CachedEdge::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 CachedEdge{a.p, aWeight}; - } - - const auto v = a.p->v; - if (eliminate[v]) { - // Lookup nodes marked for elimination in the compute table if all - // lower-level qubits are eliminated as well: if the trace has already - // been computed, return the result - 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 = getTraceComputeTable().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); - - // The resulting weight is continuously normalized to the range [0,1] for - // matrix nodes - if constexpr (std::is_same_v) { - r.w = r.w / 2.0; - } - - // Insert result into compute table if all lower-level qubits are - // eliminated as well - if (eliminateAll) { - getTraceComputeTable().insert(a.p, r); - } - r.w = r.w * aWeight; - return r; - } - - std::array, NEDGE> edge{}; - std::transform(a.p->e.cbegin(), a.p->e.cend(), edge.begin(), - [this, &eliminate, &alreadyEliminated, - &level](const Edge& e) -> CachedEdge { - 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; - } + mCachedEdge trace(const mEdge& a, const std::vector& eliminate, + std::size_t level, std::size_t alreadyEliminated = 0); /** * @brief Recursively checks if a given matrix is close to the identity @@ -1692,13 +1446,6 @@ class Package { mEdge createInitialMatrix(const std::vector& ancillary); - /// - /// Noise Operations - /// - StochasticNoiseOperationTable stochasticNoiseOperationCache{ - nqubits, config_.stochasticCacheOps}; - DensityNoiseTable densityNoise{config_.ctDmNoiseNumBucket}; - /// /// Ancillary and garbage reduction /// @@ -1921,7 +1668,7 @@ class Package { } else { std::string version; std::getline(is, version); - if (std::stoi(version) != SERIALIZATION_VERSION) { + if (std::cmp_not_equal(std::stoi(version), SERIALIZATION_VERSION)) { throw std::runtime_error( "Wrong Version of serialization file version. version of file: " + version + diff --git a/include/mqt-core/dd/Package_fwd.hpp b/include/mqt-core/dd/Package_fwd.hpp index ec58b3023f..2522437eef 100644 --- a/include/mqt-core/dd/Package_fwd.hpp +++ b/include/mqt-core/dd/Package_fwd.hpp @@ -18,5 +18,4 @@ class Package; using VectorDD = Edge; using MatrixDD = Edge; -using DensityMatrixDD = Edge; } // namespace dd diff --git a/include/mqt-core/dd/RealNumber.hpp b/include/mqt-core/dd/RealNumber.hpp index 6456fe23a0..705d574895 100644 --- a/include/mqt-core/dd/RealNumber.hpp +++ b/include/mqt-core/dd/RealNumber.hpp @@ -214,7 +214,7 @@ struct RealNumber final : LLBase { /// numerical tolerance to be used for floating point values // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables) - static inline fp eps = std::numeric_limits::epsilon() * 1024; + static inline fp eps = std::numeric_limits::epsilon() * 1024; }; static_assert(sizeof(RealNumber) == 16); diff --git a/include/mqt-core/dd/StateGeneration.hpp b/include/mqt-core/dd/StateGeneration.hpp index eac3e7114c..9a4c15ad85 100644 --- a/include/mqt-core/dd/StateGeneration.hpp +++ b/include/mqt-core/dd/StateGeneration.hpp @@ -10,6 +10,7 @@ #pragma once +#include "dd/DDDefinitions.hpp" #include "dd/Node.hpp" #include "dd/Package.hpp" @@ -87,7 +88,7 @@ VectorDD makeWState(std::size_t n, Package& dd); VectorDD makeStateFromVector(const CVec& vec, Package& dd); /// @brief The strategy to wire two layers. -enum GenerationWireStrategy : std::uint8_t { +enum class GenerationWireStrategy : std::uint8_t { ROUNDROBIN, // Choose nodes in the next layer in a round-robin fashion. RANDOM // Randomly choose nodes in the next layer. }; diff --git a/include/mqt-core/dd/StochasticNoiseOperationTable.hpp b/include/mqt-core/dd/StochasticNoiseOperationTable.hpp deleted file mode 100644 index 960829e4d8..0000000000 --- a/include/mqt-core/dd/StochasticNoiseOperationTable.hpp +++ /dev/null @@ -1,88 +0,0 @@ -/* - * 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 - */ - -#pragma once - -#include "dd/statistics/TableStatistics.hpp" -#include "ir/Definitions.hpp" -#include "ir/operations/OpType.hpp" - -#include -#include -#include -#include - -namespace dd { -template class StochasticNoiseOperationTable { -public: - explicit StochasticNoiseOperationTable( - const std::size_t nv, - const size_t numberOfStochasticOperations = qc::OpType::OpTypeEnd) - : nvars(nv), numberOfStochasticOperations_(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(numberOfStochasticOperations_)); - } - - void insert(std::uint8_t kind, qc::Qubit target, const Edge& r) { - assert(kind < - numberOfStochasticOperations_); // There are new operations in - // OpType. Increase the value of - // numberOfOperations accordingly - table.at(target).at(kind) = r; - stats.trackInsert(); - } - - Edge* lookup(std::uint8_t kind, qc::Qubit target) { - assert(kind < - numberOfStochasticOperations_); // There are new operations in - // OpType. Increase the value of - // numberOfOperations accordingly - ++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; - size_t numberOfStochasticOperations_; - std::vector> table; - TableStatistics stats{}; -}; -} // namespace dd diff --git a/include/mqt-core/dd/UniqueTable.hpp b/include/mqt-core/dd/UniqueTable.hpp index 613027a9ff..44bdfb3944 100644 --- a/include/mqt-core/dd/UniqueTable.hpp +++ b/include/mqt-core/dd/UniqueTable.hpp @@ -26,6 +26,7 @@ #include #include #include +#include #include #include @@ -89,11 +90,7 @@ class UniqueTable { template [[nodiscard]] static bool nodesAreEqual(const Node& p, const Node& q) { - if constexpr (std::is_same_v) { - return (p.e == q.e && (p.flags == q.flags)); - } else { - return p.e == q.e; - } + return p.e == q.e; } // Lookup a node in the unique table for the appropriate variable and insert @@ -158,10 +155,8 @@ class UniqueTable { static_assert(std::is_base_of_v, "Node must be derived from NodeBase"); auto q = cfg.nVars - 1U; - for (auto it = tables.rbegin(); it != tables.rend(); ++it) { - auto& table = *it; - std::cout << "\tq" << q << ":" - << "\n"; + for (const auto& table : std::ranges::reverse_view(tables)) { + std::cout << "\tq" << q << ":" << "\n"; for (std::size_t key = 0; key < table.size(); ++key) { auto* p = static_cast(table[key]); if (p != nullptr) { diff --git a/include/mqt-core/dd/statistics/PackageStatistics.hpp b/include/mqt-core/dd/statistics/PackageStatistics.hpp index a6d96758ee..c410d7c3af 100644 --- a/include/mqt-core/dd/statistics/PackageStatistics.hpp +++ b/include/mqt-core/dd/statistics/PackageStatistics.hpp @@ -39,14 +39,14 @@ namespace dd { */ [[nodiscard]] double computePeakMemoryMiB(const Package& package); -[[nodiscard]] nlohmann::basic_json<> +[[nodiscard]] nlohmann::json getStatistics(Package& package, bool includeIndividualTables = false); /** * @brief Get some key statistics about data structures used by the DD package * @return A JSON representation of the statistics */ -[[nodiscard]] nlohmann::basic_json<> getDataStructureStatistics(); +[[nodiscard]] nlohmann::json getDataStructureStatistics(); [[nodiscard]] std::string getStatisticsString(Package& package); diff --git a/include/mqt-core/dd/statistics/TableStatistics.hpp b/include/mqt-core/dd/statistics/TableStatistics.hpp index 5762acc716..6914e557af 100644 --- a/include/mqt-core/dd/statistics/TableStatistics.hpp +++ b/include/mqt-core/dd/statistics/TableStatistics.hpp @@ -18,7 +18,7 @@ namespace dd { /// A utility class for storing statistics of a table -struct TableStatistics : public Statistics { +struct TableStatistics : Statistics { /// The size of a single entry std::size_t entrySize = 0U; /// The number of buckets in the table diff --git a/include/mqt-core/dd/statistics/UniqueTableStatistics.hpp b/include/mqt-core/dd/statistics/UniqueTableStatistics.hpp index c192ac1bc2..2b05f4f5e4 100644 --- a/include/mqt-core/dd/statistics/UniqueTableStatistics.hpp +++ b/include/mqt-core/dd/statistics/UniqueTableStatistics.hpp @@ -17,7 +17,7 @@ namespace dd { /// \brief A class for storing statistics of a unique table -struct UniqueTableStatistics : public TableStatistics { +struct UniqueTableStatistics : TableStatistics { /// The number of garbage collection runs std::size_t gcRuns = 0U; diff --git a/include/mqt-core/ir/operations/OpType.inc b/include/mqt-core/ir/operations/OpType.inc index ca3a225749..fdea9c2d46 100644 --- a/include/mqt-core/ir/operations/OpType.inc +++ b/include/mqt-core/ir/operations/OpType.inc @@ -57,12 +57,6 @@ HANDLE_OP_TYPE(31, Reset, OpTypeNone, "reset") // Control flow operations HANDLE_OP_TYPE(32, IfElse, OpTypeNone, "if_else") -// Noise operations -HANDLE_OP_TYPE(33, ATrue, 0, "a_true") -HANDLE_OP_TYPE(34, AFalse, 0, "a_false") -HANDLE_OP_TYPE(35, MultiATrue, 0, "multi_a_true") -HANDLE_OP_TYPE(36, MultiAFalse, 0, "multi_a_false") - // Neutral atom shuttling operations HANDLE_OP_TYPE(37, Move, 0, "move") HANDLE_OP_TYPE(42, Bridge, 0, "bridge") diff --git a/src/dd/CachedEdge.cpp b/src/dd/CachedEdge.cpp index 58c4e25712..fa4bd223f6 100644 --- a/src/dd/CachedEdge.cpp +++ b/src/dd/CachedEdge.cpp @@ -33,11 +33,12 @@ namespace dd { ///----------------------------------------------------------------------------- template -template > auto CachedEdge::normalize(Node* p, const std::array& e, MemoryManager& mm, ComplexNumbers& cn) - -> CachedEdge { + -> CachedEdge + requires IsVector +{ assert(p != nullptr && "Node pointer passed to normalize is null."); const auto zero = std::array{e[0].w.approximatelyZero(), e[1].w.approximatelyZero()}; @@ -93,11 +94,12 @@ auto CachedEdge::normalize(Node* p, ///----------------------------------------------------------------------------- template -template > auto CachedEdge::normalize(Node* p, const std::array& e, MemoryManager& mm, ComplexNumbers& cn) - -> CachedEdge { + -> CachedEdge + requires IsMatrix +{ assert(p != nullptr && "Node pointer passed to normalize is null."); const auto zero = std::array{e[0].w.approximatelyZero(), e[1].w.approximatelyZero(), @@ -157,22 +159,6 @@ auto CachedEdge::normalize(Node* p, template struct CachedEdge; template struct CachedEdge; -template struct CachedEdge; - -template auto -CachedEdge::normalize(vNode* p, const std::array& e, - MemoryManager& mm, ComplexNumbers& cn) - -> CachedEdge; - -template auto -CachedEdge::normalize(mNode* p, const std::array& e, - MemoryManager& mm, ComplexNumbers& cn) - -> CachedEdge; - -template auto -CachedEdge::normalize(dNode* p, const std::array& e, - MemoryManager& mm, ComplexNumbers& cn) - -> CachedEdge; } // namespace dd @@ -187,5 +173,4 @@ auto hash>::operator()( template struct hash>; template struct hash>; -template struct hash>; } // namespace std diff --git a/src/dd/Edge.cpp b/src/dd/Edge.cpp index 5e01d28a63..15cf8159cb 100644 --- a/src/dd/Edge.cpp +++ b/src/dd/Edge.cpp @@ -48,10 +48,6 @@ auto Edge::getValueByPath(const std::size_t numQubits, } auto r = *this; - if constexpr (std::is_same_v) { - Edge::applyDmChangesToEdge(r); - } - auto level = numQubits; while (level > 0U) { const auto tmp = static_cast(decisions.at(level - 1U) - '0'); @@ -67,14 +63,7 @@ auto Edge::getValueByPath(const std::size_t numQubits, // node is at the expected level assert(tmp <= r.p->e.size()); - if constexpr (std::is_same_v) { - auto e = r; - Edge::applyDmChangesToEdge(r.p->e[tmp]); - r = r.p->e[tmp]; - Edge::revertDmChangesToEdge(e); - } else { - r = r.p->e[tmp]; - } + r = r.p->e[tmp]; c *= static_cast>(r.w); --level; } @@ -131,16 +120,17 @@ template void Edge::unmark() const noexcept { ///----------------------------------------------------------------------------- template -template > auto Edge::normalize(Node* p, const std::array& e, - MemoryManager& mm, ComplexNumbers& cn) -> Edge { + MemoryManager& mm, ComplexNumbers& cn) -> Edge + requires IsVector +{ assert(p != nullptr && "Node pointer passed to normalize is null."); const auto zero = std::array{e[0].w.exactlyZero(), e[1].w.exactlyZero()}; if (zero[0]) { if (zero[1]) { mm.returnEntry(*p); - return Edge::zero(); + return vEdge::zero(); } p->e = e; vEdge r{p, e[1].w}; @@ -190,16 +180,16 @@ auto Edge::normalize(Node* p, const std::array& e, if (min.w.exactlyZero()) { assert(p->e[argMax].w.exactlyOne() && "Edge weight should be one when minWeight is zero."); - min.p = Node::getTerminal(); + min.p = vNode::getTerminal(); } return r; } template -template > -auto Edge::getValueByIndex(const std::size_t i) const - -> std::complex { +auto Edge::getValueByIndex(const std::size_t i) const -> std::complex + requires IsVector +{ auto bitwidth = static_cast(std::log2(i + 1U)); if (!isTerminal()) { @@ -217,8 +207,9 @@ auto Edge::getValueByIndex(const std::size_t i) const } template -template > -auto Edge::getVector(const fp threshold) const -> CVec { +auto Edge::getVector(const fp threshold) const -> CVec + requires IsVector +{ if (isTerminal()) { return {static_cast>(w)}; } @@ -233,8 +224,9 @@ auto Edge::getVector(const fp threshold) const -> CVec { } template -template > -auto Edge::getSparseVector(const fp threshold) const -> SparseCVec { +auto Edge::getSparseVector(const fp threshold) const -> SparseCVec + requires IsVector +{ if (isTerminal()) { return {{0, static_cast>(w)}}; } @@ -248,8 +240,9 @@ auto Edge::getSparseVector(const fp threshold) const -> SparseCVec { } template -template > -auto Edge::printVector() const -> void { +auto Edge::printVector() const -> void + requires IsVector +{ constexpr auto precision = 3; const auto oldPrecision = std::cout.precision(); std::cout << std::setprecision(precision); @@ -272,8 +265,9 @@ auto Edge::printVector() const -> void { } template -template > -auto Edge::addToVector(CVec& amplitudes) const -> void { +auto Edge::addToVector(CVec& amplitudes) const -> void + requires IsVector +{ if (isTerminal()) { amplitudes[0] += static_cast>(w); return; @@ -286,10 +280,11 @@ auto Edge::addToVector(CVec& amplitudes) const -> void { } template -template > void Edge::traverseVector(const std::complex& amp, const std::size_t i, AmplitudeFunc f, - const fp threshold) const { + const fp threshold) const + requires IsVector +{ // calculate new accumulated amplitude const auto c = amp * static_cast>(w); @@ -315,9 +310,10 @@ void Edge::traverseVector(const std::complex& amp, /// \n Methods for matrix DDs \n ///----------------------------------------------------------------------------- template -template > auto Edge::normalize(Node* p, const std::array& e, - MemoryManager& mm, ComplexNumbers& cn) -> Edge { + MemoryManager& mm, ComplexNumbers& cn) -> Edge + requires IsMatrix +{ 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()}; @@ -374,10 +370,11 @@ auto Edge::normalize(Node* p, const std::array& e, } template -template > auto Edge::getValueByIndex(const std::size_t numQubits, const std::size_t i, const std::size_t j) const - -> std::complex { + -> std::complex + requires IsMatrix +{ if (isTerminal()) { return static_cast>(w); } @@ -402,60 +399,47 @@ auto Edge::getValueByIndex(const std::size_t numQubits, } template -template > auto Edge::getMatrix(const std::size_t numQubits, - const fp threshold) const -> CMat { + const fp threshold) const -> CMat + requires IsMatrix +{ if (numQubits == 0U) { return CMat{1, {static_cast>(w)}}; } - auto r = *this; - if constexpr (std::is_same_v) { - Edge::applyDmChangesToEdge(r); - } const std::size_t dim = 1ULL << numQubits; auto mat = CMat(dim, CVec(dim, 0.)); - r.traverseMatrix( + traverseMatrix( 1, 0ULL, 0ULL, [&mat](const std::size_t i, const std::size_t j, const std::complex& c) { mat.at(i).at(j) = c; }, numQubits, threshold); - if constexpr (std::is_same_v) { - Edge::revertDmChangesToEdge(r); - } return mat; } template -template > auto Edge::getSparseMatrix(const std::size_t numQubits, - const fp threshold) const -> SparseCMat { + const fp threshold) const -> SparseCMat + requires IsMatrix +{ if (numQubits == 0U) { return {{{0U, 0U}, static_cast>(w)}}; } - auto r = *this; - if constexpr (std::is_same_v) { - Edge::applyDmChangesToEdge(r); - } - auto mat = SparseCMat{}; - r.traverseMatrix( + traverseMatrix( 1, 0ULL, 0ULL, [&mat](const std::size_t i, const std::size_t j, const std::complex& c) { mat[{i, j}] = c; }, numQubits, threshold); - if constexpr (std::is_same_v) { - Edge::revertDmChangesToEdge(r); - } - return mat; } template -template > -auto Edge::printMatrix(const std::size_t numQubits) const -> void { +auto Edge::printMatrix(const std::size_t numQubits) const -> void + requires IsMatrix +{ constexpr auto precision = 3; const auto oldPrecision = std::cout.precision(); std::cout << std::setprecision(precision); @@ -464,14 +448,8 @@ auto Edge::printMatrix(const std::size_t numQubits) const -> void { std::cout << static_cast>(w) << "\n"; return; } - - auto r = *this; - if constexpr (std::is_same_v) { - Edge::alignDensityEdge(r); - } - // total number of qubits should not be lower than the highest qubit index - assert(r.isTerminal() || numQubits >= r.p->v); + assert(isTerminal() || numQubits >= p->v); const std::size_t element = 1ULL << numQubits; for (auto i = 0ULL; i < element; ++i) { for (auto j = 0ULL; j < element; ++j) { @@ -485,11 +463,12 @@ auto Edge::printMatrix(const std::size_t numQubits) const -> void { } template -template > -auto Edge::traverseMatrix(const std::complex& amp, +void Edge::traverseMatrix(const std::complex& amp, const std::size_t i, const std::size_t j, MatrixEntryFunc f, const std::size_t level, - const fp threshold) const -> void { + const fp threshold) const + requires IsMatrix +{ // calculate new accumulated amplitude const auto c = amp * static_cast>(w); @@ -516,161 +495,17 @@ auto Edge::traverseMatrix(const std::complex& amp, std::size_t k = 0U; for (const auto& [a, b] : coords) { if (auto& e = p->e[k++]; !e.w.exactlyZero()) { - if constexpr (std::is_same_v) { - Edge::applyDmChangesToEdge(e); - } e.traverseMatrix(c, a, b, f, nextLevel, threshold); - if constexpr (std::is_same_v) { - Edge::revertDmChangesToEdge(e); - } } } } -///----------------------------------------------------------------------------- -/// \n Methods for density matrix DDs \n -///----------------------------------------------------------------------------- - -template -template > -auto Edge::getSparseProbabilityVector(const std::size_t numQubits, - const fp threshold) const - -> SparsePVec { - if (numQubits == 0U) { - return {{0, static_cast>(w).real()}}; - } - - auto e = *this; - Edge::alignDensityEdge(e); - - auto probabilities = SparsePVec{}; - e.traverseDiagonal( - 1, 0, - [&probabilities](const std::size_t i, const fp& prob) { - probabilities[i] = prob; - }, - numQubits, threshold); - return probabilities; -} - -template -template > -auto Edge::getSparseProbabilityVectorStrKeys(const std::size_t numQubits, - const fp threshold) const - -> SparsePVecStrKeys { - if (numQubits == 0U) { - return {{"0", static_cast>(w).real()}}; - } - - auto e = *this; - Edge::alignDensityEdge(e); - const auto nqubits = static_cast(e.p->v) + 1U; - - auto probabilities = SparsePVecStrKeys{}; - e.traverseDiagonal( - 1, 0, - [&probabilities, &nqubits](const std::size_t i, const fp& prob) { - probabilities[intToBinaryString(i, nqubits)] = prob; - }, - numQubits, threshold); - return probabilities; -} - -template -template > -auto Edge::traverseDiagonal(const fp& prob, const std::size_t i, - ProbabilityFunc f, const std::size_t level, - const fp threshold) const -> void { - // 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 Explicit instantiations \n ///----------------------------------------------------------------------------- template struct Edge; template struct Edge; -template struct Edge; - -template auto Edge::normalize(vNode* p, const std::array& e, - MemoryManager& mm, ComplexNumbers& cn) - -> Edge; -template auto Edge::getValueByIndex(std::size_t i) const - -> std::complex; -template auto Edge::getVector(fp threshold) const -> CVec; -template auto Edge::getSparseVector(fp threshold) const -> SparseCVec; -template auto Edge::printVector() const -> void; -template auto Edge::addToVector(CVec& amplitudes) const -> void; -template auto Edge::traverseVector(const std::complex& amp, - std::size_t i, AmplitudeFunc f, - fp threshold) const -> void; - -template auto Edge::normalize(mNode* p, const std::array& e, - MemoryManager& mm, ComplexNumbers& cn) - -> Edge; -template auto Edge::getValueByIndex(std::size_t numQubits, std::size_t i, - std::size_t j) const - -> std::complex; -template auto Edge::getMatrix(std::size_t numQubits, fp threshold) const - -> CMat; -template auto Edge::getSparseMatrix(std::size_t numQubits, - fp threshold) const -> SparseCMat; -template auto Edge::printMatrix(std::size_t numQubits) const -> void; -template auto Edge::traverseMatrix(const std::complex& amp, - std::size_t i, std::size_t j, - MatrixEntryFunc f, std::size_t level, - fp threshold) const -> void; - -template auto Edge::normalize(dNode* p, const std::array& e, - MemoryManager& mm, ComplexNumbers& cn) - -> Edge; -template auto Edge::getMatrix(std::size_t numQubits, fp threshold) const - -> CMat; -template auto Edge::getSparseMatrix(std::size_t numQubits, - fp threshold) const -> SparseCMat; -template auto Edge::printMatrix(std::size_t numQubits) const -> void; -template auto Edge::getSparseProbabilityVector(std::size_t numQubits, - fp threshold) const - -> SparsePVec; -template auto Edge::getSparseProbabilityVectorStrKeys( - std::size_t numQubits, fp threshold) const -> SparsePVecStrKeys; -template auto Edge::getValueByIndex(std::size_t numQubits, - std::size_t i, - std::size_t j) const - -> std::complex; -template auto Edge::traverseMatrix( - const std::complex& amp, std::size_t i, std::size_t j, - MatrixEntryFunc f, std::size_t level, fp threshold) const -> void; -template auto Edge::traverseDiagonal(const fp& prob, std::size_t i, - ProbabilityFunc f, - std::size_t level, - fp threshold) const -> void; } // namespace dd @@ -684,19 +519,9 @@ auto hash>::operator()(const dd::Edge& e) const noexcept -> std::size_t { const auto h1 = dd::murmur64(reinterpret_cast(e.p)); const auto h2 = std::hash{}(e.w); - auto h3 = qc::combineHash(h1, h2); - if constexpr (std::is_same_v) { - if (e.isTerminal()) { - return h3; - } - assert(dd::dNode::isDensityMatrixTempFlagSet(e.p) == false); - const auto h4 = dd::dNode::getDensityMatrixTempFlags(e.p->flags); - h3 = qc::combineHash(h3, h4); - } - return h3; + return qc::combineHash(h1, h2); } template struct hash>; template struct hash>; -template struct hash>; } // namespace std diff --git a/src/dd/Node.cpp b/src/dd/Node.cpp deleted file mode 100644 index 7648abb702..0000000000 --- a/src/dd/Node.cpp +++ /dev/null @@ -1,106 +0,0 @@ -/* - * 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 "dd/Node.hpp" - -#include "dd/ComplexNumbers.hpp" - -#include -#include -#include - -namespace dd { - -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 = 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 = 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 = ComplexNumbers::conj(p->e[2].w); - return; - } - for (auto& edge : p->e) { - edge.w = 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(); - } -} - -} // namespace dd diff --git a/src/dd/NoiseFunctionality.cpp b/src/dd/NoiseFunctionality.cpp deleted file mode 100644 index 8a96a460fa..0000000000 --- a/src/dd/NoiseFunctionality.cpp +++ /dev/null @@ -1,469 +0,0 @@ -/* - * 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 "dd/NoiseFunctionality.hpp" - -#include "dd/ComplexNumbers.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" -#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::AmplitudeDamping); - break; - case 'P': - noiseOperationVector.emplace_back(dd::PhaseFlip); - break; - case 'D': - noiseOperationVector.emplace_back(dd::Depolarization); - break; - case 'I': - noiseOperationVector.emplace_back(dd::Identity); - break; - default: - throw std::runtime_error("Unknown noise operation '" + cNoiseEffects + - "'\n"); - } - } - return noiseOperationVector; -} -} // namespace - -namespace dd { -StochasticNoiseFunctionality::StochasticNoiseFunctionality( - 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(Package::makeIdent()) { - sanityCheckOfNoiseProbabilities(gateNoiseProbability, amplitudeDampingProb, - multiQubitGateFactor); - package->incRef(identityDD); -} - -double StochasticNoiseFunctionality::getNoiseProbability( - const bool multiQubitNoiseFlag) const { - return multiQubitNoiseFlag ? noiseProbabilityMulti : noiseProbability; -} - -qc::OpType StochasticNoiseFunctionality::getAmplitudeDampingOperationType( - const bool multiQubitNoiseFlag, const bool amplitudeDampingFlag) { - if (amplitudeDampingFlag) { - return multiQubitNoiseFlag ? qc::MultiATrue : qc::ATrue; - } - return multiQubitNoiseFlag ? qc::MultiAFalse : qc::AFalse; -} - -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, mEdge operation, 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 (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 = Complex::one(); - - package->incRef(tmp); - package->decRef(state); - state = tmp; - - // I only need to apply the operations once - operation = identityDD; - } -} - -mEdge StochasticNoiseFunctionality::stackOperation( - const mEdge& operation, const qc::Qubit target, - const qc::OpType noiseOperation, const GateMatrix& matrix) const { - if (const auto* op = - package->stochasticNoiseOperationCache.lookup(noiseOperation, target); - op != nullptr) { - return package->multiply(*op, operation); - } - const auto gateDD = package->makeGateDD(matrix, target); - package->stochasticNoiseOperationCache.insert(noiseOperation, target, gateDD); - return package->multiply(gateDD, operation); -} - -mEdge StochasticNoiseFunctionality::generateNoiseOperation( - 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 (qc::I): { - continue; - } - case (qc::MultiATrue): - case (qc::ATrue): { - const GateMatrix amplitudeDampingMatrix = - getAmplitudeDampingOperationMatrix(multiQubitOperation, true); - operation = - stackOperation(operation, target, effect, amplitudeDampingMatrix); - break; - } - case (qc::MultiAFalse): - case (qc::AFalse): { - const GateMatrix amplitudeDampingMatrix = - getAmplitudeDampingOperationMatrix(multiQubitOperation, false); - operation = - stackOperation(operation, target, effect, amplitudeDampingMatrix); - break; - } - case (qc::X): - case (qc::Y): - case (qc::Z): { - operation = stackOperation(operation, target, effect, - opToSingleQubitGateMatrix(effect)); - break; - } - default: { - throw std::runtime_error("Unknown noise operation '" + - std::to_string(effect) + "'\n"); - } - } - } - return operation; -} - -qc::OpType 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 qc::I, also 25 % of the time when depolarization is - // applied nothing happens - return qc::I; - } - if (prob < (getNoiseProbability(multiQubitNoiseFlag) * 0.25)) { - // if 0 < prob < 0.25 (25 % of the time when applying depolarization) - // apply qc::X - return qc::X; - } - if (prob < (getNoiseProbability(multiQubitNoiseFlag) * 0.5)) { - // if 0.25 < prob < 0.5 (25 % of the time when applying depolarization) - // apply qc::Y - return qc::Y; - } - // if 0.5 < prob < 0.75 (25 % of the time when applying depolarization) - // apply qc::Z - return qc::Z; - } - case PhaseFlip: { - if (prob > getNoiseProbability(multiQubitNoiseFlag)) { - return qc::I; - } - return qc::Z; - } - case Identity: { - return qc::I; - } - default: - throw std::runtime_error(std::string{"Unknown noise effect '"} + - std::to_string(noiseOperation) + "'"); - } -} - -DeterministicNoiseFunctionality::DeterministicNoiseFunctionality( - Package& 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->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 Qubit level) { - - const auto originalWeight = static_cast(originalEdge.w); - if (originalEdge.isZeroTerminal() || level <= *usedQubits.begin()) { - return {originalEdge.p, originalWeight}; - } - - auto originalCopy = dEdge{originalEdge.p, 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 diff --git a/src/dd/Package.cpp b/src/dd/Package.cpp index b169f8c776..b9e01dfe3a 100644 --- a/src/dd/Package.cpp +++ b/src/dd/Package.cpp @@ -17,14 +17,12 @@ #include "dd/ComputeTable.hpp" #include "dd/DDDefinitions.hpp" #include "dd/DDpackageConfig.hpp" -#include "dd/DensityNoiseTable.hpp" #include "dd/Edge.hpp" #include "dd/GateMatrixDefinitions.hpp" #include "dd/MemoryManager.hpp" #include "dd/Node.hpp" #include "dd/RealNumber.hpp" #include "dd/RealNumberUniqueTable.hpp" -#include "dd/StochasticNoiseOperationTable.hpp" #include "dd/UnaryComputeTable.hpp" #include "dd/UniqueTable.hpp" #include "ir/Definitions.hpp" @@ -68,8 +66,6 @@ void Package::resize(const std::size_t nq) { nqubits = nq; vUniqueTable.resize(nqubits); mUniqueTable.resize(nqubits); - dUniqueTable.resize(nqubits); - stochasticNoiseOperationCache.resize(nqubits); } void Package::reset() { @@ -82,24 +78,21 @@ void Package::reset() { void Package::resetMemoryManagers(const bool resizeToTotal) { vMemoryManager.reset(resizeToTotal); mMemoryManager.reset(resizeToTotal); - dMemoryManager.reset(resizeToTotal); cMemoryManager.reset(resizeToTotal); } void Package::clearUniqueTables() { vUniqueTable.clear(); mUniqueTable.clear(); - dUniqueTable.clear(); cUniqueTable.clear(); } bool Package::garbageCollect(bool force) { - using flags = std::tuple; + using flags = std::tuple; // return immediately if no table needs collection if (!force && !vUniqueTable.possiblyNeedsCollection() && !mUniqueTable.possiblyNeedsCollection() && - !dUniqueTable.possiblyNeedsCollection() && !cUniqueTable.possiblyNeedsCollection()) { return false; } @@ -109,11 +102,10 @@ bool Package::garbageCollect(bool force) { force |= invC; const bool invV = vUniqueTable.garbageCollect(force) > 0; const bool invM = mUniqueTable.garbageCollect(force) > 0; - const bool invD = dUniqueTable.garbageCollect(force) > 0; - return {invC, invV, invM, invD}; + return {invC, invV, invM}; }; - const auto [invC, invV, invM, invD] = roots.execute(sweep); + const auto [invC, invV, invM] = roots.execute(sweep); // invalidate all compute tables involving vectors if any vector node has // been collected @@ -132,15 +124,6 @@ bool Package::garbageCollect(bool force) { matrixTrace.clear(); matrixVectorMultiplication.clear(); matrixMatrixMultiplication.clear(); - stochasticNoiseOperationCache.clear(); - } - // invalidate all compute tables involving density matrices if any density - // matrix node has been collected - if (invD) { - densityAdd.clear(); - densityDensityMultiplication.clear(); - densityNoise.clear(); - densityTrace.clear(); } // invalidate all compute tables where any component of the entry contains // numbers from the complex table if any complex numbers were collected @@ -152,35 +135,19 @@ bool Package::garbageCollect(bool force) { vectorKronecker.clear(); matrixKronecker.clear(); matrixTrace.clear(); - stochasticNoiseOperationCache.clear(); - densityAdd.clear(); - densityDensityMultiplication.clear(); - densityNoise.clear(); - densityTrace.clear(); } - return invC || invV || invM || invD; + return invC || invV || invM; } Package::ActiveCounts Package::computeActiveCounts() { const auto count = [this]() -> ActiveCounts { return {.vector = vUniqueTable.countMarkedEntries(), .matrix = mUniqueTable.countMarkedEntries(), - .density = dUniqueTable.countMarkedEntries(), .reals = cUniqueTable.countMarkedEntries()}; }; return roots.execute(count); } -dEdge Package::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; -} - mEdge Package::makeGateDD(const GateMatrix& mat, const qc::Qubit target) { return makeGateDD(mat, qc::Controls{}, target); } @@ -437,12 +404,6 @@ void Package::clearComputeTables() { vectorKronecker.clear(); matrixKronecker.clear(); matrixTrace.clear(); - - stochasticNoiseOperationCache.clear(); - densityAdd.clear(); - densityDensityMultiplication.clear(); - densityNoise.clear(); - densityTrace.clear(); } std::string Package::measureAll(vEdge& rootEdge, const bool collapse, std::mt19937_64& mt, const fp epsilon) { @@ -605,39 +566,6 @@ char Package::measureOneCollapsing(vEdge& rootEdge, const Qubit index, performCollapsingMeasurement(rootEdge, index, pone, false); return '1'; } -char Package::measureOneCollapsing(dEdge& e, const 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 = makeGateDD(MEAS_ZERO_MAT, index); - - auto tmp0 = conjugateTranspose(measZeroDd); - auto tmp1 = multiply(e, densityFromMatrixEdge(tmp0), false); - auto tmp2 = multiply(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 = makeGateDD(MEAS_ONE_MAT, index); - tmp0 = conjugateTranspose(measOneDd); - tmp1 = multiply(e, densityFromMatrixEdge(tmp0), false); - tmp2 = multiply(densityFromMatrixEdge(measOneDd), tmp1, true); - measuredResult = '1'; - densityMatrixTrace = trace(tmp2, nrQubits); - } - - dEdge::alignDensityEdge(e); - tmp2.w = cn.lookup(e.w / densityMatrixTrace); // Normalize density matrix - incRef(tmp2); - decRef(e); - e = tmp2; - dEdge::setDensityMatrixTrue(e); - - return measuredResult; -} void Package::performCollapsingMeasurement(vEdge& rootEdge, const Qubit index, const fp probability, const bool measureZero) { @@ -725,17 +653,6 @@ MatrixDD Package::applyOperation(const MatrixDD& operation, const MatrixDD& e, garbageCollect(); return tmp; } -dEdge Package::applyOperationToDensity(dEdge& e, const mEdge& operation) { - const auto tmp0 = conjugateTranspose(operation); - const auto tmp1 = multiply(e, densityFromMatrixEdge(tmp0), false); - const auto tmp2 = multiply(densityFromMatrixEdge(operation), tmp1, true); - incRef(tmp2); - dEdge::alignDensityEdge(e); - decRef(e); - e = tmp2; - dEdge::setDensityMatrixTrue(e); - return e; -} ComplexValue Package::innerProduct(const vEdge& x, const vEdge& y) { if (x.isTerminal() || y.isTerminal() || x.w.approximatelyZero() || y.w.approximatelyZero()) { // the 0 case @@ -857,6 +774,13 @@ mEdge Package::partialTrace(const mEdge& a, auto r = trace(a, eliminate, eliminate.size()); return {r.p, cn.lookup(r.w)}; } +ComplexValue Package::trace(const mEdge& 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; +} bool Package::isCloseToIdentity(const mEdge& m, const fp tol, const std::vector& garbage, const bool checkCloseToOne) const { @@ -864,6 +788,71 @@ bool Package::isCloseToIdentity(const mEdge& m, const fp tol, visited.reserve(mUniqueTable.getNumEntries()); return isCloseToIdentityRecursive(m, visited, tol, garbage, checkCloseToOne); } +mCachedEdge Package::trace(const mEdge& a, const std::vector& eliminate, + std::size_t level, std::size_t alreadyEliminated) { + const auto aWeight = static_cast(a.w); + if (aWeight.approximatelyZero()) { + return mCachedEdge::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 mCachedEdge{a.p, aWeight}; + } + + const auto v = a.p->v; + if (eliminate[v]) { + // Lookup nodes marked for elimination in the compute table if all + // lower-level qubits are eliminated as well: if the trace has already + // been computed, return the result + 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 = getTraceComputeTable().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); + + // The resulting weight is continuously normalized to the range [0,1] for + // matrix nodes + r.w = r.w / 2.0; + + // Insert result into compute table if all lower-level qubits are + // eliminated as well + if (eliminateAll) { + getTraceComputeTable().insert(a.p, r); + } + r.w = r.w * aWeight; + return r; + } + + std::array edge{}; + std::ranges::transform(std::as_const(a.p->e), edge.begin(), + [this, &eliminate, &alreadyEliminated, + &level](const mEdge& e) -> mCachedEdge { + return trace(e, eliminate, level - 1, + alreadyEliminated); + }); + const auto adjustedV = static_cast( + static_cast(a.p->v) - + (static_cast(std::ranges::count(eliminate, true)) - + alreadyEliminated)); + auto r = makeDDNode(adjustedV, edge); + r.w = r.w * aWeight; + return r; +} bool Package::isCloseToIdentityRecursive( const mEdge& m, std::unordered_set& visited, const fp tol, const std::vector& garbage, const bool checkCloseToOne) { diff --git a/src/dd/StateGeneration.cpp b/src/dd/StateGeneration.cpp index d52d5db11c..578fde5971 100644 --- a/src/dd/StateGeneration.cpp +++ b/src/dd/StateGeneration.cpp @@ -275,7 +275,8 @@ VectorDD generateExponentialState(const std::size_t levels, Package& dd, std::vector nodesPerLevel(levels); // [1, 2, 4, 8, ...] std::ranges::generate(nodesPerLevel, [exp = 0]() mutable { return 1ULL << exp++; }); - return generateRandomState(levels, nodesPerLevel, ROUNDROBIN, dd, seed); + return generateRandomState(levels, nodesPerLevel, + GenerationWireStrategy::ROUNDROBIN, dd, seed); } VectorDD generateRandomState(const std::size_t levels, @@ -324,12 +325,12 @@ VectorDD generateRandomState(const std::size_t levels, std::vector indices(2 * n); // Indices for wireing. switch (strategy) { - case ROUNDROBIN: { + case GenerationWireStrategy::ROUNDROBIN: { std::ranges::generate(indices, [&m, r = 0UL]() mutable { return (r++) % m; }); break; } - case RANDOM: { + case GenerationWireStrategy::RANDOM: { IndexDistribution idxDist{0, m - 1}; // Ensure that all the nodes below have a connection upwards. diff --git a/src/dd/statistics/PackageStatistics.cpp b/src/dd/statistics/PackageStatistics.cpp index bc10b801d5..01140ea3b1 100644 --- a/src/dd/statistics/PackageStatistics.cpp +++ b/src/dd/statistics/PackageStatistics.cpp @@ -28,8 +28,6 @@ static constexpr auto V_NODE_MEMORY_MIB = static_cast(sizeof(vNode)) / static_cast(1ULL << 20U); static constexpr auto M_NODE_MEMORY_MIB = static_cast(sizeof(mNode)) / static_cast(1ULL << 20U); -static constexpr auto D_NODE_MEMORY_MIB = - static_cast(sizeof(dNode)) / static_cast(1ULL << 20U); static constexpr auto REAL_NUMBER_MEMORY_MIB = static_cast(sizeof(RealNumber)) / static_cast(1ULL << 20U); @@ -37,27 +35,20 @@ static constexpr auto V_EDGE_MEMORY_MIB = static_cast(sizeof(Edge)) / static_cast(1ULL << 20U); static constexpr auto M_EDGE_MEMORY_MIB = static_cast(sizeof(Edge)) / static_cast(1ULL << 20U); -static constexpr auto D_EDGE_MEMORY_MIB = - static_cast(sizeof(Edge)) / static_cast(1ULL << 20U); double computeActiveMemoryMiB(Package& package) { - const auto [vectorNodes, matrixNodes, densityNodes, realNumbers] = + const auto [vectorNodes, matrixNodes, realNumbers] = package.computeActiveCounts(); const auto vActiveEntries = static_cast(vectorNodes); const auto mActiveEntries = static_cast(matrixNodes); - const auto dActiveEntries = static_cast(densityNodes); const auto vMemoryForNodes = vActiveEntries * V_NODE_MEMORY_MIB; const auto mMemoryForNodes = mActiveEntries * M_NODE_MEMORY_MIB; - const auto dMemoryForNodes = dActiveEntries * D_NODE_MEMORY_MIB; - const auto memoryForNodes = - vMemoryForNodes + mMemoryForNodes + dMemoryForNodes; + const auto memoryForNodes = vMemoryForNodes + mMemoryForNodes; const auto vMemoryForEdges = vActiveEntries * V_EDGE_MEMORY_MIB; const auto mMemoryForEdges = mActiveEntries * M_EDGE_MEMORY_MIB; - const auto dMemoryForEdges = dActiveEntries * D_EDGE_MEMORY_MIB; - const auto memoryForEdges = - vMemoryForEdges + mMemoryForEdges + dMemoryForEdges; + const auto memoryForEdges = vMemoryForEdges + mMemoryForEdges; const auto activeRealNumbers = static_cast(realNumbers); const auto memoryForRealNumbers = activeRealNumbers * REAL_NUMBER_MEMORY_MIB; @@ -70,20 +61,14 @@ double computePeakMemoryMiB(const Package& package) { static_cast(package.vMemoryManager.getStats().peakNumUsed); const auto mPeakUsedEntries = static_cast(package.mMemoryManager.getStats().peakNumUsed); - const auto dPeakUsedEntries = - static_cast(package.dMemoryManager.getStats().peakNumUsed); const auto vMemoryForNodes = vPeakUsedEntries * V_NODE_MEMORY_MIB; const auto mMemoryForNodes = mPeakUsedEntries * M_NODE_MEMORY_MIB; - const auto dMemoryForNodes = dPeakUsedEntries * D_NODE_MEMORY_MIB; - const auto memoryForNodes = - vMemoryForNodes + mMemoryForNodes + dMemoryForNodes; + const auto memoryForNodes = vMemoryForNodes + mMemoryForNodes; const auto vMemoryForEdges = vPeakUsedEntries * V_EDGE_MEMORY_MIB; const auto mMemoryForEdges = mPeakUsedEntries * M_EDGE_MEMORY_MIB; - const auto dMemoryForEdges = dPeakUsedEntries * D_EDGE_MEMORY_MIB; - const auto memoryForEdges = - vMemoryForEdges + mMemoryForEdges + dMemoryForEdges; + const auto memoryForEdges = vMemoryForEdges + mMemoryForEdges; const auto peakRealNumbers = static_cast(package.cMemoryManager.getStats().peakNumUsed); @@ -98,8 +83,8 @@ nlohmann::basic_json<> getStatistics(Package& package, j["data_structure"] = getDataStructureStatistics(); - const auto [activeVectorNodes, activeMatrixNodes, activeDensityNodes, - activeRealNumbers] = package.computeActiveCounts(); + const auto [activeVectorNodes, activeMatrixNodes, activeRealNumbers] = + package.computeActiveCounts(); auto& vector = j["vector"]; auto& vectorUniqueTable = vector["unique_table"]; @@ -119,15 +104,6 @@ nlohmann::basic_json<> getStatistics(Package& package, } matrix["memory_manager"] = package.mMemoryManager.getStats().json(); - auto& densityMatrix = j["density_matrix"]; - auto& densityUniqueTable = densityMatrix["unique_table"]; - densityUniqueTable = - package.dUniqueTable.getStatsJson(includeIndividualTables); - if (densityUniqueTable != "unused") { - densityUniqueTable["total"]["num_active_entries"] = activeDensityNodes; - } - densityMatrix["memory_manager"] = package.dMemoryManager.getStats().json(); - auto& realNumbers = j["real_numbers"]; auto& realNumbersUniqueTable = realNumbers["unique_table"]; realNumbersUniqueTable = package.cUniqueTable.getStats().json(); @@ -139,23 +115,16 @@ nlohmann::basic_json<> getStatistics(Package& package, auto& computeTables = j["compute_tables"]; computeTables["vector_add"] = package.vectorAdd.getStats().json(); computeTables["matrix_add"] = package.matrixAdd.getStats().json(); - computeTables["density_matrix_add"] = package.densityAdd.getStats().json(); computeTables["matrix_conjugate_transpose"] = package.conjugateMatrixTranspose.getStats().json(); computeTables["matrix_vector_mult"] = package.matrixVectorMultiplication.getStats().json(); computeTables["matrix_matrix_mult"] = package.matrixMatrixMultiplication.getStats().json(); - computeTables["density_density_mult"] = - package.densityDensityMultiplication.getStats().json(); computeTables["vector_kronecker"] = package.vectorKronecker.getStats().json(); computeTables["matrix_kronecker"] = package.matrixKronecker.getStats().json(); computeTables["vector_inner_product"] = package.vectorInnerProduct.getStats().json(); - computeTables["stochastic_noise_operations"] = - package.stochasticNoiseOperationCache.getStats().json(); - computeTables["density_noise_operations"] = - package.densityNoise.getStats().json(); j["active_memory_mib"] = computeActiveMemoryMiB(package); j["peak_memory_mib"] = computePeakMemoryMiB(package); @@ -180,10 +149,6 @@ nlohmann::basic_json<> getDataStructureStatistics() { matrixNode["size_B"] = sizeof(mNode); matrixNode["alignment_B"] = alignof(mNode); - auto& densityNode = j["dNode"]; - densityNode["size_B"] = sizeof(dNode); - densityNode["alignment_B"] = alignof(dNode); - auto& vectorEdge = j["vEdge"]; vectorEdge["size_B"] = sizeof(Edge); vectorEdge["alignment_B"] = alignof(Edge); @@ -192,10 +157,6 @@ nlohmann::basic_json<> getDataStructureStatistics() { matrixEdge["size_B"] = sizeof(Edge); matrixEdge["alignment_B"] = alignof(Edge); - auto& densityEdge = j["dEdge"]; - densityEdge["size_B"] = sizeof(Edge); - densityEdge["alignment_B"] = alignof(Edge); - auto& realNumber = j["RealNumber"]; realNumber["size_B"] = sizeof(RealNumber); realNumber["alignment_B"] = alignof(RealNumber); @@ -225,11 +186,6 @@ nlohmann::basic_json<> getDataStructureStatistics() { matrixAdd["alignment_B"] = alignof(typename decltype(Package::matrixAdd)::Entry); - auto& densityAdd = ctEntries["density_add"]; - densityAdd["size_B"] = sizeof(typename decltype(Package::densityAdd)::Entry); - densityAdd["alignment_B"] = - alignof(typename decltype(Package::densityAdd)::Entry); - auto& conjugateMatrixTranspose = ctEntries["conjugate_matrix_transpose"]; conjugateMatrixTranspose["size_B"] = sizeof(typename decltype(Package::conjugateMatrixTranspose)::Entry); @@ -248,12 +204,6 @@ nlohmann::basic_json<> getDataStructureStatistics() { matrixMatrixMult["alignment_B"] = alignof(typename decltype(Package::matrixMatrixMultiplication)::Entry); - auto& densityDensityMult = ctEntries["density_density_mult"]; - densityDensityMult["size_B"] = - sizeof(typename decltype(Package::densityDensityMultiplication)::Entry); - densityDensityMult["alignment_B"] = - alignof(typename decltype(Package::densityDensityMultiplication)::Entry); - auto& vectorKronecker = ctEntries["vector_kronecker"]; vectorKronecker["size_B"] = sizeof(typename decltype(Package::vectorKronecker)::Entry); diff --git a/test/dd/test_approximations.cpp b/test/dd/test_approximations.cpp index 7df6415f42..a0553cebd6 100644 --- a/test/dd/test_approximations.cpp +++ b/test/dd/test_approximations.cpp @@ -29,7 +29,7 @@ namespace { /** * @brief Compare the elements of @p a and @p b with precision @p delta. */ -void vecNear(CVec a, CVec b, double delta = 1e-6) { +void vecNear(const CVec& a, const CVec& b, const double delta = 1e-6) { for (std::size_t i = 0; i < b.size(); ++i) { EXPECT_NEAR(a[i].real(), b[i].real(), delta); EXPECT_NEAR(a[i].imag(), b[i].imag(), delta); @@ -60,7 +60,7 @@ TEST(ApproximationTest, OneQubitKeepAllBudgetZero) { constexpr std::size_t nq = 1; constexpr double fidelity = 1; - auto dd = std::make_unique(nq); + auto dd = std::make_unique(nq); qc::QuantumComputation qc(nq); qc.x(0); @@ -93,7 +93,7 @@ TEST(ApproximationTest, OneQubitKeepAllBudgetTooSmall) { constexpr std::size_t nq = 1; constexpr double fidelity = 0.9; - auto dd = std::make_unique(nq); + auto dd = std::make_unique(nq); qc::QuantumComputation qc(nq); qc.x(0); @@ -126,7 +126,7 @@ TEST(ApproximationTest, OneQubitRemoveTerminalEdge) { constexpr std::size_t nq = 1; constexpr double fidelity = 1 - 0.25; - auto dd = std::make_unique(nq); + auto dd = std::make_unique(nq); qc::QuantumComputation qc(nq); qc.ry(qc::PI / 3, 0); @@ -163,7 +163,7 @@ TEST(ApproximationTest, TwoQubitRemoveNode) { constexpr std::size_t nq = 2; constexpr double fidelity = 1 - 0.2; - auto dd = std::make_unique(nq); + auto dd = std::make_unique(nq); qc::QuantumComputation qc(nq); qc.h(0); @@ -211,7 +211,7 @@ TEST(ApproximationTest, TwoQubitCorrectlyRebuilt) { constexpr std::size_t nq = 2; constexpr double fidelity = 1 - 0.25; - auto dd = std::make_unique(nq); + auto dd = std::make_unique(nq); qc::QuantumComputation qc(nq); qc.rx(qc::PI, 0); @@ -264,7 +264,7 @@ TEST(ApproximationTest, ThreeQubitRemoveNodeWithChildren) { constexpr std::size_t nq = 3; constexpr double fidelity = 1 - 0.25; - auto dd = std::make_unique(nq); + auto dd = std::make_unique(nq); qc::QuantumComputation qc(nq); qc.rx(qc::PI, 0); @@ -310,7 +310,7 @@ TEST(ApproximationTest, ThreeQubitRemoveUnconnected) { constexpr std::size_t nq = 3; constexpr double fidelity = 1 - 0.65; - auto dd = std::make_unique(nq); + auto dd = std::make_unique(nq); qc::QuantumComputation qc(nq); qc.ry(2 * qc::PI / 3, 0); @@ -339,7 +339,7 @@ TEST(ApproximationTest, NodesVisited) { for (std::size_t i = 0; i < n; ++i) { const std::size_t nq = qubits[i]; - auto dd = std::make_unique(nq); + auto dd = std::make_unique(nq); auto state = generateExponentialState(nq, *dd); const std::size_t preSize = state.size() - 1; // Minus terminal. diff --git a/test/dd/test_dd_functionality.cpp b/test/dd/test_dd_functionality.cpp index 508dc09e44..11a4f02ae5 100644 --- a/test/dd/test_dd_functionality.cpp +++ b/test/dd/test_dd_functionality.cpp @@ -328,11 +328,10 @@ TEST_F(DDFunctionality, BuildCircuit) { dd->decRef(dd2); dd->garbageCollect(true); - const auto counts = dd->computeActiveCounts(); - EXPECT_EQ(counts.vector, 0); - EXPECT_EQ(counts.density, 0); - EXPECT_EQ(counts.matrix, 0); - EXPECT_EQ(counts.reals, 0); + const auto [vector, matrix, reals] = dd->computeActiveCounts(); + EXPECT_EQ(vector, 0); + EXPECT_EQ(matrix, 0); + EXPECT_EQ(reals, 0); } TEST_F(DDFunctionality, NonUnitary) { @@ -386,11 +385,10 @@ TEST_F(DDFunctionality, CircuitEquivalence) { dd->decRef(dd2); dd->garbageCollect(true); - const auto counts = dd->computeActiveCounts(); - EXPECT_EQ(counts.vector, 0); - EXPECT_EQ(counts.density, 0); - EXPECT_EQ(counts.matrix, 0); - EXPECT_EQ(counts.reals, 0); + const auto [vector, matrix, reals] = dd->computeActiveCounts(); + EXPECT_EQ(vector, 0); + EXPECT_EQ(matrix, 0); + EXPECT_EQ(reals, 0); } TEST_F(DDFunctionality, ChangePermutation) { @@ -443,11 +441,10 @@ TEST_F(DDFunctionality, FuseTwoSingleQubitGates) { dd->decRef(optDD); dd->garbageCollect(true); - const auto counts = dd->computeActiveCounts(); - EXPECT_EQ(counts.vector, 0); - EXPECT_EQ(counts.density, 0); - EXPECT_EQ(counts.matrix, 0); - EXPECT_EQ(counts.reals, 0); + const auto [vector, matrix, reals] = dd->computeActiveCounts(); + EXPECT_EQ(vector, 0); + EXPECT_EQ(matrix, 0); + EXPECT_EQ(reals, 0); } TEST_F(DDFunctionality, FuseThreeSingleQubitGates) { @@ -478,11 +475,10 @@ TEST_F(DDFunctionality, FuseThreeSingleQubitGates) { dd->decRef(optDD); dd->garbageCollect(true); - const auto counts = dd->computeActiveCounts(); - EXPECT_EQ(counts.vector, 0); - EXPECT_EQ(counts.density, 0); - EXPECT_EQ(counts.matrix, 0); - EXPECT_EQ(counts.reals, 0); + const auto [vector, matrix, reals] = dd->computeActiveCounts(); + EXPECT_EQ(vector, 0); + EXPECT_EQ(matrix, 0); + EXPECT_EQ(reals, 0); } TEST_F(DDFunctionality, FuseNoSingleQubitGates) { @@ -513,11 +509,10 @@ TEST_F(DDFunctionality, FuseNoSingleQubitGates) { dd->decRef(optDD); dd->garbageCollect(true); - const auto counts = dd->computeActiveCounts(); - EXPECT_EQ(counts.vector, 0); - EXPECT_EQ(counts.density, 0); - EXPECT_EQ(counts.matrix, 0); - EXPECT_EQ(counts.reals, 0); + const auto [vector, matrix, reals] = dd->computeActiveCounts(); + EXPECT_EQ(vector, 0); + EXPECT_EQ(matrix, 0); + EXPECT_EQ(reals, 0); } TEST_F(DDFunctionality, FuseSingleQubitGatesAcrossOtherGates) { @@ -547,11 +542,10 @@ TEST_F(DDFunctionality, FuseSingleQubitGatesAcrossOtherGates) { dd->decRef(optDD); dd->garbageCollect(true); - const auto counts = dd->computeActiveCounts(); - EXPECT_EQ(counts.vector, 0); - EXPECT_EQ(counts.density, 0); - EXPECT_EQ(counts.matrix, 0); - EXPECT_EQ(counts.reals, 0); + const auto [vector, matrix, reals] = dd->computeActiveCounts(); + EXPECT_EQ(vector, 0); + EXPECT_EQ(matrix, 0); + EXPECT_EQ(reals, 0); } TEST_F(DDFunctionality, IfElseOperationConditions) { @@ -612,11 +606,10 @@ TEST_F(DDFunctionality, VectorKroneckerWithTerminal) { dd->decRef(zeroState); dd->garbageCollect(true); - const auto counts = dd->computeActiveCounts(); - EXPECT_EQ(counts.vector, 0); - EXPECT_EQ(counts.density, 0); - EXPECT_EQ(counts.matrix, 0); - EXPECT_EQ(counts.reals, 0); + const auto [vector, matrix, reals] = dd->computeActiveCounts(); + EXPECT_EQ(vector, 0); + EXPECT_EQ(matrix, 0); + EXPECT_EQ(reals, 0); } TEST_F(DDFunctionality, DynamicCircuitSimulationWithSWAP) { diff --git a/test/dd/test_dd_noise_functionality.cpp b/test/dd/test_dd_noise_functionality.cpp deleted file mode 100644 index b87780730b..0000000000 --- a/test/dd/test_dd_noise_functionality.cpp +++ /dev/null @@ -1,344 +0,0 @@ -/* - * 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 "dd/DDDefinitions.hpp" -#include "dd/DDpackageConfig.hpp" -#include "dd/NoiseFunctionality.hpp" -#include "dd/Operations.hpp" -#include "dd/Package.hpp" -#include "dd/StateGeneration.hpp" -#include "ir/Definitions.hpp" -#include "ir/QuantumComputation.hpp" -#include "ir/operations/IfElseOperation.hpp" -#include "ir/operations/OpType.hpp" - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -using namespace qc; - -class DDNoiseFunctionalityTest : public ::testing::Test { -protected: - void SetUp() override { - using namespace qc::literals; - - // 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; - 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::DENSITY_MATRIX_SIMULATOR_DD_PACKAGE_CONFIG); - - auto rootEdge = dd->makeZeroDensityOperator(qc.getNqubits()); - - const auto* const noiseEffects = "APDI"; - - auto deterministicNoiseFunctionality = dd::DeterministicNoiseFunctionality( - *dd, qc.getNqubits(), 0.01, 0.02, 0.02, 0.04, noiseEffects); - - for (auto const& op : qc) { - dd->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 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::DENSITY_MATRIX_SIMULATOR_DD_PACKAGE_CONFIG); - - auto rootEdge = dd->makeZeroDensityOperator(qc.getNqubits()); - - const auto* const noiseEffects = "D"; - - auto deterministicNoiseFunctionality = dd::DeterministicNoiseFunctionality( - *dd, qc.getNqubits(), 0.01, 0.02, 0.02, 0.04, noiseEffects); - - for (auto const& op : qc) { - dd->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 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::DENSITY_MATRIX_SIMULATOR_DD_PACKAGE_CONFIG); - - auto rootEdge = dd->makeZeroDensityOperator(qcOp.getNqubits()); - - auto deterministicNoiseFunctionality = dd::DeterministicNoiseFunctionality( - *dd, qcOp.getNqubits(), 0.01, 0.02, 0.02, 0.04, {}); - - for (auto const& op : qcOp) { - dd->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); - - dd->measureOneCollapsing(rootEdge, 0, qc.getGenerator()); - - auto tmp0 = rootEdge.getSparseProbabilityVectorStrKeys(qc.getNqubits()); - prob = 0.25; - - EXPECT_TRUE(fabs(tmp0["000"] + tmp0["001"] - prob) < tolerance); - EXPECT_TRUE(fabs(tmp0["010"] + tmp0["011"] - prob) < tolerance); - EXPECT_TRUE(fabs(tmp0["100"] + tmp0["101"] - prob) < tolerance); - EXPECT_TRUE(fabs(tmp0["110"] + tmp0["111"] - prob) < tolerance); - - dd->measureOneCollapsing(rootEdge, 1, qc.getGenerator()); - - auto tmp1 = rootEdge.getSparseProbabilityVectorStrKeys(qc.getNqubits()); - prob = 0.5; - EXPECT_TRUE(fabs(tmp0["000"] + tmp0["001"] + tmp0["010"] + tmp0["011"] - - prob) < tolerance); - EXPECT_TRUE(fabs(tmp0["100"] + tmp0["101"] + tmp0["110"] + tmp0["111"] - - prob) < tolerance); - - dd->measureOneCollapsing(rootEdge, 2, qc.getGenerator()); - - auto tmp2 = rootEdge.getSparseProbabilityVectorStrKeys(qc.getNqubits()); - EXPECT_TRUE( - fabs(tmp2["000"] - 1) < tolerance || fabs(tmp2["001"] - 1) < tolerance || - fabs(tmp2["010"] - 1) < tolerance || fabs(tmp2["011"] - 1) < tolerance || - fabs(tmp2["100"] - 1) < tolerance || fabs(tmp2["101"] - 1) < tolerance || - fabs(tmp2["111"] - 1) < tolerance); -} - -TEST_F(DDNoiseFunctionalityTest, StochSimulateAdder4TrackAPD) { - auto dd = std::make_unique( - qc.getNqubits(), dd::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::StochasticNoiseFunctionality( - *dd, qc.getNqubits(), 0.01, 0.02, 2., noiseEffects); - - for (size_t i = 0U; i < stochRuns; i++) { - auto rootEdge = 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 (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::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::StochasticNoiseFunctionality( - *dd, qc.getNqubits(), 0.01, 0.02, 2., noiseEffects); - - for (size_t i = 0U; i < stochRuns; i++) { - auto rootEdge = makeZeroState(qc.getNqubits(), *dd); - dd->incRef(rootEdge); - - for (auto const& op : qc) { - auto operation = dd::getDD(*op, *dd); - auto usedQubits = op->getUsedQubits(); - stochasticNoiseFunctionality.applyNoiseOperation( - op->getUsedQubits(), operation, rootEdge, qc.getGenerator()); - } - - const auto amplitudes = rootEdge.getVector(); - for (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, testingUsedQubits) { - auto standardOp = StandardOperation(1, qc::Z); - EXPECT_EQ(standardOp.getUsedQubits().size(), 1); - EXPECT_TRUE(standardOp.getUsedQubits().count(1)); - - auto nonUnitaryOp = NonUnitaryOperation(0, 0); - EXPECT_EQ(nonUnitaryOp.getUsedQubits().size(), 1); - EXPECT_TRUE(nonUnitaryOp.getUsedQubits().count(0) == 1U); - - auto compoundOp = qc::CompoundOperation(); - compoundOp.emplace_back(0, qc::Z); - compoundOp.emplace_back(1, qc::H); - compoundOp.emplace_back(0, qc::X); - EXPECT_EQ(compoundOp.getUsedQubits().size(), 2); - EXPECT_TRUE(compoundOp.getUsedQubits().count(0)); - EXPECT_TRUE(compoundOp.getUsedQubits().count(1)); - - auto ifElseOp = qc::IfElseOperation( - std::make_unique(0, qc::X), nullptr, 0); - auto* thenOp = ifElseOp.getThenOp(); - EXPECT_EQ(thenOp->getUsedQubits().size(), 1); - EXPECT_TRUE(thenOp->getUsedQubits().count(0) == 1U); -} - -TEST_F(DDNoiseFunctionalityTest, invalidNoiseEffect) { - auto dd = std::make_unique( - qc.getNqubits(), dd::STOCHASTIC_NOISE_SIMULATOR_DD_PACKAGE_CONFIG); - EXPECT_THROW(dd::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::STOCHASTIC_NOISE_SIMULATOR_DD_PACKAGE_CONFIG); - EXPECT_THROW(dd::StochasticNoiseFunctionality(*dd, qc.getNqubits(), 0.3, 0.6, - 2, "APD"), - std::runtime_error); -} diff --git a/test/dd/test_edge_functionality.cpp b/test/dd/test_edge_functionality.cpp index 40875384e0..118cdc1b7d 100644 --- a/test/dd/test_edge_functionality.cpp +++ b/test/dd/test_edge_functionality.cpp @@ -10,20 +10,14 @@ #include "dd/DDDefinitions.hpp" #include "dd/Node.hpp" -#include "dd/Operations.hpp" #include "dd/Package.hpp" #include "dd/RealNumber.hpp" #include "dd/StateGeneration.hpp" -#include "ir/operations/OpType.hpp" -#include "ir/operations/StandardOperation.hpp" #include #include #include -#include #include -#include -#include namespace dd { @@ -360,113 +354,4 @@ TEST(MatrixFunctionality, SizeBellState) { EXPECT_EQ(bell.size(), 3); } -///----------------------------------------------------------------------------- -/// \n Tests for density matrix DDs \n -///----------------------------------------------------------------------------- - -TEST(DensityMatrixFunctionality, GetValueByPathTerminal) { - EXPECT_EQ(dEdge::zero().getValueByPath(0, "0"), 0.); - EXPECT_EQ(dEdge::one().getValueByPath(0, "0"), 1.); -} - -TEST(DensityMatrixFunctionality, GetValueByIndexTerminal) { - EXPECT_EQ(dEdge::zero().getValueByIndex(0, 0, 0), 0.); - EXPECT_EQ(dEdge::one().getValueByIndex(0, 0, 0), 1.); -} - -TEST(DensityMatrixFunctionality, GetValueByIndexProperDensityMatrix) { - const auto nqubits = 1U; - auto dd = std::make_unique(nqubits); - auto zero = dd->makeZeroDensityOperator(nqubits); - const auto op1 = getDD(qc::StandardOperation(0U, qc::H), *dd); - const auto op2 = getDD(qc::StandardOperation(0, qc::RZ, {PI_4}), *dd); - auto state = dd->applyOperationToDensity(zero, op1); - state = dd->applyOperationToDensity(state, op2); - - const auto diagValRef = 0.5; - const auto offDiagValRef = 0.25 * std::numbers::sqrt2; - - const CMat dmRef = {{{diagValRef, 0.}, {offDiagValRef, -offDiagValRef}}, - {{offDiagValRef, offDiagValRef}, {diagValRef, 0.}}}; - - const auto dm = state.getMatrix(nqubits); - - for (std::size_t i = 0U; i < dm.size(); ++i) { - for (std::size_t j = 0U; j < dm.size(); ++j) { - const auto val = state.getValueByIndex(nqubits, i, j); - const auto ref = dmRef[i][j]; - EXPECT_NEAR(ref.real(), val.real(), 1e-10); - EXPECT_NEAR(ref.imag(), val.imag(), 1e-10); - } - } -} - -TEST(DensityMatrixFunctionality, GetSparseMatrixTerminal) { - const auto zero = SparseCMat{{{0, 0}, 0.}}; - EXPECT_EQ(dEdge::zero().getSparseMatrix(0), zero); - const auto one = SparseCMat{{{0, 0}, 1.}}; - EXPECT_EQ(dEdge::one().getSparseMatrix(0), one); -} - -TEST(DensityMatrixFunctionality, GetSparseMatrixConsistency) { - const auto nqubits = 1U; - auto dd = std::make_unique(nqubits); - auto zero = dd->makeZeroDensityOperator(nqubits); - const auto op1 = getDD(qc::StandardOperation(0U, qc::H), *dd); - const auto op2 = getDD(qc::StandardOperation(0, qc::RZ, {PI_4}), *dd); - auto state = dd->applyOperationToDensity(zero, op1); - state = dd->applyOperationToDensity(state, op2); - - const auto dm = state.getSparseMatrix(nqubits); - const auto dmDense = state.getMatrix(nqubits); - - for (const auto& [index, value] : dm) { - const auto val = dmDense.at(index.first).at(index.second); - EXPECT_NEAR(value.real(), val.real(), 1e-10); - EXPECT_NEAR(value.imag(), val.imag(), 1e-10); - } -} - -TEST(DensityMatrixFunctionality, PrintMatrixTerminal) { - testing::internal::CaptureStdout(); - dEdge::zero().printMatrix(0); - const auto zeroStr = testing::internal::GetCapturedStdout(); - EXPECT_EQ(zeroStr, "(0,0)\n"); - testing::internal::CaptureStdout(); - dEdge::one().printMatrix(0); - const auto oneStr = testing::internal::GetCapturedStdout(); - EXPECT_EQ(oneStr, "(1,0)\n"); -} - -TEST(DensityMatrixFunctionality, PrintMatrix) { - const auto nqubits = 1U; - auto dd = std::make_unique(nqubits); - auto zero = dd->makeZeroDensityOperator(nqubits); - const auto op1 = getDD(qc::StandardOperation(0U, qc::H), *dd); - const auto op2 = getDD(qc::StandardOperation(0, qc::RZ, {PI_4}), *dd); - auto state = dd->applyOperationToDensity(zero, op1); - state = dd->applyOperationToDensity(state, op2); - - const auto diagValRef = 0.5; - const auto offDiagValRef = 0.25 * std::numbers::sqrt2; - - const CMat dmRef = {{{diagValRef, 0.}, {offDiagValRef, -offDiagValRef}}, - {{offDiagValRef, offDiagValRef}, {diagValRef, 0.}}}; - - testing::internal::CaptureStdout(); - state.printMatrix(nqubits); - const auto matStr = testing::internal::GetCapturedStdout(); - - std::stringstream ss{}; - constexpr auto prec = 3U; - for (std::size_t i = 0U; i < dmRef.size(); ++i) { - for (std::size_t j = 0U; j < dmRef.size(); ++j) { - ss << std::setprecision(prec) << dmRef[i][j] << " "; - } - ss << "\n"; - } - - EXPECT_EQ(matStr, ss.str()); -} - } // namespace dd diff --git a/test/dd/test_package.cpp b/test/dd/test_package.cpp index e34251667b..355340e3b1 100644 --- a/test/dd/test_package.cpp +++ b/test/dd/test_package.cpp @@ -9,7 +9,6 @@ */ #include "dd/DDDefinitions.hpp" -#include "dd/DDpackageConfig.hpp" #include "dd/Export.hpp" #include "dd/GateMatrixDefinitions.hpp" #include "dd/MemoryManager.hpp" @@ -381,7 +380,7 @@ TEST(DDPackageTest, TraceComplexity) { // instead of paths in the DD due to the usage of a compute table for (std::size_t numQubits = 1; numQubits <= 10; ++numQubits) { auto dd = std::make_unique(numQubits); - auto& computeTable = dd->getTraceComputeTable(); + auto& computeTable = dd->getTraceComputeTable(); const auto hGate = getDD(qc::StandardOperation(0, qc::H), *dd); auto hKron = hGate; for (std::size_t i = 0; i < numQubits - 1; ++i) { @@ -944,27 +943,22 @@ TEST(DDPackageTest, PackageReset) { TEST(DDPackageTest, ResetClearsRoots) { auto dd = std::make_unique(2); - auto& vRoots = dd->getRootSet(); - auto& mRoots = dd->getRootSet(); - auto& dRoots = dd->getRootSet(); + const auto& vRoots = dd->getRootSet(); + const auto& mRoots = dd->getRootSet(); - auto vec = dd::makeZeroState(2, *dd); - auto mat = getDD(qc::StandardOperation(0, qc::X), *dd); - auto dens = dd->makeZeroDensityOperator(2); + const auto vec = dd::makeZeroState(2, *dd); + const auto mat = getDD(qc::StandardOperation(0, qc::X), *dd); dd->incRef(vec); dd->incRef(mat); - dd->incRef(dens); EXPECT_EQ(vRoots.size(), 1U); EXPECT_EQ(mRoots.size(), 1U); - EXPECT_EQ(dRoots.size(), 1U); (*dd).reset(); EXPECT_TRUE(vRoots.empty()); EXPECT_TRUE(mRoots.empty()); - EXPECT_TRUE(dRoots.empty()); } TEST(DDPackageTest, DuplicatetrackDoesNotLeaveStaleRoot) { @@ -972,7 +966,6 @@ TEST(DDPackageTest, DuplicatetrackDoesNotLeaveStaleRoot) { auto& vRoots = dd->getRootSet(); auto& mRoots = dd->getRootSet(); - auto& dRoots = dd->getRootSet(); // vector root auto vec = dd::makeZeroState(1, *dd); @@ -996,17 +989,6 @@ TEST(DDPackageTest, DuplicatetrackDoesNotLeaveStaleRoot) { dd->decRef(mat); EXPECT_TRUE(mRoots.empty()); EXPECT_THROW(dd->decRef(mat), std::invalid_argument); - - // density root - auto dens = dd->makeZeroDensityOperator(1); - EXPECT_EQ(dRoots.size(), 1U); - dd->incRef(dens); - EXPECT_EQ(dRoots.at(dens), 2U); - dd->decRef(dens); - EXPECT_EQ(dRoots.at(dens), 1U); - dd->decRef(dens); - EXPECT_TRUE(dRoots.empty()); - EXPECT_THROW(dd->decRef(dens), std::invalid_argument); } TEST(DDPackageTest, Inverse) { @@ -1568,168 +1550,6 @@ TEST(DDPackageTest, CloseToIdentityWithGarbageInTheMiddle) { {true, false, true}, false)); } -TEST(DDPackageTest, dNodeMultiply) { - // Multiply dNode with mNode (MxMxM) - constexpr auto nrQubits = 3U; - const auto dd = std::make_unique( - nrQubits, DENSITY_MATRIX_SIMULATOR_DD_PACKAGE_CONFIG); - // Make zero density matrix - auto state = dd->makeZeroDensityOperator(dd->qubits()); - std::vector operations = {}; - operations.emplace_back(getDD(qc::StandardOperation(0, qc::H), *dd)); - operations.emplace_back(getDD(qc::StandardOperation(1, qc::H), *dd)); - operations.emplace_back(getDD(qc::StandardOperation(2, qc::H), *dd)); - operations.emplace_back(getDD(qc::StandardOperation(2, qc::Z), *dd)); - - for (const auto& op : operations) { - dd->applyOperationToDensity(state, op); - } - - const auto stateDensityMatrix = state.getMatrix(dd->qubits()); - - for (const auto& stateVector : stateDensityMatrix) { - for (const auto& cValue : stateVector) { - std::cout << "r:" << cValue.real() << " i:" << cValue.imag(); - } - std::cout << "\n"; - } - - for (std::size_t i = 0; i < (1 << nrQubits); i++) { - for (std::size_t j = 0; j < (1 << nrQubits); j++) { - EXPECT_EQ(std::abs(stateDensityMatrix[i][j].imag()), 0); - if ((i < 4 && j < 4) || (i >= 4 && j >= 4)) { - EXPECT_TRUE(stateDensityMatrix[i][j].real() > 0); - } else { - EXPECT_TRUE(stateDensityMatrix[i][j].real() < 0); - } - EXPECT_TRUE(std::abs(std::abs(stateDensityMatrix[i][j]) - 0.125) < - 0.000001); - } - } - - const auto probVector = state.getSparseProbabilityVector(nrQubits, 0.001); - for (const auto& [s, prob] : probVector) { - constexpr auto tolerance = 1e-10; - std::cout << s << ": " << prob << "\n"; - EXPECT_NEAR(prob, 0.125, tolerance); - } -} - -TEST(DDPackageTest, dNodeMultiply2) { - // Multiply dNode with mNode (MxMxM) - constexpr auto nrQubits = 3U; - const auto dd = std::make_unique( - nrQubits, DENSITY_MATRIX_SIMULATOR_DD_PACKAGE_CONFIG); - // Make zero density matrix - auto state = dd->makeZeroDensityOperator(dd->qubits()); - std::vector operations = {}; - operations.emplace_back(getDD(qc::StandardOperation(0, qc::H), *dd)); - operations.emplace_back(getDD(qc::StandardOperation(1, qc::H), *dd)); - operations.emplace_back(getDD(qc::StandardOperation(2, qc::H), *dd)); - operations.emplace_back(getDD(qc::StandardOperation(2, qc::Z), *dd)); - - for (const auto& op : operations) { - dd->applyOperationToDensity(state, op); - } - operations[0].printMatrix(dd->qubits()); - - const auto stateDensityMatrix = state.getMatrix(dd->qubits()); - - for (std::size_t i = 0; i < (1 << nrQubits); i++) { - for (std::size_t j = 0; j < (1 << nrQubits); j++) { - EXPECT_TRUE(std::abs(stateDensityMatrix[i][j].imag()) == 0); - if ((i < 4 && j < 4) || (i >= 4 && j >= 4)) { - EXPECT_TRUE(stateDensityMatrix[i][j].real() > 0); - } else { - EXPECT_TRUE(stateDensityMatrix[i][j].real() < 0); - } - EXPECT_TRUE(std::abs(std::abs(stateDensityMatrix[i][j]) - 0.125) < - 0.000001); - } - } - const auto probVector = state.getSparseProbabilityVector(nrQubits, 0.001); - for (const auto& [s, prob] : probVector) { - constexpr double tolerance = 1e-10; - std::cout << s << ": " << prob << "\n"; - EXPECT_NEAR(prob, 0.125, tolerance); - } -} - -TEST(DDPackageTest, dNodeMulCache1) { - // Make caching test with dNodes - constexpr auto nrQubits = 1U; - const auto dd = std::make_unique( - nrQubits, DENSITY_MATRIX_SIMULATOR_DD_PACKAGE_CONFIG); - // Make zero density matrix - auto state = dd->makeZeroDensityOperator(nrQubits); - - const auto operation = getDD(qc::StandardOperation(0, qc::H), *dd); - dd->applyOperationToDensity(state, operation); - - state = dd->makeZeroDensityOperator(nrQubits); - auto& computeTable = dd->getMultiplicationComputeTable(); - - const auto& densityMatrix0 = - densityFromMatrixEdge(dd->conjugateTranspose(operation)); - - const auto* cachedResult = - computeTable.lookup(state.p, densityMatrix0.p, false); - ASSERT_NE(cachedResult, nullptr); - ASSERT_NE(cachedResult->p, nullptr); - state = dd->multiply(state, densityMatrix0, false); - ASSERT_NE(state.p, nullptr); - ASSERT_EQ(state.p, cachedResult->p); - - const auto densityMatrix1 = densityFromMatrixEdge(operation); - const auto* cachedResult1 = - computeTable.lookup(densityMatrix1.p, state.p, true); - ASSERT_NE(cachedResult1, nullptr); - ASSERT_NE(cachedResult1->p, nullptr); - const auto state2 = dd->multiply(densityMatrix1, state, true); - ASSERT_NE(state2.p, nullptr); - ASSERT_EQ(state2.p, cachedResult1->p); - - // try a repeated lookup - const auto* cachedResult2 = - computeTable.lookup(densityMatrix1.p, state.p, true); - ASSERT_NE(cachedResult2, nullptr); - ASSERT_NE(cachedResult2->p, nullptr); - ASSERT_EQ(cachedResult2->p, cachedResult1->p); - - computeTable.clear(); - const auto* cachedResult3 = - computeTable.lookup(densityMatrix1.p, state.p, true); - ASSERT_EQ(cachedResult3, nullptr); -} - -TEST(DDPackageTest, dNoiseCache) { - // Test the flags for dnode, vnode and mnodes - constexpr auto nrQubits = 1U; - const auto dd = std::make_unique(nrQubits); - // Make zero density matrix - const auto initialState = dd->makeZeroDensityOperator(nrQubits); - - // nothing pre-cached - const std::vector target = {0}; - const auto cachedNoise = dd->densityNoise.lookup(initialState, target); - ASSERT_EQ(cachedNoise.p, nullptr); - - auto state = initialState; - const auto operation = getDD(qc::StandardOperation(0, qc::X), *dd); - dd->applyOperationToDensity(state, operation); - dd->densityNoise.insert(initialState, state, target); - - // noise pre-cached - const auto cachedNoise1 = dd->densityNoise.lookup(initialState, target); - ASSERT_NE(cachedNoise1.p, nullptr); - ASSERT_EQ(cachedNoise1.p, state.p); - - // no noise pre-cached after clear - dd->densityNoise.clear(); - const auto cachedNoise2 = dd->densityNoise.lookup(initialState, target); - ASSERT_EQ(cachedNoise2.p, nullptr); -} - TEST(DDPackageTest, calCulpDistance) { constexpr auto nrQubits = 1U; auto dd = std::make_unique(nrQubits); @@ -1739,46 +1559,6 @@ TEST(DDPackageTest, calCulpDistance) { EXPECT_EQ(tmp1, 0); } -TEST(DDPackageTest, dStochCache) { - constexpr auto nrQubits = 4U; - const auto dd = std::make_unique( - nrQubits, STOCHASTIC_NOISE_SIMULATOR_DD_PACKAGE_CONFIG); - - std::vector operations = {}; - operations.emplace_back(getDD(qc::StandardOperation(0, qc::X), *dd)); - operations.emplace_back(getDD(qc::StandardOperation(1, qc::Z), *dd)); - operations.emplace_back(getDD(qc::StandardOperation(2, qc::Y), *dd)); - operations.emplace_back(getDD(qc::StandardOperation(3, qc::H), *dd)); - - dd->stochasticNoiseOperationCache.insert( - 0, 0, operations[0]); // insert X operations with target 0 - dd->stochasticNoiseOperationCache.insert( - 1, 1, operations[1]); // insert Z operations with target 1 - dd->stochasticNoiseOperationCache.insert( - 2, 2, operations[2]); // insert Y operations with target 2 - dd->stochasticNoiseOperationCache.insert( - 3, 3, operations[3]); // insert H operations with target 3 - - for (std::uint8_t i = 0; i < 4; i++) { - for (Qubit j = 0; j < 4; j++) { - const auto* op = dd->stochasticNoiseOperationCache.lookup(i, j); - if (static_cast(i) == j) { - EXPECT_TRUE(op != nullptr && op->p == operations[i].p); - } else { - EXPECT_EQ(op, nullptr); - } - } - } - - dd->stochasticNoiseOperationCache.clear(); - for (std::uint8_t i = 0; i < 4; i++) { - for (Qubit j = 0; j < 4; j++) { - auto* op = dd->stochasticNoiseOperationCache.lookup(i, j); - EXPECT_EQ(op, nullptr); - } - } -} - TEST(DDPackageTest, stateFromVectorBell) { const auto dd = std::make_unique(2); const auto v = std::vector>{SQRT2_2, 0, 0, SQRT2_2}; @@ -2566,14 +2346,10 @@ TEST(DDPackageTest, DataStructureStatistics) { EXPECT_EQ(stats["vNode"]["alignment_B"], 8U); EXPECT_EQ(stats["mNode"]["size_B"], 112U); EXPECT_EQ(stats["mNode"]["alignment_B"], 8U); - EXPECT_EQ(stats["dNode"]["size_B"], 112U); - EXPECT_EQ(stats["dNode"]["alignment_B"], 8U); EXPECT_EQ(stats["vEdge"]["size_B"], 24U); EXPECT_EQ(stats["vEdge"]["alignment_B"], 8U); EXPECT_EQ(stats["mEdge"]["size_B"], 24U); EXPECT_EQ(stats["mEdge"]["alignment_B"], 8U); - EXPECT_EQ(stats["dEdge"]["size_B"], 24U); - EXPECT_EQ(stats["dEdge"]["alignment_B"], 8U); EXPECT_EQ(stats["RealNumber"]["size_B"], 16U); EXPECT_EQ(stats["RealNumber"]["alignment_B"], 8U); } @@ -2588,7 +2364,6 @@ TEST(DDPackageTest, DDStatistics) { std::cout << stats.dump(2) << "\n"; EXPECT_TRUE(stats.contains("vector")); ASSERT_TRUE(stats.contains("matrix")); - EXPECT_TRUE(stats.contains("density_matrix")); EXPECT_TRUE(stats.contains("real_numbers")); EXPECT_TRUE(stats.contains("compute_tables")); const auto& matrixStats = stats["matrix"]; diff --git a/test/dd/test_state_generation.cpp b/test/dd/test_state_generation.cpp index e8f2120d4d..cd2d111b8c 100644 --- a/test/dd/test_state_generation.cpp +++ b/test/dd/test_state_generation.cpp @@ -397,7 +397,8 @@ TEST(StateGenerationTest, GenerateRandomOneQubit) { const std::vector nodesPerLevel{1}; const auto dd = std::make_unique(nq); - const auto state = generateRandomState(nq, nodesPerLevel, ROUNDROBIN, *dd); + const auto state = generateRandomState( + nq, nodesPerLevel, GenerationWireStrategy::ROUNDROBIN, *dd); EXPECT_NEAR(norm(state.getVector()), 1., 1e-6); EXPECT_EQ(state.size(), 2); // Node plus Terminal. @@ -425,7 +426,8 @@ TEST(StateGenerationTest, GenerateRandomRoundRobin) { std::accumulate(nodesPerLevel.begin(), nodesPerLevel.end(), 0UL); const auto dd = std::make_unique(nq); - const auto state = generateRandomState(nq, nodesPerLevel, ROUNDROBIN, *dd); + const auto state = generateRandomState( + nq, nodesPerLevel, GenerationWireStrategy::ROUNDROBIN, *dd); const auto rebuild = makeStateFromVector(state.getVector(), *dd); EXPECT_NEAR(norm(state.getVector()), 1., 1e-6); @@ -457,8 +459,8 @@ TEST(StateGenerationTest, GenerateRandomRoundRobinWithSeed) { std::accumulate(nodesPerLevel.begin(), nodesPerLevel.end(), 0UL); const auto dd = std::make_unique(nq); - const auto state = - generateRandomState(nq, nodesPerLevel, ROUNDROBIN, *dd, 72U); + const auto state = generateRandomState( + nq, nodesPerLevel, GenerationWireStrategy::ROUNDROBIN, *dd, 72U); const auto rebuild = makeStateFromVector(state.getVector(), *dd); EXPECT_NEAR(norm(state.getVector()), 1., 1e-6); @@ -489,7 +491,8 @@ TEST(StateGenerationTest, GenerateRandomRandom) { std::accumulate(nodesPerLevel.begin(), nodesPerLevel.end(), 0UL); const auto dd = std::make_unique(nq); - const auto state = generateRandomState(nq, nodesPerLevel, RANDOM, *dd); + const auto state = generateRandomState(nq, nodesPerLevel, + GenerationWireStrategy::RANDOM, *dd); const auto rebuild = makeStateFromVector(state.getVector(), *dd); EXPECT_NEAR(norm(state.getVector()), 1., 1e-6); @@ -521,7 +524,8 @@ TEST(StateGenerationTest, GenerateRandomRandomWithSeed) { std::accumulate(nodesPerLevel.begin(), nodesPerLevel.end(), 0UL); const auto dd = std::make_unique(nq); - const auto state = generateRandomState(nq, nodesPerLevel, RANDOM, *dd, 1337U); + const auto state = generateRandomState( + nq, nodesPerLevel, GenerationWireStrategy::RANDOM, *dd, 1337U); const auto rebuild = makeStateFromVector(state.getVector(), *dd); EXPECT_NEAR(norm(state.getVector()), 1., 1e-6); @@ -548,18 +552,29 @@ TEST(StateGenerationTest, GenerateRandomInvalidArguments) { auto dd = std::make_unique(nq); EXPECT_THROW( - { generateRandomState(nq + 1, {}, RANDOM, *dd, 1337U); }, + { + generateRandomState(nq + 1, {}, GenerationWireStrategy::RANDOM, *dd, + 1337U); + }, std::invalid_argument); EXPECT_THROW( - { generateRandomState(0, {0}, RANDOM, *dd, 1337U); }, + { + generateRandomState(0, {0}, GenerationWireStrategy::RANDOM, *dd, 1337U); + }, std::invalid_argument); EXPECT_THROW( - { generateRandomState(nq, {0}, RANDOM, *dd, 1337U); }, + { + generateRandomState(nq, {0}, GenerationWireStrategy::RANDOM, *dd, + 1337U); + }, std::invalid_argument); EXPECT_THROW( - { generateRandomState(nq, {1, 2, 5}, RANDOM, *dd, 1337U); }, + { + generateRandomState(nq, {1, 2, 5}, GenerationWireStrategy::RANDOM, *dd, + 1337U); + }, std::invalid_argument); } From a7b5606188da1d30f46c1c7ec1e4a64714382bf7 Mon Sep 17 00:00:00 2001 From: burgholzer Date: Sat, 17 Jan 2026 01:48:14 +0100 Subject: [PATCH 2/2] =?UTF-8?q?=F0=9F=91=8C=20Address=20CodeRabbit=20comme?= =?UTF-8?q?nts?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: burgholzer --- include/mqt-core/dd/DDDefinitions.hpp | 8 +++----- src/dd/Edge.cpp | 4 ++-- 2 files changed, 5 insertions(+), 7 deletions(-) diff --git a/include/mqt-core/dd/DDDefinitions.hpp b/include/mqt-core/dd/DDDefinitions.hpp index da936fa2d6..c1a7e76ca2 100644 --- a/include/mqt-core/dd/DDDefinitions.hpp +++ b/include/mqt-core/dd/DDDefinitions.hpp @@ -56,13 +56,11 @@ enum class BasisStates : std::uint8_t { left // NOLINT(readability-identifier-naming) }; -static constexpr fp SQRT2_2 = static_cast( +static constexpr auto SQRT2_2 = static_cast( 0.707106781186547524400844362104849039284835937688474036588L); static constexpr fp PI = std::numbers::pi; -static constexpr fp PI_2 = static_cast( - 1.570796326794896619231321691639751442098584699687552910487L); -static constexpr fp PI_4 = static_cast( - 0.785398163397448309615660845819875721049292349843776455243L); +static constexpr auto PI_2 = PI / 2; +static constexpr fp PI_4 = PI / 4; static constexpr std::uint64_t SERIALIZATION_VERSION = 1; diff --git a/src/dd/Edge.cpp b/src/dd/Edge.cpp index 15cf8159cb..952056035e 100644 --- a/src/dd/Edge.cpp +++ b/src/dd/Edge.cpp @@ -62,7 +62,7 @@ auto Edge::getValueByPath(const std::size_t numQubits, } // node is at the expected level - assert(tmp <= r.p->e.size()); + assert(tmp < r.p->e.size()); r = r.p->e[tmp]; c *= static_cast>(r.w); --level; @@ -449,7 +449,7 @@ auto Edge::printMatrix(const std::size_t numQubits) const -> void return; } // total number of qubits should not be lower than the highest qubit index - assert(isTerminal() || numQubits >= p->v); + assert(isTerminal() || numQubits > p->v); const std::size_t element = 1ULL << numQubits; for (auto i = 0ULL; i < element; ++i) { for (auto j = 0ULL; j < element; ++j) {