diff --git a/eval/eval_dd_package.cpp b/eval/eval_dd_package.cpp index 85f0f6aa0..6281caed2 100644 --- a/eval/eval_dd_package.cpp +++ b/eval/eval_dd_package.cpp @@ -77,7 +77,7 @@ benchmarkSimulate(const qc::QuantumComputation& qc) { const auto nq = qc.getNqubits(); exp->dd = std::make_unique(nq); const auto start = std::chrono::high_resolution_clock::now(); - const auto in = exp->dd->makeZeroState(nq); + const auto in = exp->dd->vectors().makeZeroState(nq); exp->sim = simulate(qc, in, *(exp->dd)); const auto end = std::chrono::high_resolution_clock::now(); exp->runtime = diff --git a/include/mqt-core/dd/ComputeTable.hpp b/include/mqt-core/dd/ComputeTable.hpp index e442baef1..b4722d639 100644 --- a/include/mqt-core/dd/ComputeTable.hpp +++ b/include/mqt-core/dd/ComputeTable.hpp @@ -35,11 +35,12 @@ namespace dd { template class ComputeTable { public: + static constexpr std::size_t DEFAULT_NUM_BUCKETS = 16384U; /** * Default constructor * @param numBuckets Number of hash table buckets. Must be a power of two. */ - explicit ComputeTable(const size_t numBuckets = 16384U) { + explicit ComputeTable(const size_t numBuckets = DEFAULT_NUM_BUCKETS) { // 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."); diff --git a/include/mqt-core/dd/DensityDDContainer.hpp b/include/mqt-core/dd/DensityDDContainer.hpp new file mode 100644 index 000000000..45a2d5f58 --- /dev/null +++ b/include/mqt-core/dd/DensityDDContainer.hpp @@ -0,0 +1,31 @@ +/* + * Copyright (c) 2025 Chair for Design Automation, TUM + * All rights reserved. + * + * SPDX-License-Identifier: MIT + * + * Licensed under the MIT License + */ + +#pragma once + +#include "dd/DesicionDiagramContainer.hpp" + +namespace dd { + +class DensityDDContainer : public DDContainer { +public: + using DDContainer::DDContainer; + + /** + * @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); + + char measureOneCollapsing(dEdge& e, Qubit index, std::mt19937_64& mt); +}; + +} // namespace dd diff --git a/include/mqt-core/dd/MatrixDDContainer.hpp b/include/mqt-core/dd/MatrixDDContainer.hpp new file mode 100644 index 000000000..0ac2a2597 --- /dev/null +++ b/include/mqt-core/dd/MatrixDDContainer.hpp @@ -0,0 +1,275 @@ +/* + * Copyright (c) 2025 Chair for Design Automation, TUM + * All rights reserved. + * + * SPDX-License-Identifier: MIT + * + * Licensed under the MIT License + */ + +#pragma once + +#include "dd/DesicionDiagramContainer.hpp" +#include "dd/UnaryComputeTable.hpp" + +namespace dd { + +class MatrixDDContainer : public DDContainer { +public: + struct Config : public DDContainer::Config { + std::size_t ctMatConjTransposeNumBucket = + UnaryComputeTable::DEFAULT_NUM_BUCKETS; + }; + + MatrixDDContainer(std::size_t nqubits, RealNumberUniqueTable& cUt, + ComplexNumbers& cn, const Config& config) + : DDContainer(nqubits, cUt, cn, config), + conjugateMatrixTranspose(config.ctMatConjTransposeNumBucket) {} + + inline void reset() { + DDContainer::reset(); + conjugateMatrixTranspose.clear(); + } + + inline bool garbageCollect(bool force) { + const bool collect = DDContainer::garbageCollect(force); + if (collect) { + conjugateMatrixTranspose.clear(); + } + return collect; + } + + /** + * @brief Construct the DD for a single-qubit gate + * @param mat The matrix representation of the gate + * @param target The target qubit + * @return A decision diagram for the gate + */ + [[nodiscard]] mEdge makeGateDD(const GateMatrix& mat, qc::Qubit target); + + /** + * @brief Construct the DD for a single-qubit controlled gate + * @param mat The matrix representation of the gate + * @param control The control qubit + * @param target The target qubit + * @return A decision diagram for the gate + */ + [[nodiscard]] mEdge makeGateDD(const GateMatrix& mat, + const qc::Control& control, qc::Qubit target); + + /** + * @brief Construct the DD for a multi-controlled single-qubit gate + * @param mat The matrix representation of the gate + * @param controls The control qubits + * @param target The target qubit + * @return A decision diagram for the gate + */ + [[nodiscard]] mEdge makeGateDD(const GateMatrix& mat, + const qc::Controls& controls, + qc::Qubit target); + + /** + * @brief Creates the DD for a two-qubit gate + * @param mat Matrix representation of the gate + * @param target0 First target qubit + * @param target1 Second target qubit + * @return DD representing the gate + * @throws std::runtime_error if the number of qubits is larger than the + * package configuration + */ + [[nodiscard]] mEdge makeTwoQubitGateDD(const TwoQubitGateMatrix& mat, + qc::Qubit target0, qc::Qubit target1); + + /** + * @brief Creates the DD for a two-qubit gate + * @param mat Matrix representation of the gate + * @param control Control qubit of the two-qubit gate + * @param target0 First target qubit + * @param target1 Second target qubit + * @return DD representing the gate + * @throws std::runtime_error if the number of qubits is larger than the + * package configuration + */ + [[nodiscard]] mEdge makeTwoQubitGateDD(const TwoQubitGateMatrix& mat, + const qc::Control& control, + qc::Qubit target0, qc::Qubit target1); + + /** + * @brief Creates the DD for a two-qubit gate + * @param mat Matrix representation of the gate + * @param controls Control qubits of the two-qubit gate + * @param target0 First target qubit + * @param target1 Second target qubit + * @return DD representing the gate + * @throws std::runtime_error if the number of qubits is larger than the + * package configuration + */ + [[nodiscard]] mEdge makeTwoQubitGateDD(const TwoQubitGateMatrix& mat, + const qc::Controls& controls, + qc::Qubit target0, qc::Qubit target1); + + /** + * @brief Converts a given matrix to a decision diagram + * @param matrix A complex matrix to convert to a DD. + * @return A decision diagram representing the matrix. + * @throws std::invalid_argument If the given matrix is not square or its + * length is not a power of two. + */ + [[nodiscard]] mEdge makeDDFromMatrix(const CMat& matrix); + +private: + /** + * @brief Constructs a decision diagram (DD) from a complex matrix using a + * recursive algorithm. + * + * @param matrix The complex matrix from which to create the DD. + * @param level The current level of recursion. Starts at the highest level of + * the matrix (log base 2 of the matrix size - 1). + * @param rowStart The starting row of the quadrant being processed. + * @param rowEnd The ending row of the quadrant being processed. + * @param colStart The starting column of the quadrant being processed. + * @param colEnd The ending column of the quadrant being processed. + * @return An mCachedEdge representing the root node of the created DD. + * + * @details This function recursively breaks down the matrix into quadrants + * until each quadrant has only one element. At each level of recursion, four + * new edges are created, one for each quadrant of the matrix. The four + * resulting decision diagram edges are used to create a new decision diagram + * node at the current level, and this node is returned as the result of the + * current recursive call. At the base case of recursion, the matrix has only + * one element, which is converted into a terminal node of the decision + * diagram. + * + * @note This function assumes that the matrix size is a power of two. + */ + [[nodiscard]] mCachedEdge makeDDFromMatrix(const CMat& matrix, Qubit level, + std::size_t rowStart, + std::size_t rowEnd, + std::size_t colStart, + std::size_t colEnd); + +public: + /** + * @brief Computes the conjugate transpose of a given matrix edge. + * + * @param a The matrix edge to conjugate transpose. + * @return The conjugated transposed matrix edge. + */ + [[nodiscard]] mEdge conjugateTranspose(const mEdge& a); + + /** + * @brief Recursively computes the conjugate transpose of a given matrix edge. + * + * @param a The matrix edge to conjugate transpose. + * @return The conjugated transposed matrix edge. + */ + [[nodiscard]] mCachedEdge conjugateTransposeRec(const mEdge& a); + + /** + * @brief Checks if a given matrix is close to the identity matrix. + * @details This function checks if a given matrix is close to the identity + * matrix, while ignoring any potential garbage qubits and ignoring the + * diagonal weights if `checkCloseToOne` is set to false. + * @param m An mEdge that represents the DD of the matrix. + * @param tol The accepted tolerance for the edge weights of the DD. + * @param garbage A vector of boolean values that defines which qubits are + * considered garbage qubits. If it's empty, then no qubit is considered to be + * a garbage qubit. + * @param checkCloseToOne If false, the function only checks if the matrix is + * close to a diagonal matrix. + */ + [[nodiscard]] bool isCloseToIdentity(const mEdge& m, fp tol = 1e-10, + const std::vector& garbage = {}, + bool checkCloseToOne = true) const; + + /** + * @brief Recursively checks if a given matrix is close to the identity + * matrix. + * + * @param m The matrix edge to check. + * @param visited A set of visited nodes to avoid redundant checks. + * @param tol The tolerance for comparing edge weights. + * @param garbage A vector of boolean values indicating which qubits are + * considered garbage. + * @param checkCloseToOne A flag to indicate whether to check if diagonal + * elements are close to one. + * @return True if the matrix is close to the identity matrix, false + * otherwise. + */ + static bool isCloseToIdentityRecursive( + const mEdge& m, std::unordered_set& visited, fp tol, + const std::vector& garbage, bool checkCloseToOne); + + /** + * @brief Reduces the decision diagram by handling ancillary qubits. + * + * @param e The matrix decision diagram edge to be reduced. + * @param ancillary A boolean vector indicating which qubits are ancillary + * (true) or not (false). + * @param regular Flag indicating whether to perform regular (true) or inverse + * (false) reduction. + * @return The reduced matrix decision diagram edge. + * + * @details This function modifies the decision diagram to account for + * ancillary qubits by: + * 1. Early returning if there are no ancillary qubits or if the edge is zero + * 2. Special handling for identity matrices by creating appropriate zero + * nodes + * 3. Finding the lowest ancillary qubit as a starting point + * 4. Recursively reducing nodes starting from the lowest ancillary qubit + * 5. Adding zero nodes for any remaining higher ancillary qubits + * + * The function maintains proper reference counting by incrementing the + * reference count of the result and decrementing the reference count of the + * input edge. + */ + mEdge reduceAncillae(mEdge e, const std::vector& ancillary, + bool regular = true); + + /// Create identity DD represented by the one-terminal. + [[nodiscard]] static mEdge makeIdent(); + + [[nodiscard]] mEdge createInitialMatrix(const std::vector& ancillary); + + /** + * @brief Reduces garbage qubits in a matrix decision diagram. + * + * @param e The matrix decision diagram edge to be reduced. + * @param garbage A boolean vector indicating which qubits are garbage (true) + * or not (false). + * @param regular Flag indicating whether to apply regular (true) or inverse + * (false) reduction. In regular mode, garbage entries are summed in the first + * two components, in inverse mode, they are summed in the first and third + * components. + * @param normalizeWeights Flag indicating whether to normalize weights to + * their magnitudes. When true, all weights in the DD are changed to their + * magnitude, also for non-garbage qubits. This is used for checking partial + * equivalence where only measurement probabilities matter. + * @return The reduced matrix decision diagram edge. + * + * @details For each garbage qubit q, this function sums all the entries for + * q=0 and q=1, setting the entry for q=0 to the sum and the entry for q=1 to + * zero. To maintain proper probabilities, the function computes sqrt(|a|^2 + + * |b|^2) for two entries a and b. The function handles special cases like + * zero terminals and identity matrices separately and maintains proper + * reference counting throughout the reduction process. + */ + [[nodiscard]] mEdge reduceGarbage(const mEdge& e, + const std::vector& garbage, + bool regular = true, + bool normalizeWeights = false); + +private: + [[nodiscard]] mCachedEdge + reduceAncillaeRecursion(mNode* p, const std::vector& ancillary, + Qubit lowerbound, bool regular = true); + + [[nodiscard]] mCachedEdge + reduceGarbageRecursion(mNode* p, const std::vector& garbage, + Qubit lowerbound, bool regular = true, + bool normalizeWeights = false); + + UnaryComputeTable conjugateMatrixTranspose; +}; + +} // namespace dd diff --git a/include/mqt-core/dd/MemoryManager.hpp b/include/mqt-core/dd/MemoryManager.hpp index 97ac1919e..c49cf1f77 100644 --- a/include/mqt-core/dd/MemoryManager.hpp +++ b/include/mqt-core/dd/MemoryManager.hpp @@ -39,9 +39,14 @@ struct LLBase; * edge weights, etc. */ class MemoryManager { - MemoryManager(size_t entrySize, std::size_t initialAllocationSize); - public: + struct Config { + std::size_t initialAllocationSize = INITIAL_ALLOCATION_SIZE; + std::size_t entrySize; + }; + + MemoryManager(const Config& config); + // delete copy construction and assignment MemoryManager(const MemoryManager&) = delete; MemoryManager& operator=(const MemoryManager&) = delete; @@ -71,7 +76,7 @@ class MemoryManager { template static MemoryManager create(const std::size_t initialAllocationSize = INITIAL_ALLOCATION_SIZE) { - return {sizeof(T), initialAllocationSize}; + return {{sizeof(T), initialAllocationSize}}; } /// default destructor diff --git a/include/mqt-core/dd/Node.hpp b/include/mqt-core/dd/Node.hpp index fd3779a2c..fd6848d1f 100644 --- a/include/mqt-core/dd/Node.hpp +++ b/include/mqt-core/dd/Node.hpp @@ -84,7 +84,6 @@ struct vNode final : NodeBase { // NOLINT(readability-identifier-naming) }; using vEdge = Edge; using vCachedEdge = CachedEdge; -using VectorDD = vEdge; /** * @brief A matrix DD node @@ -103,7 +102,6 @@ struct mNode final : NodeBase { // NOLINT(readability-identifier-naming) }; using mEdge = Edge; using mCachedEdge = CachedEdge; -using MatrixDD = mEdge; /** * @brief A density matrix DD node diff --git a/include/mqt-core/dd/NoiseFunctionality.hpp b/include/mqt-core/dd/NoiseFunctionality.hpp index 7ea3180ce..c441e0e14 100644 --- a/include/mqt-core/dd/NoiseFunctionality.hpp +++ b/include/mqt-core/dd/NoiseFunctionality.hpp @@ -53,7 +53,7 @@ class StochasticNoiseFunctionality { double multiQubitGateFactor, const std::string& cNoiseEffects); - ~StochasticNoiseFunctionality() { package->decRef(identityDD); } + ~StochasticNoiseFunctionality() { package->matrices().decRef(identityDD); } protected: Package* package; diff --git a/include/mqt-core/dd/Operations.hpp b/include/mqt-core/dd/Operations.hpp index e18e60229..a0f8ffb99 100644 --- a/include/mqt-core/dd/Operations.hpp +++ b/include/mqt-core/dd/Operations.hpp @@ -254,8 +254,8 @@ void changePermutation(DDType& on, qc::Permutation& from, // swap i and j auto saved = on; - const auto swapDD = dd.makeTwoQubitGateDD(opToTwoQubitGateMatrix(qc::SWAP), - from.at(i), from.at(j)); + const auto swapDD = dd.matrices().makeTwoQubitGateDD( + opToTwoQubitGateMatrix(qc::SWAP), from.at(i), from.at(j)); if constexpr (std::is_same_v) { on = dd.multiply(swapDD, on); } else { diff --git a/include/mqt-core/dd/Package.hpp b/include/mqt-core/dd/Package.hpp index c4ac637be..957010755 100644 --- a/include/mqt-core/dd/Package.hpp +++ b/include/mqt-core/dd/Package.hpp @@ -16,8 +16,10 @@ #include "dd/ComputeTable.hpp" #include "dd/DDDefinitions.hpp" #include "dd/DDpackageConfig.hpp" +#include "dd/DensityDDContainer.hpp" #include "dd/DensityNoiseTable.hpp" #include "dd/Edge.hpp" +#include "dd/MatrixDDContainer.hpp" #include "dd/MemoryManager.hpp" #include "dd/Node.hpp" #include "dd/Package_fwd.hpp" // IWYU pragma: export @@ -26,6 +28,7 @@ #include "dd/StochasticNoiseOperationTable.hpp" #include "dd/UnaryComputeTable.hpp" #include "dd/UniqueTable.hpp" +#include "dd/VectorDDContainer.hpp" #include "ir/Definitions.hpp" #include "ir/Permutation.hpp" #include "ir/operations/Control.hpp" @@ -107,876 +110,151 @@ class Package { /// Get the number of qubits [[nodiscard]] auto qubits() const { return nqubits; } -private: +public: std::size_t nqubits; DDPackageConfig config_; -public: - /// The memory manager for vector nodes - MemoryManager vMemoryManager{ - MemoryManager::create(config_.utVecInitialAllocationSize)}; - /// 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 - * separately. Hence, it suffices for the manager to only manage real numbers. - */ - MemoryManager cMemoryManager{MemoryManager::create()}; - - /** - * @brief Get the memory manager for a given type - * @tparam T The type to get the manager for - * @return A reference to the manager - */ - template [[nodiscard]] auto& getMemoryManager() { - if constexpr (std::is_same_v) { - 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; - } - } - - /** - * @brief Reset all memory managers - * @arg resizeToTotal If set to true, each manager allocates one chunk of - * memory as large as all chunks combined before the reset. - * @see MemoryManager::reset - */ - void resetMemoryManagers(bool resizeToTotal = false); - - /// The unique table used for vector nodes - UniqueTable vUniqueTable{vMemoryManager, {0U, 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}}; - /** - * @brief The unique table used for complex numbers - * @note The table actually only stores real numbers in the interval [0, 1], - * but is used to manages all complex numbers throughout the package. - * @see RealNumberUniqueTable - */ - RealNumberUniqueTable cUniqueTable{cMemoryManager}; - ComplexNumbers cn{cUniqueTable}; - - /** - * @brief Get the unique table for a given type - * @tparam T The type to get the unique table for - * @return A reference to the unique table - */ - template [[nodiscard]] auto& getUniqueTable() { - if constexpr (std::is_same_v) { - 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; - } - } - - /** - * @brief Clear all unique tables - * @see UniqueTable::clear - * @see RealNumberUniqueTable::clear - */ - void clearUniqueTables(); - - /** - * @brief Increment the reference count of an edge - * @details This is the main function for increasing reference counts within - * the DD package. It increases the reference count of the complex edge weight - * as well as the DD node itself. If the node newly becomes active, meaning - * that it had a reference count of zero beforehand, the reference count of - * all children is recursively increased. - * @tparam Node The node type of the edge. - * @param e The edge to increase the reference count of - */ - template void incRef(const Edge& e) noexcept { - cn.incRef(e.w); - const auto& p = e.p; - const auto inc = getUniqueTable().incRef(p); - if (inc && p->ref == 1U) { - for (const auto& child : p->e) { - incRef(child); - } - } - } - - /** - * @brief Decrement the reference count of an edge - * @details This is the main function for decreasing reference counts within - * the DD package. It decreases the reference count of the complex edge weight - * as well as the DD node itself. If the node newly becomes dead, meaning - * that its reference count hit zero, the reference count of all children is - * recursively decreased. - * @tparam Node The node type of the edge. - * @param e The edge to decrease the reference count of - */ - template void decRef(const Edge& e) noexcept { - cn.decRef(e.w); - const auto& p = e.p; - const auto dec = getUniqueTable().decRef(p); - if (dec && p->ref == 0U) { - for (const auto& child : p->e) { - decRef(child); - } - } - } - - /** - * @brief Trigger garbage collection in all unique tables - * - * @details Garbage collection is the process of removing all nodes from the - * unique tables that have a reference count of zero. - * Such nodes are considered "dead" and they can be safely removed from the - * unique tables. This process is necessary to free up memory that is no - * longer needed. By default, garbage collection is only triggered if the - * unique table indicates that it possibly needs collection. Whenever some - * nodes are recollected, some compute tables need to be invalidated as well. - * - * @param force - * @return - */ - bool garbageCollect(bool force = false); - - /// - /// 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); - - /** - * @brief Construct the all-zero state \f$|0...0\rangle\f$ - * @param n The number of qubits - * @param start The starting qubit index. Default is 0. - * @return A decision diagram for the all-zero state - */ - vEdge makeZeroState(std::size_t n, std::size_t start = 0); - - /** - * @brief Construct a computational basis state \f$|b_{n-1}...b_0\rangle\f$ - * @param n The number of qubits - * @param state The state to construct - * @param start The starting qubit index. Default is 0. - * @return A decision diagram for the computational basis state - */ - vEdge makeBasisState(std::size_t n, const std::vector& state, - std::size_t start = 0); - - /** - * @brief Construct a product state out of - * \f$\{0, 1, +, -, R, L\}^{\otimes n}\f$. - * @param n The number of qubits - * @param state The state to construct - * @param start The starting qubit index. Default is 0. - * @return A decision diagram for the product state - */ - vEdge makeBasisState(std::size_t n, const std::vector& state, - std::size_t start = 0); - - /** - * @brief Construct a GHZ state \f$|0...0\rangle + |1...1\rangle\f$ - * @param n The number of qubits - * @return A decision diagram for the GHZ state - */ - vEdge makeGHZState(std::size_t n); - - /** - * @brief Construct a W state - * @details The W state is defined as - * \f[ - * |0...01\rangle + |0...10\rangle + |10...0\rangle - * \f] - * @param n The number of qubits - * @return A decision diagram for the W state - */ - vEdge makeWState(std::size_t n); - - /** - * @brief Construct a decision diagram from an arbitrary state vector - * @param stateVector The state vector to convert to a DD - * @return A decision diagram for the state - */ - vEdge makeStateFromVector(const CVec& stateVector); - - /// - /// Matrix nodes, edges and quantum gates - /// - - /** - * @brief Construct the DD for a single-qubit gate - * @param mat The matrix representation of the gate - * @param target The target qubit - * @return A decision diagram for the gate - */ - mEdge makeGateDD(const GateMatrix& mat, qc::Qubit target); - - /** - * @brief Construct the DD for a single-qubit controlled gate - * @param mat The matrix representation of the gate - * @param control The control qubit - * @param target The target qubit - * @return A decision diagram for the gate - */ - mEdge makeGateDD(const GateMatrix& mat, const qc::Control& control, - qc::Qubit target); - - /** - * @brief Construct the DD for a multi-controlled single-qubit gate - * @param mat The matrix representation of the gate - * @param controls The control qubits - * @param target The target qubit - * @return A decision diagram for the gate - */ - mEdge makeGateDD(const GateMatrix& mat, const qc::Controls& controls, - qc::Qubit target); - - /** - * @brief Creates the DD for a two-qubit gate - * @param mat Matrix representation of the gate - * @param target0 First target qubit - * @param target1 Second target qubit - * @return DD representing the gate - * @throws std::runtime_error if the number of qubits is larger than the - * package configuration - */ - mEdge makeTwoQubitGateDD(const TwoQubitGateMatrix& mat, qc::Qubit target0, - qc::Qubit target1); - - /** - * @brief Creates the DD for a two-qubit gate - * @param mat Matrix representation of the gate - * @param control Control qubit of the two-qubit gate - * @param target0 First target qubit - * @param target1 Second target qubit - * @return DD representing the gate - * @throws std::runtime_error if the number of qubits is larger than the - * package configuration - */ - mEdge makeTwoQubitGateDD(const TwoQubitGateMatrix& mat, - const qc::Control& control, qc::Qubit target0, - qc::Qubit target1); - - /** - * @brief Creates the DD for a two-qubit gate - * @param mat Matrix representation of the gate - * @param controls Control qubits of the two-qubit gate - * @param target0 First target qubit - * @param target1 Second target qubit - * @return DD representing the gate - * @throws std::runtime_error if the number of qubits is larger than the - * package configuration - */ - mEdge makeTwoQubitGateDD(const TwoQubitGateMatrix& mat, - const qc::Controls& controls, qc::Qubit target0, - qc::Qubit target1); - - /** - * @brief Converts a given matrix to a decision diagram - * @param matrix A complex matrix to convert to a DD. - * @return A decision diagram representing the matrix. - * @throws std::invalid_argument If the given matrix is not square or its - * length is not a power of two. - */ - mEdge makeDDFromMatrix(const CMat& matrix); - -private: - /** - * @brief Constructs a decision diagram (DD) from a state vector using a - * recursive algorithm. - * - * @param begin Iterator pointing to the beginning of the state vector. - * @param end Iterator pointing to the end of the state vector. - * @param level The current level of recursion. Starts at the highest level of - * the state vector (log base 2 of the vector size - 1). - * @return A vCachedEdge representing the root node of the created DD. - * - * @details This function recursively breaks down the state vector into halves - * until each half has only one element. At each level of recursion, two new - * edges are created, one for each half of the state vector. The two resulting - * decision diagram edges are used to create a new decision diagram node at - * the current level, and this node is returned as the result of the current - * recursive call. At the base case of recursion, the state vector has only - * two elements, which are converted into terminal nodes of the decision - * diagram. - * - * @note This function assumes that the state vector size is a power of two. - */ - vCachedEdge makeStateFromVector(const CVec::const_iterator& begin, - const CVec::const_iterator& end, Qubit level); - - /** - * @brief Constructs a decision diagram (DD) from a complex matrix using a - * recursive algorithm. - * - * @param matrix The complex matrix from which to create the DD. - * @param level The current level of recursion. Starts at the highest level of - * the matrix (log base 2 of the matrix size - 1). - * @param rowStart The starting row of the quadrant being processed. - * @param rowEnd The ending row of the quadrant being processed. - * @param colStart The starting column of the quadrant being processed. - * @param colEnd The ending column of the quadrant being processed. - * @return An mCachedEdge representing the root node of the created DD. - * - * @details This function recursively breaks down the matrix into quadrants - * until each quadrant has only one element. At each level of recursion, four - * new edges are created, one for each quadrant of the matrix. The four - * resulting decision diagram edges are used to create a new decision diagram - * node at the current level, and this node is returned as the result of the - * current recursive call. At the base case of recursion, the matrix has only - * one element, which is converted into a terminal node of the decision - * diagram. - * - * @note This function assumes that the matrix size is a power of two. - */ - mCachedEdge makeDDFromMatrix(const CMat& matrix, Qubit level, - std::size_t rowStart, std::size_t rowEnd, - std::size_t colStart, std::size_t colEnd); - -public: - /** - * @brief Create a normalized DD node and return an edge pointing to it. - * - * @details The node is not recreated if it already exists. This function - * retrieves a node from the memory manager, sets its variable, and normalizes - * the edges. If the node resembles the identity, it is skipped. The function - * then looks up the node in the unique table and returns an edge pointing to - * it. - * - * @tparam Node The type of the node. - * @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) { - auto& memoryManager = getMemoryManager(); - auto p = memoryManager.template get(); - assert(p->ref == 0U); - - p->v = var; - if constexpr (std::is_same_v || std::is_same_v) { - 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 (!e.isTerminal()) { - const auto& es = e.p->e; - // Check if node resembles the identity. If so, skip it. - if ((es[0].p == es[3].p) && - (es[0].w.exactlyOne() && es[1].w.exactlyZero() && - es[2].w.exactlyZero() && es[3].w.exactlyOne())) { - auto* ptr = es[0].p; - memoryManager.returnEntry(*e.p); - return EdgeType{ptr, e.w}; - } - } - } - - // look it up in the unique tables - auto& uniqueTable = getUniqueTable(); - auto* l = uniqueTable.lookup(e.p); - - return EdgeType{l, e.w}; - } - - /** - * @brief Delete an edge from the decision diagram. - * - * @tparam Node The type of the node. - * @param e The edge to delete. - * @param v The variable associated with the edge. - * @param edgeIdx The index of the edge to delete. - * @return The modified edge after deletion. - */ - template - Edge deleteEdge(const Edge& e, const Qubit v, - const std::size_t edgeIdx) { - std::unordered_map> nodes{}; - return deleteEdge(e, v, edgeIdx, nodes); - } - - /** - * @brief Helper function to delete an edge from the decision diagram. - * - * @tparam Node The type of the node. - * @param e The edge to delete. - * @param v The variable associated with the edge. - * @param edgeIdx The index of the edge to delete. - * @param nodes A map to keep track of processed nodes. - * @return The modified edge after deletion. - */ - template - Edge deleteEdge(const Edge& e, const Qubit v, - const std::size_t edgeIdx, - std::unordered_map>& nodes) { - if (e.isTerminal()) { - return e; - } - - const auto& nodeIt = nodes.find(e.p); - Edge r{}; - if (nodeIt != nodes.end()) { - r = nodeIt->second; - } else { - constexpr std::size_t n = std::tuple_size_ve)>; - std::array, n> edges{}; - if (e.p->v == v) { - for (std::size_t i = 0; i < n; i++) { - edges[i] = i == edgeIdx - ? Edge::zero() - : e.p->e[i]; // optimization -> node cannot occur below - // again, since dd is assumed to be free - } - } else { - for (std::size_t i = 0; i < n; i++) { - edges[i] = deleteEdge(e.p->e[i], v, edgeIdx, nodes); - } - } - - r = makeDDNode(e.p->v, edges); - nodes[e.p] = r; - } - r.w = cn.lookup(r.w * e.w); - return r; - } - - /// - /// Compute table definitions - /// - - /** - * @brief Clear all compute tables. - * - * @details This method clears all entries in the compute tables used for - * various operations. It resets the state of the compute tables, making them - * ready for new computations. - */ - void clearComputeTables(); - - /// - /// Measurements from state decision diagrams - /// - - /** - * @brief Measure all qubits in the given decision diagram. - * - * @details This function measures all qubits in the decision diagram - * represented by `rootEdge`. It checks for numerical instabilities and - * collapses the state if requested. - * - * @param rootEdge The decision diagram to measure. - * @param collapse If true, the state is collapsed after measurement. - * @param mt A random number generator. - * @param epsilon The tolerance for numerical instabilities. - * @return A string representing the measurement result. - * @throws std::runtime_error If numerical instabilities are detected or if - * probabilities do not sum to 1. - */ - std::string measureAll(vEdge& rootEdge, bool collapse, std::mt19937_64& mt, - fp epsilon = 0.001); + MemoryManager cMm{MemoryManager::create()}; + RealNumberUniqueTable cUt{cMm}; + ComplexNumbers cn{cUt}; private: - /** - * @brief Assigns probabilities to nodes in a decision diagram. - * - * @details This function recursively assigns probabilities to nodes in a - * decision diagram. It calculates the probability of reaching each node and - * stores the result in a map. - * - * @param edge The edge to start the probability assignment from. - * @param probs A map to store the probabilities of each node. - * @return The probability of the given edge. - */ - static fp assignProbabilities(const vEdge& edge, - std::unordered_map& probs); + VectorDDContainer vectors_{nqubits, cUt, cn, {/*todo: add config here */}}; + MatrixDDContainer matrices_{nqubits, cUt, cn, {/*todo: add config here */}}; + DensityDDContainer densities_{nqubits, cUt, cn, {/*todo: add config here */}}; public: - /** - * @brief Determine the measurement probabilities for a given qubit index. - * - * @param rootEdge The root edge of the decision diagram. - * @param index The qubit index to determine the measurement probabilities - * for. - * @return A pair of floating-point values representing the probabilities of - * measuring 0 and 1, respectively. - * - * @details This function calculates the probabilities of measuring 0 and 1 - * for a given qubit index in the decision diagram. It uses a breadth-first - * search to traverse the decision diagram and accumulate the measurement - * probabilities. The function maintains a map of measurement probabilities - * for each node and a set of visited nodes to avoid redundant calculations. - * It also uses a queue to process nodes level by level. - */ - static std::pair - determineMeasurementProbabilities(const vEdge& rootEdge, Qubit index); + VectorDDContainer& vectors() { return vectors_; } + const VectorDDContainer& vectors() const { return vectors_; } - /** - * @brief Measures the qubit with the given index in the given state vector - * decision diagram. Collapses the state according to the measurement result. - * @param rootEdge the root edge of the state vector decision diagram - * @param index the index of the qubit to be measured - * @param mt the random number generator - * @param epsilon the numerical precision used for checking the normalization - * of the state vector decision diagram - * @return the measurement result ('0' or '1') - * @throws std::runtime_error if a numerical instability is detected during - * the measurement. - */ - char measureOneCollapsing(vEdge& rootEdge, Qubit index, std::mt19937_64& mt, - fp epsilon = 0.001); + MatrixDDContainer& matrices() { return matrices_; } + const MatrixDDContainer& matrices() const { return matrices_; } - char measureOneCollapsing(dEdge& e, Qubit index, std::mt19937_64& mt); + DensityDDContainer& densities() { return densities_; } + const DensityDDContainer& densities() const { return densities_; } - /** - * @brief Performs a specific measurement on the given state vector decision - * diagram. Collapses the state according to the measurement result. - * @param rootEdge the root edge of the state vector decision diagram - * @param index the index of the qubit to be measured - * @param probability the probability of the measurement result (required for - * normalization) - * @param measureZero whether or not to measure '0' (otherwise '1' is - * measured) - */ - void performCollapsingMeasurement(vEdge& rootEdge, Qubit index, - fp probability, bool measureZero); + void incRef(const vEdge& e) { vectors().incRef(e); } + void incRef(const mEdge& e) { matrices().incRef(e); } + void incRef(const dEdge& e) { densities().incRef(e); } - /// - /// Addition - /// - ComputeTable vectorAdd{ - config_.ctVecAddNumBucket}; - ComputeTable matrixAdd{ - config_.ctMatAddNumBucket}; - ComputeTable densityAdd{ - config_.ctDmAddNumBucket}; + void decRef(const vEdge& e) { vectors().decRef(e); } + void decRef(const mEdge& e) { matrices().decRef(e); } + void decRef(const dEdge& e) { densities().decRef(e); } - /** - * @brief Get the compute table for addition operations. - * - * @tparam Node The type of the node. - * @return A reference to the appropriate compute table for the given node - * type. - */ - template [[nodiscard]] auto& getAddComputeTable() { - if constexpr (std::is_same_v) { - return vectorAdd; - } else if constexpr (std::is_same_v) { - return matrixAdd; - } else if constexpr (std::is_same_v) { - return densityAdd; - } + [[nodiscard]] vEdge reduceGarbage(vEdge& e, const std::vector& garbage, + bool normalizeWeights = false) { + return vectors().reduceGarbage(e, garbage, normalizeWeights); } - ComputeTable vectorAddMagnitudes{ - config_.ctVecAddMagNumBucket}; - ComputeTable matrixAddMagnitudes{ - config_.ctMatAddMagNumBucket}; - - /** - * @brief Get the compute table for addition operations with magnitudes. - * - * @tparam Node The type of the node. - * @return A reference to the appropriate compute table for the given node - * type. - */ - template [[nodiscard]] auto& getAddMagnitudesComputeTable() { - if constexpr (std::is_same_v) { - return vectorAddMagnitudes; - } else if constexpr (std::is_same_v) { - return matrixAddMagnitudes; - } + [[nodiscard]] mEdge reduceGarbage(const mEdge& e, + const std::vector& garbage, + bool regular = true, + bool normalizeWeights = false) { + return matrices().reduceGarbage(e, garbage, regular, normalizeWeights); } - /** - * @brief Add two decision diagrams. - * - * @tparam Node The type of the node. - * @param x The first DD. - * @param y The second DD. - * @return The resulting DD after addition. - * - * @details This function performs the addition of two decision diagrams - * (DDs). It uses a compute table to cache intermediate results and avoid - * redundant computations. The addition is conducted recursively, where the - * function traverses the nodes of the DDs, adds corresponding edges, and - * normalizes the resulting edges. If the nodes are terminal, their weights - * are directly added. The function ensures that the resulting DD is properly - * normalized and stored in the unique table to maintain the canonical form. - */ - template - Edge add(const Edge& x, const Edge& y) { - Qubit var{}; - if (!x.isTerminal()) { - var = x.p->v; - } - if (!y.isTerminal() && (y.p->v) > var) { - var = y.p->v; - } - - const auto result = add2(CachedEdge{x.p, x.w}, {y.p, y.w}, var); - return cn.lookup(result); - } - - /** - * @brief Internal function to add two decision diagrams. - * - * This function is used internally to add two decision diagrams (DDs) of type - * Node. It is not intended to be called directly. - * - * @tparam Node The type of the node. - * @param x The first DD. - * @param y The second DD. - * @param var The variable associated with the current level of recursion. - * @return The resulting DD after addition. - */ - template - CachedEdge add2(const CachedEdge& x, const CachedEdge& y, - const Qubit var) { - if (x.w.exactlyZero()) { - if (y.w.exactlyZero()) { - return CachedEdge::zero(); - } - return y; - } - if (y.w.exactlyZero()) { - return x; - } - if (x.p == y.p) { - const auto rWeight = x.w + y.w; - return {x.p, rWeight}; - } - - auto& computeTable = getAddComputeTable(); - if (const auto* r = computeTable.lookup(x, y); r != nullptr) { - return *r; - } - - constexpr std::size_t n = std::tuple_size_ve)>; - 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 (x.isIdentity() || x.p->v < var) { - // [ 0 | 1 ] [ x | 0 ] - // --------- = --------- - // [ 2 | 3 ] [ 0 | x ] - if (i == 0 || i == 3) { - e1 = x; - } - } else { - auto& xSuccessor = x.p->e[i]; - e1 = {xSuccessor.p, 0}; - if (!xSuccessor.w.exactlyZero()) { - e1.w = x.w * xSuccessor.w; - } - } - } else { - auto& xSuccessor = x.p->e[i]; - e1 = {xSuccessor.p, 0}; - if (!xSuccessor.w.exactlyZero()) { - e1.w = x.w * xSuccessor.w; - } - } - CachedEdge e2{}; - if constexpr (std::is_same_v || - std::is_same_v) { - if (y.isIdentity() || y.p->v < var) { - // [ 0 | 1 ] [ y | 0 ] - // --------- = --------- - // [ 2 | 3 ] [ 0 | y ] - if (i == 0 || i == 3) { - e2 = y; - } - } else { - auto& ySuccessor = y.p->e[i]; - e2 = {ySuccessor.p, 0}; - if (!ySuccessor.w.exactlyZero()) { - e2.w = y.w * ySuccessor.w; - } - } - } else { - auto& ySuccessor = y.p->e[i]; - e2 = {ySuccessor.p, 0}; - if (!ySuccessor.w.exactlyZero()) { - 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); - } - } - auto r = makeDDNode(var, edge); - computeTable.insert(x, y, r); - return r; + vEdge makeDDNode( + const Qubit var, + const std::array>& edges, + [[maybe_unused]] const bool generateDensityMatrix = false) { + return vectors().makeDDNode(var, edges, generateDensityMatrix); } - /** - * @brief Compute the element-wise magnitude sum of two vectors or matrices. - * - * For two vectors (or matrices) \p x and \p y, this function returns a result - * \p r such that for each index \p i: - * \f$ r[i] = \sqrt{|x[i]|^2 + |y[i]|^2} \f$ - * - * @param x DD representation of the first operand. - * @param y DD representation of the second operand. - * @param var Number of qubits in the DD. - * @return DD representing the result. - */ - template - CachedEdge addMagnitudes(const CachedEdge& x, - const CachedEdge& y, const Qubit var) { - if (x.w.exactlyZero()) { - if (y.w.exactlyZero()) { - return CachedEdge::zero(); - } - const auto rWeight = y.w.mag(); - return {y.p, rWeight}; - } - if (y.w.exactlyZero()) { - const auto rWeight = x.w.mag(); - return {x.p, rWeight}; - } - if (x.p == y.p) { - const auto rWeight = std::sqrt(x.w.mag2() + y.w.mag2()); - return {x.p, rWeight}; - } - - auto& computeTable = getAddMagnitudesComputeTable(); - if (const auto* r = computeTable.lookup(x, y); r != nullptr) { - return *r; - } - - constexpr std::size_t n = std::tuple_size_ve)>; - 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 (x.isIdentity() || x.p->v < var) { - if (i == 0 || i == 3) { - e1 = x; - } - } else { - auto& xSuccessor = x.p->e[i]; - e1 = {xSuccessor.p, 0}; - if (!xSuccessor.w.exactlyZero()) { - e1.w = x.w * xSuccessor.w; - } - } - } else { - auto& xSuccessor = x.p->e[i]; - e1 = {xSuccessor.p, 0}; - if (!xSuccessor.w.exactlyZero()) { - e1.w = x.w * xSuccessor.w; - } - } - CachedEdge e2{}; - if constexpr (std::is_same_v || - std::is_same_v) { - if (y.isIdentity() || y.p->v < var) { - if (i == 0 || i == 3) { - e2 = y; - } - } else { - auto& ySuccessor = y.p->e[i]; - e2 = {ySuccessor.p, 0}; - if (!ySuccessor.w.exactlyZero()) { - e2.w = y.w * ySuccessor.w; - } - } - } else { - auto& ySuccessor = y.p->e[i]; - e2 = {ySuccessor.p, 0}; - if (!ySuccessor.w.exactlyZero()) { - e2.w = y.w * ySuccessor.w; - } - } - edge[i] = addMagnitudes(e1, e2, var - 1); - } - auto r = makeDDNode(var, edge); - computeTable.insert(x, y, r); - return r; + vCachedEdge + makeDDNode(const Qubit var, + const std::array>& edges, + [[maybe_unused]] const bool generateDensityMatrix = false) { + return vectors().makeDDNode(var, edges, generateDensityMatrix); } - /// - /// Vector conjugation - /// - UnaryComputeTable conjugateVector{ - config_.ctVecConjNumBucket}; + mEdge makeDDNode( + const Qubit var, + const std::array>& edges, + [[maybe_unused]] const bool generateDensityMatrix = false) { + return matrices().makeDDNode(var, edges, generateDensityMatrix); + } + mCachedEdge + makeDDNode(const Qubit var, + const std::array>& edges, + [[maybe_unused]] const bool generateDensityMatrix = false) { + return matrices().makeDDNode(var, edges, generateDensityMatrix); + } + + dEdge makeDDNode( + const Qubit var, + const std::array>& edges, + [[maybe_unused]] const bool generateDensityMatrix = false) { + return densities().makeDDNode(var, edges, generateDensityMatrix); + } + dCachedEdge + makeDDNode(const Qubit var, + const std::array>& edges, + [[maybe_unused]] const bool generateDensityMatrix = false) { + return densities().makeDDNode(var, edges, generateDensityMatrix); + } + + vCachedEdge add(const vCachedEdge& x, const vCachedEdge& y, Qubit var) { + return vectors().add(x, y, var); + } + mCachedEdge add(const mCachedEdge& x, const mCachedEdge& y, Qubit var) { + return matrices().add(x, y, var); + } + + dCachedEdge add(const dCachedEdge& x, const dCachedEdge& y, Qubit var) { + return densities().add(x, y, var); + } +public: /** - * @brief Conjugates a given decision diagram edge. + * @brief Trigger garbage collection in all unique tables * - * @param a The decision diagram edge to conjugate. - * @return The conjugated decision diagram edge. - */ - vEdge conjugate(const vEdge& a); - /** - * @brief Recursively conjugates a given decision diagram edge. + * @details Garbage collection is the process of removing all nodes from the + * unique tables that have a reference count of zero. + * Such nodes are considered "dead" and they can be safely removed from the + * unique tables. This process is necessary to free up memory that is no + * longer needed. By default, garbage collection is only triggered if the + * unique table indicates that it possibly needs collection. Whenever some + * nodes are recollected, some compute tables need to be invalidated as well. * - * @param a The decision diagram edge to conjugate. - * @return The conjugated decision diagram edge. + * @param force + * @return */ - vCachedEdge conjugateRec(const vEdge& a); - - /// - /// Matrix (conjugate) transpose - /// - UnaryComputeTable conjugateMatrixTranspose{ - config_.ctMatConjTransNumBucket}; + bool garbageCollect(bool force = false); /** - * @brief Computes the conjugate transpose of a given matrix edge. - * - * @param a The matrix edge to conjugate transpose. - * @return The conjugated transposed matrix edge. + * @brief Performs a specific measurement on the given state vector decision + * diagram. Collapses the state according to the measurement result. + * @param rootEdge the root edge of the state vector decision diagram + * @param index the index of the qubit to be measured + * @param probability the probability of the measurement result (required for + * normalization) + * @param measureZero whether or not to measure '0' (otherwise '1' is + * measured) */ - mEdge conjugateTranspose(const mEdge& a); + void performCollapsingMeasurement(vEdge& rootEdge, Qubit index, + fp probability, bool measureZero); + /** - * @brief Recursively computes the conjugate transpose of a given matrix edge. - * - * @param a The matrix edge to conjugate transpose. - * @return The conjugated transposed matrix edge. + * @brief Measures the qubit with the given index in the given state vector + * decision diagram. Collapses the state according to the measurement result. + * @param rootEdge the root edge of the state vector decision diagram + * @param index the index of the qubit to be measured + * @param mt the random number generator + * @param epsilon the numerical precision used for checking the normalization + * of the state vector decision diagram + * @return the measurement result ('0' or '1') + * @throws std::runtime_error if a numerical instability is detected during + * the measurement. */ - mCachedEdge conjugateTransposeRec(const mEdge& a); + char measureOneCollapsing(vEdge& rootEdge, Qubit index, std::mt19937_64& mt, + fp epsilon = 0.001); + + char measureOneCollapsing(dEdge& e, Qubit index, std::mt19937_64& mt); +public: /// /// Multiplication /// @@ -1230,7 +508,7 @@ class Package { } else if (!m.w.exactlyZero()) { dNode::applyDmChangesToNode(edge[idx].p); dNode::applyDmChangesToNode(m.p); - edge[idx] = add2(edge[idx], m, v); + edge[idx] = add(edge[idx], m, v); dNode::revertDmChangesToNode(m.p); dNode::revertDmChangesToNode(edge[idx].p); } @@ -1242,7 +520,7 @@ class Package { if (k == 0 || edge[idx].w.exactlyZero()) { edge[idx] = m; } else if (!m.w.exactlyZero()) { - edge[idx] = add2(edge[idx], m, v); + edge[idx] = add(edge[idx], m, v); } } } @@ -1256,83 +534,6 @@ class Package { return e; } - /// - /// Inner product, fidelity, expectation value - /// -public: - ComputeTable vectorInnerProduct{ - config_.ctVecInnerProdNumBucket}; - - /** - * @brief Calculates the inner product of two vector decision diagrams. - * - * @param x A vector DD representing a quantum state. - * @param y A vector DD representing a quantum state. - * @return A complex number representing the scalar product of the DDs. - */ - ComplexValue innerProduct(const vEdge& x, const vEdge& y); - - /** - * @brief Calculates the fidelity between two vector decision diagrams. - * - * @param x A vector DD representing a quantum state. - * @param y A vector DD representing a quantum state. - * @return The fidelity between the two quantum states. - */ - fp fidelity(const vEdge& x, const vEdge& y); - - /** - * @brief Calculates the fidelity between a vector decision diagram and a - * sparse probability vector. - * - * @details This function computes the fidelity between a quantum state - * represented by a vector decision diagram and a sparse probability vector. - * The optional permutation of qubits can be provided to match the qubit - * ordering. - * - * @param e The root edge of the decision diagram. - * @param probs A map of probabilities for each measurement outcome. - * @param permutation An optional permutation of qubits. - * @return The fidelity of the measurement outcomes. - */ - static fp - fidelityOfMeasurementOutcomes(const vEdge& e, const SparsePVec& probs, - const qc::Permutation& permutation = {}); - -private: - /** - * @brief Recursively calculates the inner product of two vector decision - * diagrams. - * - * @param x A vector DD representing a quantum state. - * @param y A vector DD representing a quantum state. - * @param var The number of levels contained in each vector DD. - * @return A complex number representing the scalar product of the DDs. - * @note This function is called recursively such that the number of levels - * decreases each time to traverse the DDs. - */ - ComplexValue innerProduct(const vEdge& x, const vEdge& y, Qubit var); - - /** - * @brief Recursively calculates the fidelity of measurement outcomes. - * - * @details This function computes the fidelity between a quantum state - * represented by a vector decision diagram and a sparse probability vector. - * It traverses the decision diagram recursively, calculating the contribution - * of each path to the overall fidelity. An optional permutation of qubits can - * be provided to match the qubit ordering. - * - * @param e The root edge of the decision diagram. - * @param probs A map of probabilities for each measurement outcome. - * @param i The current index in the decision diagram traversal. - * @param permutation An optional permutation of qubits. - * @param nQubits The number of qubits in the decision diagram. - * @return The fidelity of the measurement outcomes. - */ - static fp fidelityOfMeasurementOutcomesRecursive( - const vEdge& e, const SparsePVec& probs, std::size_t i, - const qc::Permutation& permutation, std::size_t nQubits); - public: /** * @brief Calculates the expectation value of an operator with respect to a @@ -1353,141 +554,6 @@ class Package { */ fp expectationValue(const mEdge& x, const vEdge& y); - /// - /// Kronecker/tensor product - /// - - ComputeTable vectorKronecker{ - config_.ctVecKronNumBucket}; - ComputeTable matrixKronecker{ - config_.ctMatKronNumBucket}; - - /** - * @brief Get the compute table for Kronecker product operations. - * - * @tparam Node The type of the node. - * @return A reference to the appropriate compute table for the given node - * type. - */ - template [[nodiscard]] auto& getKroneckerComputeTable() { - if constexpr (std::is_same_v) { - return vectorKronecker; - } else { - return matrixKronecker; - } - } - - /** - * @brief Computes the Kronecker product of two decision diagrams. - * - * @tparam Node The type of the node. - * @param x The first decision diagram. - * @param y The second decision diagram. - * @param yNumQubits The number of qubits in the second decision diagram. - * @param incIdx Whether to increment the index of the nodes in the second - * 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); - } - -private: - /** - * @brief Internal function to compute the Kronecker product of two decision - * diagrams. - * - * This function is used internally to compute the Kronecker product of two - * decision diagrams (DDs) of type Node. It is not intended to be called - * directly. - * - * @tparam Node The type of the node. - * @param x The first decision diagram. - * @param y The second decision diagram. - * @param yNumQubits The number of qubits in the second decision diagram. - * @param incIdx Whether to increment the qubit index. - * @return The resulting decision diagram after the Kronecker product. - */ - template - CachedEdge kronecker2(const Edge& x, const Edge& y, - const std::size_t yNumQubits, - const bool incIdx = true) { - if (x.w.exactlyZero() || y.w.exactlyZero()) { - return CachedEdge::zero(); - } - const auto xWeight = static_cast(x.w); - if (xWeight.approximatelyZero()) { - return CachedEdge::zero(); - } - const auto yWeight = static_cast(y.w); - if (yWeight.approximatelyZero()) { - return CachedEdge::zero(); - } - const auto rWeight = xWeight * yWeight; - if (rWeight.approximatelyZero()) { - return CachedEdge::zero(); - } - - if (x.isTerminal() && y.isTerminal()) { - return {x.p, rWeight}; - } - - if constexpr (std::is_same_v || std::is_same_v) { - if (x.isIdentity()) { - return {y.p, rWeight}; - } - } else { - if (x.isTerminal()) { - return {y.p, rWeight}; - } - if (y.isTerminal()) { - return {x.p, rWeight}; - } - } - - // check if we already computed the product before and return the result - auto& computeTable = getKroneckerComputeTable(); - if (const auto* r = computeTable.lookup(x.p, y.p); r != nullptr) { - return {r->p, rWeight}; - } - - constexpr std::size_t n = std::tuple_size_ve)>; - std::array, n> edge{}; - for (auto i = 0U; i < n; ++i) { - edge[i] = kronecker2(x.p->e[i], y, yNumQubits, incIdx); - } - - // Increase the qubit index - 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 (y.isIdentity()) { - idx += static_cast(yNumQubits); - } else { - idx += static_cast(y.p->v + 1U); - } - } else { - idx += static_cast(y.p->v + 1U); - } - } - auto e = makeDDNode(idx, edge, true); - computeTable.insert(x.p, y.p, {e.p, e.w}); - return {e.p, rWeight}; - } - /// /// (Partial) trace /// @@ -1538,23 +604,6 @@ class Package { return trace(a, eliminate, numQubits).w; } - /** - * @brief Checks if a given matrix is close to the identity matrix. - * @details This function checks if a given matrix is close to the identity - * matrix, while ignoring any potential garbage qubits and ignoring the - * diagonal weights if `checkCloseToOne` is set to false. - * @param m An mEdge that represents the DD of the matrix. - * @param tol The accepted tolerance for the edge weights of the DD. - * @param garbage A vector of boolean values that defines which qubits are - * considered garbage qubits. If it's empty, then no qubit is considered to be - * a garbage qubit. - * @param checkCloseToOne If false, the function only checks if the matrix is - * close to a diagonal matrix. - */ - [[nodiscard]] bool isCloseToIdentity(const mEdge& m, fp tol = 1e-10, - const std::vector& garbage = {}, - bool checkCloseToOne = true) const; - private: /** * @brief Computes the normalized (partial) trace using a compute table to @@ -1611,8 +660,8 @@ class Package { } 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); + auto r = add(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 @@ -1645,373 +694,13 @@ class Package { return r; } - /** - * @brief Recursively checks if a given matrix is close to the identity - * matrix. - * - * @param m The matrix edge to check. - * @param visited A set of visited nodes to avoid redundant checks. - * @param tol The tolerance for comparing edge weights. - * @param garbage A vector of boolean values indicating which qubits are - * considered garbage. - * @param checkCloseToOne A flag to indicate whether to check if diagonal - * elements are close to one. - * @return True if the matrix is close to the identity matrix, false - * otherwise. - */ - static bool isCloseToIdentityRecursive( - const mEdge& m, std::unordered_set& visited, fp tol, - const std::vector& garbage, bool checkCloseToOne); - public: - /// - /// Identity matrices - /// - - /// Create identity DD represented by the one-terminal. - static mEdge makeIdent(); - - mEdge createInitialMatrix(const std::vector& ancillary); - /// /// Noise Operations /// StochasticNoiseOperationTable stochasticNoiseOperationCache{ nqubits, config_.stochasticCacheOps}; DensityNoiseTable densityNoise{config_.ctDmNoiseNumBucket}; - - /// - /// Ancillary and garbage reduction - /// - - /** - * @brief Reduces the decision diagram by handling ancillary qubits. - * - * @param e The matrix decision diagram edge to be reduced. - * @param ancillary A boolean vector indicating which qubits are ancillary - * (true) or not (false). - * @param regular Flag indicating whether to perform regular (true) or inverse - * (false) reduction. - * @return The reduced matrix decision diagram edge. - * - * @details This function modifies the decision diagram to account for - * ancillary qubits by: - * 1. Early returning if there are no ancillary qubits or if the edge is zero - * 2. Special handling for identity matrices by creating appropriate zero - * nodes - * 3. Finding the lowest ancillary qubit as a starting point - * 4. Recursively reducing nodes starting from the lowest ancillary qubit - * 5. Adding zero nodes for any remaining higher ancillary qubits - * - * The function maintains proper reference counting by incrementing the - * reference count of the result and decrementing the reference count of the - * input edge. - */ - mEdge reduceAncillae(mEdge e, const std::vector& ancillary, - bool regular = true); - - /** - * @brief Reduces the given decision diagram by summing entries for garbage - * qubits. - * - * For each garbage qubit q, this function sums all the entries for q = 0 and - * q = 1, setting the entry for q = 0 to the sum and the entry for q = 1 to - * zero. To ensure that the probabilities of the resulting state are the sum - * of the probabilities of the initial state, the function computes - * `sqrt(|a|^2 + |b|^2)` for two entries `a` and `b`. - * - * @param e DD representation of the matrix/vector. - * @param garbage Vector that describes which qubits are garbage and which - * ones are not. If garbage[i] = true, then qubit q_i is considered garbage. - * @param normalizeWeights By default set to `false`. If set to `true`, the - * function changes all weights in the DD to their magnitude, also for - * non-garbage qubits. This is used for checking - * partial equivalence of circuits. For partial equivalence, only the - * measurement probabilities are considered, so we - * need to consider only the magnitudes of each entry. - * @return DD representing the reduced matrix/vector. - */ - vEdge reduceGarbage(vEdge& e, const std::vector& garbage, - bool normalizeWeights = false); - - /** - * @brief Reduces garbage qubits in a matrix decision diagram. - * - * @param e The matrix decision diagram edge to be reduced. - * @param garbage A boolean vector indicating which qubits are garbage (true) - * or not (false). - * @param regular Flag indicating whether to apply regular (true) or inverse - * (false) reduction. In regular mode, garbage entries are summed in the first - * two components, in inverse mode, they are summed in the first and third - * components. - * @param normalizeWeights Flag indicating whether to normalize weights to - * their magnitudes. When true, all weights in the DD are changed to their - * magnitude, also for non-garbage qubits. This is used for checking partial - * equivalence where only measurement probabilities matter. - * @return The reduced matrix decision diagram edge. - * - * @details For each garbage qubit q, this function sums all the entries for - * q=0 and q=1, setting the entry for q=0 to the sum and the entry for q=1 to - * zero. To maintain proper probabilities, the function computes sqrt(|a|^2 + - * |b|^2) for two entries a and b. The function handles special cases like - * zero terminals and identity matrices separately and maintains proper - * reference counting throughout the reduction process. - */ - mEdge reduceGarbage(const mEdge& e, const std::vector& garbage, - bool regular = true, bool normalizeWeights = false); - -private: - mCachedEdge reduceAncillaeRecursion(mNode* p, - const std::vector& ancillary, - Qubit lowerbound, bool regular = true); - - vCachedEdge reduceGarbageRecursion(vNode* p, const std::vector& garbage, - Qubit lowerbound, - bool normalizeWeights = false); - mCachedEdge reduceGarbageRecursion(mNode* p, const std::vector& garbage, - Qubit lowerbound, bool regular = true, - bool normalizeWeights = false); - - /// - /// Vector and matrix extraction from DDs - /// -public: - /// transfers a decision diagram from another package to this package - template Edge transfer(Edge& original) { - if (original.isTerminal()) { - return {original.p, cn.lookup(original.w)}; - } - - // POST ORDER TRAVERSAL USING ONE STACK - // https://www.geeksforgeeks.org/iterative-postorder-traversal-using-stack/ - Edge root{}; - std::stack*> stack; - - std::unordered_map mappedNode{}; - - Edge* currentEdge = &original; - constexpr std::size_t n = std::tuple_size_ve)>; - // NOLINTNEXTLINE(cppcoreguidelines-avoid-do-while) - do { - while (currentEdge != nullptr && !currentEdge->isTerminal()) { - for (std::size_t i = n - 1; i > 0; --i) { - auto& edge = currentEdge->p->e[i]; - if (edge.isTerminal()) { - continue; - } - if (edge.w.approximatelyZero()) { - continue; - } - if (mappedNode.find(edge.p) != mappedNode.end()) { - continue; - } - - // non-zero edge to be included - stack.push(&edge); - } - stack.push(currentEdge); - currentEdge = ¤tEdge->p->e[0]; - } - currentEdge = stack.top(); - stack.pop(); - - bool hasChild = false; - for (std::size_t i = 1; i < n && !hasChild; ++i) { - auto& edge = currentEdge->p->e[i]; - if (edge.w.approximatelyZero()) { - continue; - } - if (mappedNode.find(edge.p) != mappedNode.end()) { - continue; - } - hasChild = edge.p == stack.top()->p; - } - - if (hasChild) { - Edge* temp = stack.top(); - stack.pop(); - stack.push(currentEdge); - currentEdge = temp; - } else { - if (mappedNode.find(currentEdge->p) != mappedNode.end()) { - currentEdge = nullptr; - continue; - } - std::array, n> edges{}; - for (std::size_t i = 0; i < n; i++) { - if (currentEdge->p->e[i].isTerminal()) { - edges[i].p = currentEdge->p->e[i].p; - } else { - edges[i].p = mappedNode[currentEdge->p->e[i].p]; - } - edges[i].w = cn.lookup(currentEdge->p->e[i].w); - } - root = makeDDNode(currentEdge->p->v, edges); - mappedNode[currentEdge->p] = root.p; - currentEdge = nullptr; - } - } while (!stack.empty()); - root.w = cn.lookup(original.w * root.w); - return root; - } - - /// - /// Deserialization - /// Note: do not rely on the binary format being portable across different - /// architectures/platforms - /// - - template , - std::size_t N = std::tuple_size_v> - Edge deserialize(std::istream& is, const bool readBinary = false) { - auto result = CachedEdge{}; - ComplexValue rootweight{}; - - std::unordered_map nodes{}; - std::int64_t nodeIndex{}; - Qubit v{}; - std::array edgeWeights{}; - std::array edgeIndices{}; - edgeIndices.fill(-2); - - if (readBinary) { - std::remove_const_t version{}; - is.read(reinterpret_cast(&version), - sizeof(decltype(SERIALIZATION_VERSION))); - if (version != SERIALIZATION_VERSION) { - throw std::runtime_error( - "Wrong Version of serialization file version. version of file: " + - std::to_string(version) + - "; current version: " + std::to_string(SERIALIZATION_VERSION)); - } - - if (!is.eof()) { - rootweight.readBinary(is); - } - - while (is.read(reinterpret_cast(&nodeIndex), - sizeof(decltype(nodeIndex)))) { - is.read(reinterpret_cast(&v), sizeof(decltype(v))); - for (std::size_t i = 0U; i < N; i++) { - is.read(reinterpret_cast(&edgeIndices[i]), - sizeof(decltype(edgeIndices[i]))); - edgeWeights[i].readBinary(is); - } - result = deserializeNode(nodeIndex, v, edgeIndices, edgeWeights, nodes); - } - } else { - std::string version; - std::getline(is, version); - if (std::stoi(version) != SERIALIZATION_VERSION) { - throw std::runtime_error( - "Wrong Version of serialization file version. version of file: " + - version + - "; current version: " + std::to_string(SERIALIZATION_VERSION)); - } - - const std::string complexRealRegex = - R"(([+-]?(?:\d+(?:\.\d*)?|\.\d+)(?:[eE][+-]?\d+)?(?![ \d\.]*(?:[eE][+-])?\d*[iI]))?)"; - const std::string complexImagRegex = - R"(( ?[+-]? ?(?:(?:\d+(?:\.\d*)?|\.\d+)(?:[eE][+-]?\d+)?)?[iI])?)"; - const std::string edgeRegex = - " \\(((-?\\d+) (" + complexRealRegex + complexImagRegex + "))?\\)"; - const std::regex complexWeightRegex(complexRealRegex + complexImagRegex); - - std::string lineConstruct = "(\\d+) (\\d+)"; - for (std::size_t i = 0U; i < N; ++i) { - lineConstruct += "(?:" + edgeRegex + ")"; - } - lineConstruct += " *(?:#.*)?"; - const std::regex lineRegex(lineConstruct); - std::smatch m; - - std::string line; - if (std::getline(is, line)) { - if (!std::regex_match(line, m, complexWeightRegex)) { - throw std::runtime_error("Regex did not match second line: " + line); - } - rootweight.fromString(m.str(1), m.str(2)); - } - - while (std::getline(is, line)) { - if (line.empty() || line.size() == 1) { - continue; - } - - if (!std::regex_match(line, m, lineRegex)) { - throw std::runtime_error("Regex did not match line: " + line); - } - - // match 1: node_idx - // match 2: qubit_idx - - // repeats for every edge - // match 3: edge content - // match 4: edge_target_idx - // match 5: real + imag (without i) - // match 6: real - // match 7: imag (without i) - nodeIndex = std::stoi(m.str(1)); - v = static_cast(std::stoi(m.str(2))); - - for (auto edgeIdx = 3U, i = 0U; i < N; i++, edgeIdx += 5) { - if (m.str(edgeIdx).empty()) { - continue; - } - - edgeIndices[i] = std::stoi(m.str(edgeIdx + 1)); - edgeWeights[i].fromString(m.str(edgeIdx + 3), m.str(edgeIdx + 4)); - } - - result = deserializeNode(nodeIndex, v, edgeIndices, edgeWeights, nodes); - } - } - return {result.p, cn.lookup(result.w * rootweight)}; - } - - template > - Edge deserialize(const std::string& inputFilename, const bool readBinary) { - auto ifs = std::ifstream(inputFilename, std::ios::binary); - - if (!ifs.good()) { - throw std::invalid_argument("Cannot open serialized file: " + - inputFilename); - } - - return deserialize(ifs, readBinary); - } - -private: - template > - CachedEdge - deserializeNode(const std::int64_t index, const Qubit v, - std::array& edgeIdx, - const std::array& edgeWeight, - std::unordered_map& nodes) { - if (index == -1) { - return CachedEdge::zero(); - } - - std::array, N> edges{}; - for (auto i = 0U; i < N; ++i) { - if (edgeIdx[i] == -2) { - edges[i] = CachedEdge::zero(); - } else { - if (edgeIdx[i] == -1) { - edges[i] = CachedEdge::one(); - } else { - edges[i].p = nodes[edgeIdx[i]]; - } - edges[i].w = edgeWeight[i]; - } - } - // reset - edgeIdx.fill(-2); - - auto r = makeDDNode(v, edges); - nodes[index] = r.p; - return r; - } }; } // namespace dd diff --git a/include/mqt-core/dd/SmallSpan.hpp b/include/mqt-core/dd/SmallSpan.hpp new file mode 100644 index 000000000..f5facbf35 --- /dev/null +++ b/include/mqt-core/dd/SmallSpan.hpp @@ -0,0 +1,89 @@ +/* + * Copyright (c) 2025 Chair for Design Automation, TUM + * All rights reserved. + * + * SPDX-License-Identifier: MIT + * + * Licensed under the MIT License + */ + +#pragma once + +#include +#include +#include + +namespace dd { + +template class SmallSpan { +public: + SmallSpan(T* begin, std::uint8_t size) : begin_(begin), size_(size) {} + + template + explicit SmallSpan(std::array& arr) + : begin_(arr.data()), size_(static_cast(N)) { + assert(N <= std::numeric_limits::max()); + } + + explicit SmallSpan(std::vector& vec) + : begin_(vec.data()), size_(static_cast(vec.size())) { + assert(vec.size() <= std::numeric_limits::max()); + } + + bool operator==(const SmallSpan& other) const { + if (size() != other.size()) { + return false; + } + return std::equal(begin(), end(), other.begin()); + } + + [[nodiscard]] inline constexpr const T* begin() const noexcept { + return begin_; + } + [[nodiscard]] inline constexpr const T* end() const noexcept { + return begin_ + size_; + } + + [[nodiscard]] inline constexpr const T* cbegin() const noexcept { + return begin(); + } + [[nodiscard]] inline constexpr const T* cend() const noexcept { + return end(); + } + + [[nodiscard]] inline constexpr T* begin() noexcept { return begin_; } + [[nodiscard]] inline constexpr T* end() noexcept { return begin_ + size_; } + + [[nodiscard]] inline constexpr std::uint8_t size() const noexcept { + return size_; + } + + [[nodiscard]] inline constexpr const T& + operator[](std::size_t i) const noexcept { + return begin_[i]; + } + + [[nodiscard]] inline constexpr T& operator[](std::size_t i) noexcept { + return begin_[i]; + } + + [[nodiscard]] inline constexpr const T& at(std::size_t i) const { + if (i >= size()) { + throw std::out_of_range("Index out of range"); + } + return begin_[i]; + } + + [[nodiscard]] inline constexpr T& at(std::size_t i) { + if (i >= size()) { + throw std::out_of_range("Index out of range"); + } + return begin_[i]; + } + +private: + T* begin_; + std::uint8_t size_; +}; + +} // namespace dd diff --git a/include/mqt-core/dd/UnaryComputeTable.hpp b/include/mqt-core/dd/UnaryComputeTable.hpp index 7e62b3628..1536c5d6f 100644 --- a/include/mqt-core/dd/UnaryComputeTable.hpp +++ b/include/mqt-core/dd/UnaryComputeTable.hpp @@ -31,6 +31,8 @@ namespace dd { template class UnaryComputeTable { public: + static constexpr std::size_t DEFAULT_NUM_BUCKETS = NBUCKET; + /// Default constructor explicit UnaryComputeTable(const size_t numBuckets = 32768U) { // numBuckets must be a power of two diff --git a/include/mqt-core/dd/UniqueTable.hpp b/include/mqt-core/dd/UniqueTable.hpp index 4681a3ddb..ffd4a07a4 100644 --- a/include/mqt-core/dd/UniqueTable.hpp +++ b/include/mqt-core/dd/UniqueTable.hpp @@ -178,7 +178,7 @@ class UniqueTable { void clear(); - template void print() { + template void print() const { static_assert(std::is_base_of_v, "Node must be derived from NodeBase"); auto q = cfg.nVars - 1U; diff --git a/include/mqt-core/dd/VectorDDContainer.hpp b/include/mqt-core/dd/VectorDDContainer.hpp new file mode 100644 index 000000000..4b5194f2e --- /dev/null +++ b/include/mqt-core/dd/VectorDDContainer.hpp @@ -0,0 +1,289 @@ +/* + * Copyright (c) 2025 Chair for Design Automation, TUM + * All rights reserved. + * + * SPDX-License-Identifier: MIT + * + * Licensed under the MIT License + */ + +#pragma once + +#include "dd/DesicionDiagramContainer.hpp" +#include "dd/UnaryComputeTable.hpp" +#include "ir/Permutation.hpp" + +namespace dd { + +class VectorDDContainer : public DDContainer { +public: + struct Config : public DDContainer::Config { + std::size_t ctVecConjNumBucket = + UnaryComputeTable::DEFAULT_NUM_BUCKETS; + std::size_t vectorInnerProduct = + ComputeTable::DEFAULT_NUM_BUCKETS; + }; + + VectorDDContainer(std::size_t nqubits, RealNumberUniqueTable& cUt, + ComplexNumbers& cn, const Config& config) + : DDContainer(nqubits, cUt, cn, config), + conjugateVector(config.ctVecConjNumBucket), + vectorInnerProduct(config.vectorInnerProduct) {} + + void reset(); + + bool garbageCollect(bool force); + + /** + * @brief Construct the all-zero state \f$|0...0\rangle\f$ + * @param n The number of qubits + * @param start The starting qubit index. Default is 0. + * @return A decision diagram for the all-zero state + */ + vEdge makeZeroState(std::size_t n, std::size_t start = 0); + + /** + * @brief Construct a computational basis state \f$|b_{n-1}...b_0\rangle\f$ + * @param n The number of qubits + * @param state The state to construct + * @param start The starting qubit index. Default is 0. + * @return A decision diagram for the computational basis state + */ + vEdge makeBasisState(std::size_t n, const std::vector& state, + std::size_t start = 0); + + /** + * @brief Construct a product state out of + * \f$\{0, 1, +, -, R, L\}^{\otimes n}\f$. + * @param n The number of qubits + * @param state The state to construct + * @param start The starting qubit index. Default is 0. + * @return A decision diagram for the product state + */ + vEdge makeBasisState(std::size_t n, const std::vector& state, + std::size_t start = 0); + + /** + * @brief Construct a GHZ state \f$|0...0\rangle + |1...1\rangle\f$ + * @param n The number of qubits + * @return A decision diagram for the GHZ state + */ + vEdge makeGHZState(std::size_t n); + + /** + * @brief Construct a W state + * @details The W state is defined as + * \f[ + * |0...01\rangle + |0...10\rangle + |10...0\rangle + * \f] + * @param n The number of qubits + * @return A decision diagram for the W state + */ + vEdge makeWState(std::size_t n); + + /** + * @brief Construct a decision diagram from an arbitrary state vector + * @param stateVector The state vector to convert to a DD + * @return A decision diagram for the state + */ + vEdge makeStateFromVector(const CVec& stateVector); + + /** + * @brief Measure all qubits in the given decision diagram. + * + * @details This function measures all qubits in the decision diagram + * represented by `rootEdge`. It checks for numerical instabilities and + * collapses the state if requested. + * + * @param rootEdge The decision diagram to measure. + * @param collapse If true, the state is collapsed after measurement. + * @param mt A random number generator. + * @param epsilon The tolerance for numerical instabilities. + * @return A string representing the measurement result. + * @throws std::runtime_error If numerical instabilities are detected or if + * probabilities do not sum to 1. + */ + std::string measureAll(vEdge& rootEdge, bool collapse, std::mt19937_64& mt, + fp epsilon = 0.001); + + /** + * @brief Determine the measurement probabilities for a given qubit index. + * + * @param rootEdge The root edge of the decision diagram. + * @param index The qubit index to determine the measurement probabilities + * for. + * @return A pair of floating-point values representing the probabilities of + * measuring 0 and 1, respectively. + * + * @details This function calculates the probabilities of measuring 0 and 1 + * for a given qubit index in the decision diagram. It uses a breadth-first + * search to traverse the decision diagram and accumulate the measurement + * probabilities. The function maintains a map of measurement probabilities + * for each node and a set of visited nodes to avoid redundant calculations. + * It also uses a queue to process nodes level by level. + */ + static std::pair + determineMeasurementProbabilities(const vEdge& rootEdge, Qubit index); + +private: + /** + * @brief Assigns probabilities to nodes in a decision diagram. + * + * @details This function recursively assigns probabilities to nodes in a + * decision diagram. It calculates the probability of reaching each node and + * stores the result in a map. + * + * @param edge The edge to start the probability assignment from. + * @param probs A map to store the probabilities of each node. + * @return The probability of the given edge. + */ + static fp assignProbabilities(const vEdge& edge, + std::unordered_map& probs); + + /** + * @brief Constructs a decision diagram (DD) from a state vector using a + * recursive algorithm. + * + * @param begin Iterator pointing to the beginning of the state vector. + * @param end Iterator pointing to the end of the state vector. + * @param level The current level of recursion. Starts at the highest level of + * the state vector (log base 2 of the vector size - 1). + * @return A vCachedEdge representing the root node of the created DD. + * + * @details This function recursively breaks down the state vector into halves + * until each half has only one element. At each level of recursion, two new + * edges are created, one for each half of the state vector. The two resulting + * decision diagram edges are used to create a new decision diagram node at + * the current level, and this node is returned as the result of the current + * recursive call. At the base case of recursion, the state vector has only + * two elements, which are converted into terminal nodes of the decision + * diagram. + * + * @note This function assumes that the state vector size is a power of two. + */ + vCachedEdge makeStateFromVector(const CVec::const_iterator& begin, + const CVec::const_iterator& end, Qubit level); + +public: + /** + * @brief Conjugates a given decision diagram edge. + * + * @param a The decision diagram edge to conjugate. + * @return The conjugated decision diagram edge. + */ + vEdge conjugate(const vEdge& a); + + /** + * @brief Recursively conjugates a given decision diagram edge. + * + * @param a The decision diagram edge to conjugate. + * @return The conjugated decision diagram edge. + */ + vCachedEdge conjugateRec(const vEdge& a); + + /** + * @brief Calculates the inner product of two vector decision diagrams. + * + * @param x A vector DD representing a quantum state. + * @param y A vector DD representing a quantum state. + * @return A complex number representing the scalar product of the DDs. + */ + ComplexValue innerProduct(const vEdge& x, const vEdge& y); + + /** + * @brief Calculates the fidelity between two vector decision diagrams. + * + * @param x A vector DD representing a quantum state. + * @param y A vector DD representing a quantum state. + * @return The fidelity between the two quantum states. + */ + fp fidelity(const vEdge& x, const vEdge& y); + + /** + * @brief Calculates the fidelity between a vector decision diagram and a + * sparse probability vector. + * + * @details This function computes the fidelity between a quantum state + * represented by a vector decision diagram and a sparse probability vector. + * The optional permutation of qubits can be provided to match the qubit + * ordering. + * + * @param e The root edge of the decision diagram. + * @param probs A map of probabilities for each measurement outcome. + * @param permutation An optional permutation of qubits. + * @return The fidelity of the measurement outcomes. + */ + static fp + fidelityOfMeasurementOutcomes(const vEdge& e, const SparsePVec& probs, + const qc::Permutation& permutation = {}); + +private: + /** + * @brief Recursively calculates the inner product of two vector decision + * diagrams. + * + * @param x A vector DD representing a quantum state. + * @param y A vector DD representing a quantum state. + * @param var The number of levels contained in each vector DD. + * @return A complex number representing the scalar product of the DDs. + * @note This function is called recursively such that the number of levels + * decreases each time to traverse the DDs. + */ + ComplexValue innerProduct(const vEdge& x, const vEdge& y, Qubit var); + + /** + * @brief Recursively calculates the fidelity of measurement outcomes. + * + * @details This function computes the fidelity between a quantum state + * represented by a vector decision diagram and a sparse probability vector. + * It traverses the decision diagram recursively, calculating the contribution + * of each path to the overall fidelity. An optional permutation of qubits can + * be provided to match the qubit ordering. + * + * @param e The root edge of the decision diagram. + * @param probs A map of probabilities for each measurement outcome. + * @param i The current index in the decision diagram traversal. + * @param permutation An optional permutation of qubits. + * @param nQubits The number of qubits in the decision diagram. + * @return The fidelity of the measurement outcomes. + */ + static fp fidelityOfMeasurementOutcomesRecursive( + const vEdge& e, const SparsePVec& probs, std::size_t i, + const qc::Permutation& permutation, std::size_t nQubits); + +public: + /** + * @brief Reduces the given decision diagram by summing entries for garbage + * qubits. + * + * For each garbage qubit q, this function sums all the entries for q = 0 and + * q = 1, setting the entry for q = 0 to the sum and the entry for q = 1 to + * zero. To ensure that the probabilities of the resulting state are the sum + * of the probabilities of the initial state, the function computes + * `sqrt(|a|^2 + |b|^2)` for two entries `a` and `b`. + * + * @param e DD representation of the matrix/vector. + * @param garbage Vector that describes which qubits are garbage and which + * ones are not. If garbage[i] = true, then qubit q_i is considered garbage. + * @param normalizeWeights By default set to `false`. If set to `true`, the + * function changes all weights in the DD to their magnitude, also for + * non-garbage qubits. This is used for checking + * partial equivalence of circuits. For partial equivalence, only the + * measurement probabilities are considered, so we + * need to consider only the magnitudes of each entry. + * @return DD representing the reduced matrix/vector. + */ + [[nodiscard]] vEdge reduceGarbage(vEdge& e, const std::vector& garbage, + bool normalizeWeights = false); + +private: + vCachedEdge reduceGarbageRecursion(vNode* p, const std::vector& garbage, + Qubit lowerbound, + bool normalizeWeights = false); + + UnaryComputeTable conjugateVector; + + ComputeTable vectorInnerProduct; +}; + +} // namespace dd diff --git a/src/dd/DensityDDContainer.cpp b/src/dd/DensityDDContainer.cpp new file mode 100644 index 000000000..86ec5d79b --- /dev/null +++ b/src/dd/DensityDDContainer.cpp @@ -0,0 +1,26 @@ +/* + * Copyright (c) 2025 Chair for Design Automation, TUM + * All rights reserved. + * + * SPDX-License-Identifier: MIT + * + * Licensed under the MIT License + */ + +#include "dd/DensityDDContainer.hpp" + +#include "dd/GateMatrixDefinitions.hpp" + +namespace dd { + +dEdge DensityDDContainer::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; +} + +} // namespace dd diff --git a/src/dd/FunctionalityConstruction.cpp b/src/dd/FunctionalityConstruction.cpp index 6884b76d8..0bc25d36a 100644 --- a/src/dd/FunctionalityConstruction.cpp +++ b/src/dd/FunctionalityConstruction.cpp @@ -26,7 +26,7 @@ MatrixDD buildFunctionality(const qc::QuantumComputation& qc, Package& dd) { } auto permutation = qc.initialLayout; - auto e = dd.createInitialMatrix(qc.getAncillary()); + auto e = dd.matrices().createInitialMatrix(qc.getAncillary()); for (const auto& op : qc) { // SWAP gates can be executed virtually by changing the permutation @@ -40,7 +40,7 @@ MatrixDD buildFunctionality(const qc::QuantumComputation& qc, Package& dd) { } // correct permutation if necessary changePermutation(e, permutation, qc.outputPermutation, dd); - e = dd.reduceAncillae(e, qc.getAncillary()); + e = dd.matrices().reduceAncillae(e, qc.getAncillary()); e = dd.reduceGarbage(e, qc.getGarbage()); return e; @@ -53,7 +53,7 @@ bool buildFunctionalityRecursive(const qc::QuantumComputation& qc, qc::Permutation& permutation, Package& dd) { // base case if (depth == 1U) { - auto e = Package::makeIdent(); + auto e = MatrixDDContainer::makeIdent(); if (const auto& op = qc.at(opIdx); op->getType() == qc::OpType::SWAP && !op->isControlled()) { const auto& targets = op->getTargets(); @@ -68,7 +68,7 @@ bool buildFunctionalityRecursive(const qc::QuantumComputation& qc, dd.incRef(e); return false; } - auto f = Package::makeIdent(); + auto f = MatrixDDContainer::makeIdent(); if (const auto& op = qc.at(opIdx); op->getType() == qc::OpType::SWAP && !op->isControlled()) { const auto& targets = op->getTargets(); @@ -134,8 +134,8 @@ MatrixDD buildFunctionalityRecursive(const qc::QuantumComputation& qc, // correct permutation if necessary changePermutation(e, permutation, qc.outputPermutation, dd); - e = dd.reduceAncillae(e, qc.getAncillary()); - e = dd.reduceGarbage(e, qc.getGarbage()); + e = dd.matrices().reduceAncillae(e, qc.getAncillary()); + e = dd.matrices().reduceGarbage(e, qc.getGarbage()); return e; } diff --git a/src/dd/MatrixDDContainer.cpp b/src/dd/MatrixDDContainer.cpp new file mode 100644 index 000000000..a65884df5 --- /dev/null +++ b/src/dd/MatrixDDContainer.cpp @@ -0,0 +1,667 @@ +/* + * Copyright (c) 2025 Chair for Design Automation, TUM + * All rights reserved. + * + * SPDX-License-Identifier: MIT + * + * Licensed under the MIT License + */ + +#include "dd/MatrixDDContainer.hpp" + +namespace dd { + +mEdge MatrixDDContainer::makeGateDD(const GateMatrix& mat, + const qc::Qubit target) { + return makeGateDD(mat, qc::Controls{}, target); +} + +mEdge MatrixDDContainer::makeGateDD(const GateMatrix& mat, + const qc::Control& control, + const qc::Qubit target) { + return makeGateDD(mat, qc::Controls{control}, target); +} + +mEdge MatrixDDContainer::makeGateDD(const GateMatrix& mat, + const qc::Controls& controls, + const qc::Qubit target) { + if (std::any_of(controls.begin(), controls.end(), + [this](const auto& c) { + return c.qubit > static_cast(nqubits - 1U); + }) || + target > static_cast(nqubits - 1U)) { + throw std::runtime_error{ + "Requested gate acting on qubit(s) with index larger than " + + std::to_string(nqubits - 1U) + + " while the package configuration only supports up to " + + std::to_string(nqubits) + + " qubits. Please allocate a larger package instance."}; + } + std::array em{}; + for (auto i = 0U; i < NEDGE; ++i) { + em[i] = mCachedEdge::terminal(mat[i]); + } + + if (controls.empty()) { + // Single qubit operation + const auto e = makeDDNode(static_cast(target), em); + return {e.p, getCn().lookup(e.w)}; + } + + auto it = controls.begin(); + auto edges = std::array{mCachedEdge::zero(), mCachedEdge::zero(), + mCachedEdge::zero(), mCachedEdge::zero()}; + + // process lines below target + for (; it != controls.end() && it->qubit < target; ++it) { + for (auto i1 = 0U; i1 < RADIX; ++i1) { + for (auto i2 = 0U; i2 < RADIX; ++i2) { + const auto i = (i1 * RADIX) + i2; + if (it->type == qc::Control::Type::Neg) { // neg. control + edges[0] = em[i]; + edges[3] = (i1 == i2) ? mCachedEdge::one() : mCachedEdge::zero(); + } else { // pos. control + edges[0] = (i1 == i2) ? mCachedEdge::one() : mCachedEdge::zero(); + edges[3] = em[i]; + } + em[i] = makeDDNode(static_cast(it->qubit), edges); + } + } + } + + // target line + auto e = makeDDNode(static_cast(target), em); + + // process lines above target + for (; it != controls.end(); ++it) { + if (it->type == qc::Control::Type::Neg) { // neg. control + edges[0] = e; + edges[3] = mCachedEdge::one(); + e = makeDDNode(static_cast(it->qubit), edges); + } else { // pos. control + edges[0] = mCachedEdge::one(); + edges[3] = e; + e = makeDDNode(static_cast(it->qubit), edges); + } + } + return {e.p, getCn().lookup(e.w)}; +} + +mEdge MatrixDDContainer::makeTwoQubitGateDD(const TwoQubitGateMatrix& mat, + const qc::Qubit target0, + const qc::Qubit target1) { + return makeTwoQubitGateDD(mat, qc::Controls{}, target0, target1); +} + +mEdge MatrixDDContainer::makeTwoQubitGateDD(const TwoQubitGateMatrix& mat, + const qc::Control& control, + const qc::Qubit target0, + const qc::Qubit target1) { + return makeTwoQubitGateDD(mat, qc::Controls{control}, target0, target1); +} + +mEdge MatrixDDContainer::makeTwoQubitGateDD(const TwoQubitGateMatrix& mat, + const qc::Controls& controls, + const qc::Qubit target0, + const qc::Qubit target1) { + // sanity check + if (std::any_of(controls.begin(), controls.end(), + [this](const auto& c) { + return c.qubit > static_cast(nqubits - 1U); + }) || + target0 > static_cast(nqubits - 1U) || + target1 > static_cast(nqubits - 1U)) { + throw std::runtime_error{ + "Requested gate acting on qubit(s) with index larger than " + + std::to_string(nqubits - 1U) + + " while the package configuration only supports up to " + + std::to_string(nqubits) + + " qubits. Please allocate a larger package instance."}; + } + + // create terminal edge matrix + std::array, NEDGE> em{}; + for (auto i1 = 0U; i1 < NEDGE; i1++) { + const auto& matRow = mat.at(i1); + auto& emRow = em.at(i1); + for (auto i2 = 0U; i2 < NEDGE; i2++) { + emRow.at(i2) = mCachedEdge::terminal(matRow.at(i2)); + } + } + + // process lines below smaller target + auto it = controls.begin(); + const auto smallerTarget = std::min(target0, target1); + + auto edges = std::array{mCachedEdge::zero(), mCachedEdge::zero(), + mCachedEdge::zero(), mCachedEdge::zero()}; + + for (; it != controls.end() && it->qubit < smallerTarget; ++it) { + for (auto row = 0U; row < NEDGE; ++row) { + for (auto col = 0U; col < NEDGE; ++col) { + if (it->type == qc::Control::Type::Neg) { // negative control + edges[0] = em[row][col]; + edges[3] = (row == col) ? mCachedEdge::one() : mCachedEdge::zero(); + } else { // positive control + edges[0] = (row == col) ? mCachedEdge::one() : mCachedEdge::zero(); + edges[3] = em[row][col]; + } + em[row][col] = makeDDNode(static_cast(it->qubit), edges); + } + } + } + + // process the smaller target by taking the 16 submatrices and appropriately + // combining them into four DDs. + std::array em0{}; + for (std::size_t row = 0; row < RADIX; ++row) { + for (std::size_t col = 0; col < RADIX; ++col) { + std::array local{}; + if (target0 > target1) { + for (std::size_t i = 0; i < RADIX; ++i) { + for (std::size_t j = 0; j < RADIX; ++j) { + local.at((i * RADIX) + j) = + em.at((row * RADIX) + i).at((col * RADIX) + j); + } + } + } else { + for (std::size_t i = 0; i < RADIX; ++i) { + for (std::size_t j = 0; j < RADIX; ++j) { + local.at((i * RADIX) + j) = + em.at((i * RADIX) + row).at((j * RADIX) + col); + } + } + } + em0.at((row * RADIX) + col) = + makeDDNode(static_cast(smallerTarget), local); + } + } + + // process lines between the two targets + const auto largerTarget = std::max(target0, target1); + for (; it != controls.end() && it->qubit < largerTarget; ++it) { + for (auto i = 0U; i < NEDGE; ++i) { + if (it->type == qc::Control::Type::Neg) { // negative control + edges[0] = em0[i]; + edges[3] = + (i == 0 || i == 3) ? mCachedEdge::one() : mCachedEdge::zero(); + } else { // positive control + edges[0] = + (i == 0 || i == 3) ? mCachedEdge::one() : mCachedEdge::zero(); + edges[3] = em0[i]; + } + em0[i] = makeDDNode(static_cast(it->qubit), edges); + } + } + + // process the larger target by combining the four DDs from the smaller + // target + auto e = makeDDNode(static_cast(largerTarget), em0); + + // process lines above the larger target + for (; it != controls.end(); ++it) { + if (it->type == qc::Control::Type::Neg) { // negative control + edges[0] = e; + edges[3] = mCachedEdge::one(); + } else { // positive control + edges[0] = mCachedEdge::one(); + edges[3] = e; + } + e = makeDDNode(static_cast(it->qubit), edges); + } + + return {e.p, getCn().lookup(e.w)}; +} + +mEdge MatrixDDContainer::makeDDFromMatrix(const CMat& matrix) { + if (matrix.empty()) { + return mEdge::one(); + } + + const auto& length = matrix.size(); + if ((length & (length - 1)) != 0) { + throw std::invalid_argument("Matrix must have a length of a power of two."); + } + + const auto& width = matrix[0].size(); + if (length != width) { + throw std::invalid_argument("Matrix must be square."); + } + + if (length == 1) { + return mEdge::terminal(getCn().lookup(matrix[0][0])); + } + + const auto level = static_cast(std::log2(length) - 1); + const auto MatrixDDContainer = + makeDDFromMatrix(matrix, level, 0, length, 0, width); + return {MatrixDDContainer.p, getCn().lookup(MatrixDDContainer.w)}; +} + +mCachedEdge MatrixDDContainer::makeDDFromMatrix(const CMat& matrix, + const Qubit level, + const std::size_t rowStart, + const std::size_t rowEnd, + const std::size_t colStart, + const std::size_t colEnd) { + // base case + if (level == 0U) { + assert(rowEnd - rowStart == 2); + assert(colEnd - colStart == 2); + return makeDDNode( + 0U, {mCachedEdge::terminal(matrix[rowStart][colStart]), + mCachedEdge::terminal(matrix[rowStart][colStart + 1]), + mCachedEdge::terminal(matrix[rowStart + 1][colStart]), + mCachedEdge::terminal(matrix[rowStart + 1][colStart + 1])}); + } + + // recursively call the function on all quadrants + const auto rowMid = (rowStart + rowEnd) / 2; + const auto colMid = (colStart + colEnd) / 2; + const auto l = static_cast(level - 1U); + + return makeDDNode( + level, {makeDDFromMatrix(matrix, l, rowStart, rowMid, colStart, colMid), + makeDDFromMatrix(matrix, l, rowStart, rowMid, colMid, colEnd), + makeDDFromMatrix(matrix, l, rowMid, rowEnd, colStart, colMid), + makeDDFromMatrix(matrix, l, rowMid, rowEnd, colMid, colEnd)}); +} + +mEdge MatrixDDContainer::conjugateTranspose(const mEdge& a) { + const auto r = conjugateTransposeRec(a); + return {r.p, getCn().lookup(r.w)}; +} + +mCachedEdge MatrixDDContainer::conjugateTransposeRec(const mEdge& a) { + if (a.isTerminal()) { // terminal case + return {a.p, ComplexNumbers::conj(a.w)}; + } + + // check if in compute table + if (const auto* r = conjugateMatrixTranspose.lookup(a.p); r != nullptr) { + return {r->p, r->w * ComplexNumbers::conj(a.w)}; + } + + std::array e{}; + // conjugate transpose submatrices and rearrange as required + for (auto i = 0U; i < RADIX; ++i) { + for (auto j = 0U; j < RADIX; ++j) { + e[(RADIX * i) + j] = conjugateTransposeRec(a.p->e[(RADIX * j) + i]); + } + } + // create new top node + auto res = makeDDNode(a.p->v, e); + + // put it in the compute table + conjugateMatrixTranspose.insert(a.p, res); + + // adjust top weight including conjugate + return {res.p, res.w * ComplexNumbers::conj(a.w)}; +} + +bool MatrixDDContainer::isCloseToIdentity(const mEdge& m, const fp tol, + const std::vector& garbage, + const bool checkCloseToOne) const { + std::unordered_set visited{}; + visited.reserve(ut.getNumActiveEntries()); + return isCloseToIdentityRecursive(m, visited, tol, garbage, checkCloseToOne); +} +bool MatrixDDContainer::isCloseToIdentityRecursive( + const mEdge& m, std::unordered_set& visited, const fp tol, + const std::vector& garbage, const bool checkCloseToOne) { + // immediately return if this node is identical to the identity or zero + if (m.isTerminal()) { + return true; + } + + // immediately return if this node has already been visited + if (visited.find(m.p) != visited.end()) { + return true; + } + + const auto n = m.p->v; + + if (garbage.size() > n && garbage[n]) { + return isCloseToIdentityRecursive(m.p->e[0U], visited, tol, garbage, + checkCloseToOne) && + isCloseToIdentityRecursive(m.p->e[1U], visited, tol, garbage, + checkCloseToOne) && + isCloseToIdentityRecursive(m.p->e[2U], visited, tol, garbage, + checkCloseToOne) && + isCloseToIdentityRecursive(m.p->e[3U], visited, tol, garbage, + checkCloseToOne); + } + + // check whether any of the middle successors is non-zero, i.e., m = [ x 0 0 + // y ] + const auto mag1 = dd::ComplexNumbers::mag2(m.p->e[1U].w); + const auto mag2 = dd::ComplexNumbers::mag2(m.p->e[2U].w); + if (mag1 > tol || mag2 > tol) { + return false; + } + + if (checkCloseToOne) { + // check whether m = [ ~1 0 0 y ] + const auto mag0 = dd::ComplexNumbers::mag2(m.p->e[0U].w); + if (std::abs(mag0 - 1.0) > tol) { + return false; + } + const auto arg0 = dd::ComplexNumbers::arg(m.p->e[0U].w); + if (std::abs(arg0) > tol) { + return false; + } + + // check whether m = [ x 0 0 ~1 ] or m = [ x 0 0 ~0 ] (the last case is + // true for an ancillary qubit) + const auto mag3 = dd::ComplexNumbers::mag2(m.p->e[3U].w); + if (mag3 > tol) { + if (std::abs(mag3 - 1.0) > tol) { + return false; + } + const auto arg3 = dd::ComplexNumbers::arg(m.p->e[3U].w); + if (std::abs(arg3) > tol) { + return false; + } + } + } + // m either has the form [ ~1 0 0 ~1 ] or [ ~1 0 0 ~0 ] + const auto ident0 = isCloseToIdentityRecursive(m.p->e[0U], visited, tol, + garbage, checkCloseToOne); + + if (!ident0) { + return false; + } + // m either has the form [ I 0 0 ~1 ] or [ I 0 0 ~0 ] + const auto ident3 = isCloseToIdentityRecursive(m.p->e[3U], visited, tol, + garbage, checkCloseToOne); + + visited.insert(m.p); + return ident3; +} +mEdge MatrixDDContainer::makeIdent() { return mEdge::one(); } +mEdge MatrixDDContainer::createInitialMatrix( + const std::vector& ancillary) { + return reduceAncillae(makeIdent(), ancillary); +} +mEdge MatrixDDContainer::reduceAncillae(mEdge e, + const std::vector& ancillary, + const bool regular) { + // return if no more ancillaries left + if (std::none_of(ancillary.begin(), ancillary.end(), + [](const bool v) { return v; }) || + e.isZeroTerminal()) { + return e; + } + + // if we have only identities and no other nodes + if (e.isIdentity()) { + auto g = e; + for (auto i = 0U; i < ancillary.size(); ++i) { + if (ancillary[i]) { + g = makeDDNode( + static_cast(i), + std::array{g, mEdge::zero(), mEdge::zero(), mEdge::zero()}); + } + } + incRef(g); + return g; + } + + Qubit lowerbound = 0; + for (auto i = 0U; i < ancillary.size(); ++i) { + if (ancillary[i]) { + lowerbound = static_cast(i); + break; + } + } + + auto g = CachedEdge{e.p, 1.}; + if (e.p->v >= lowerbound) { + g = reduceAncillaeRecursion(e.p, ancillary, lowerbound, regular); + } + + for (std::size_t i = e.p->v + 1; i < ancillary.size(); ++i) { + if (ancillary[i]) { + g = makeDDNode(static_cast(i), + std::array{g, mCachedEdge::zero(), mCachedEdge::zero(), + mCachedEdge::zero()}); + } + } + const auto res = mEdge{g.p, getCn().lookup(g.w * e.w)}; + incRef(res); + decRef(e); + return res; +} + +mEdge MatrixDDContainer::reduceGarbage(const mEdge& e, + const std::vector& garbage, + const bool regular, + const bool normalizeWeights) { + // return if no more garbage left + if (!normalizeWeights && + (std::none_of(garbage.begin(), garbage.end(), [](bool v) { return v; }) || + e.isZeroTerminal())) { + return e; + } + + // if we have only identities and no other nodes + if (e.isIdentity()) { + auto g = e; + for (auto i = 0U; i < garbage.size(); ++i) { + if (garbage[i]) { + if (regular) { + g = makeDDNode(static_cast(i), + std::array{g, g, mEdge::zero(), mEdge::zero()}); + } else { + g = makeDDNode(static_cast(i), + std::array{g, mEdge::zero(), g, mEdge::zero()}); + } + } + } + incRef(g); + return g; + } + + Qubit lowerbound = 0; + for (auto i = 0U; i < garbage.size(); ++i) { + if (garbage[i]) { + lowerbound = static_cast(i); + break; + } + } + + auto g = CachedEdge{e.p, 1.}; + if (e.p->v >= lowerbound || normalizeWeights) { + g = reduceGarbageRecursion(e.p, garbage, lowerbound, regular, + normalizeWeights); + } + + for (std::size_t i = e.p->v + 1; i < garbage.size(); ++i) { + if (garbage[i]) { + if (regular) { + g = makeDDNode( + static_cast(i), + std::array{g, g, mCachedEdge::zero(), mCachedEdge::zero()}); + } else { + g = makeDDNode( + static_cast(i), + std::array{g, mCachedEdge::zero(), g, mCachedEdge::zero()}); + } + } + } + + auto weight = g.w * e.w; + if (normalizeWeights) { + weight = weight.mag(); + } + const auto res = mEdge{g.p, getCn().lookup(weight)}; + + incRef(res); + decRef(e); + return res; +} +mCachedEdge MatrixDDContainer::reduceAncillaeRecursion( + mNode* p, const std::vector& ancillary, const Qubit lowerbound, + const bool regular) { + if (p->v < lowerbound) { + return {p, 1.}; + } + + std::array edges{}; + std::bitset handled{}; + for (auto i = 0U; i < NEDGE; ++i) { + if (ancillary[p->v]) { + // no need to reduce ancillaries for entries that will be zeroed anyway + if ((i == 3) || (i == 1 && regular) || (i == 2 && !regular)) { + continue; + } + } + if (handled.test(i)) { + continue; + } + + if (p->e[i].isZeroTerminal()) { + edges[i] = {p->e[i].p, p->e[i].w}; + handled.set(i); + continue; + } + + if (p->e[i].isIdentity()) { + auto g = mCachedEdge::one(); + for (auto j = lowerbound; j < p->v; ++j) { + if (ancillary[j]) { + g = makeDDNode(j, + std::array{g, mCachedEdge::zero(), mCachedEdge::zero(), + mCachedEdge::zero()}); + } + } + edges[i] = {g.p, p->e[i].w}; + handled.set(i); + continue; + } + + edges[i] = + reduceAncillaeRecursion(p->e[i].p, ancillary, lowerbound, regular); + for (Qubit j = p->e[i].p->v + 1U; j < p->v; ++j) { + if (ancillary[j]) { + edges[i] = + makeDDNode(j, std::array{edges[i], mCachedEdge::zero(), + mCachedEdge::zero(), mCachedEdge::zero()}); + } + } + + for (auto j = i + 1U; j < NEDGE; ++j) { + if (p->e[i].p == p->e[j].p) { + edges[j] = edges[i]; + edges[j].w = edges[j].w * p->e[j].w; + handled.set(j); + } + } + edges[i].w = edges[i].w * p->e[i].w; + handled.set(i); + } + if (!ancillary[p->v]) { + return makeDDNode(p->v, edges); + } + + // something to reduce for this qubit + if (regular) { + return makeDDNode(p->v, std::array{edges[0], mCachedEdge::zero(), edges[2], + mCachedEdge::zero()}); + } + return makeDDNode(p->v, std::array{edges[0], edges[1], mCachedEdge::zero(), + mCachedEdge::zero()}); +} + +mCachedEdge MatrixDDContainer::reduceGarbageRecursion( + mNode* p, const std::vector& garbage, const Qubit lowerbound, + const bool regular, const bool normalizeWeights) { + if (!normalizeWeights && p->v < lowerbound) { + return {p, 1.}; + } + + std::array edges{}; + std::bitset handled{}; + for (auto i = 0U; i < NEDGE; ++i) { + if (handled.test(i)) { + continue; + } + + if (p->e[i].isZeroTerminal()) { + edges[i] = mCachedEdge::zero(); + handled.set(i); + continue; + } + + if (p->e[i].isIdentity()) { + edges[i] = mCachedEdge::one(); + for (auto j = lowerbound; j < p->v; ++j) { + if (garbage[j]) { + if (regular) { + edges[i] = makeDDNode(j, std::array{edges[i], edges[i], + mCachedEdge::zero(), + mCachedEdge::zero()}); + } else { + edges[i] = makeDDNode(j, std::array{edges[i], mCachedEdge::zero(), + edges[i], mCachedEdge::zero()}); + } + } + } + if (normalizeWeights) { + edges[i].w = edges[i].w * ComplexNumbers::mag(p->e[i].w); + } else { + edges[i].w = edges[i].w * p->e[i].w; + } + handled.set(i); + continue; + } + + edges[i] = reduceGarbageRecursion(p->e[i].p, garbage, lowerbound, regular, + normalizeWeights); + for (Qubit j = p->e[i].p->v + 1U; j < p->v; ++j) { + if (garbage[j]) { + if (regular) { + edges[i] = + makeDDNode(j, std::array{edges[i], edges[i], mCachedEdge::zero(), + mCachedEdge::zero()}); + } else { + edges[i] = makeDDNode(j, std::array{edges[i], mCachedEdge::zero(), + edges[i], mCachedEdge::zero()}); + } + } + } + + for (auto j = i + 1; j < NEDGE; ++j) { + if (p->e[i].p == p->e[j].p) { + edges[j] = edges[i]; + edges[j].w = edges[j].w * p->e[j].w; + if (normalizeWeights) { + edges[j].w = edges[j].w.mag(); + } + handled.set(j); + } + } + edges[i].w = edges[i].w * p->e[i].w; + if (normalizeWeights) { + edges[i].w = edges[i].w.mag(); + } + handled.set(i); + } + if (!garbage[p->v]) { + return makeDDNode(p->v, edges); + } + + if (regular) { + return makeDDNode(p->v, + std::array{addMagnitudes(edges[0], edges[2], p->v - 1), + addMagnitudes(edges[1], edges[3], p->v - 1), + mCachedEdge::zero(), mCachedEdge::zero()}); + } + return makeDDNode(p->v, + std::array{addMagnitudes(edges[0], edges[1], p->v - 1), + mCachedEdge::zero(), + addMagnitudes(edges[2], edges[3], p->v - 1), + mCachedEdge::zero()}); +} + +} // namespace dd diff --git a/src/dd/MemoryManager.cpp b/src/dd/MemoryManager.cpp index e0e188f36..fcb407df3 100644 --- a/src/dd/MemoryManager.cpp +++ b/src/dd/MemoryManager.cpp @@ -16,14 +16,13 @@ namespace dd { -MemoryManager::MemoryManager(size_t entrySize, - const std::size_t initialAllocationSize) - : entrySize_(entrySize), available(nullptr), - chunks(1, Chunk(initialAllocationSize * entrySize)), +MemoryManager::MemoryManager(const Config& config) + : entrySize_(config.entrySize), available(nullptr), + chunks(1, Chunk(config.initialAllocationSize * config.entrySize)), chunkIt(chunks[0].begin()), chunkEndIt(chunks[0].end()), - stats(entrySize) { + stats(config.entrySize) { stats.numAllocations = 1U; - stats.numAllocated = initialAllocationSize; + stats.numAllocated = config.initialAllocationSize; } LLBase* MemoryManager::get() { diff --git a/src/dd/NoiseFunctionality.cpp b/src/dd/NoiseFunctionality.cpp index 9d4974075..493a628b2 100644 --- a/src/dd/NoiseFunctionality.cpp +++ b/src/dd/NoiseFunctionality.cpp @@ -79,7 +79,7 @@ StochasticNoiseFunctionality::StochasticNoiseFunctionality( ampDampingFalseMulti( {1, 0, 0, oneMinusSqrtAmplitudeDampingProbabilityMulti}), noiseEffects(initializeNoiseEffects(cNoiseEffects)), - identityDD(Package::makeIdent()) { + identityDD(MatrixDDContainer::makeIdent()) { sanityCheckOfNoiseProbabilities(gateNoiseProbability, amplitudeDampingProb, multiQubitGateFactor); package->incRef(identityDD); @@ -145,7 +145,7 @@ mEdge StochasticNoiseFunctionality::stackOperation( op != nullptr) { return package->multiply(*op, operation); } - const auto gateDD = package->makeGateDD(matrix, target); + const auto gateDD = package->matrices().makeGateDD(matrix, target); package->stochasticNoiseOperationCache.insert(noiseOperation, target, gateDD); return package->multiply(gateDD, operation); } @@ -364,7 +364,7 @@ void DeterministicNoiseFunctionality::applyAmplitudeDampingToEdges( {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); + e[0] = package->add(e[0], {e[3].p, e[3].w * probability}, var); } else { e[0] = {e[3].p, e[3].w * probability}; } @@ -415,7 +415,7 @@ void DeterministicNoiseFunctionality::applyDepolarisationToEdges( } // e[0] = helperEdge[0] + helperEdge[1] - e[0] = package->add2(helperEdge[0], helperEdge[1], var); + e[0] = package->add(helperEdge[0], helperEdge[1], var); } // e[1]=(1-p)*e[1] @@ -444,7 +444,7 @@ void DeterministicNoiseFunctionality::applyDepolarisationToEdges( } else { helperEdge[1].w = 0; } - e[3] = package->add2(helperEdge[0], helperEdge[1], var); + e[3] = package->add(helperEdge[0], helperEdge[1], var); } } diff --git a/src/dd/Operations.cpp b/src/dd/Operations.cpp index 7b4d8a27b..de939fb14 100644 --- a/src/dd/Operations.cpp +++ b/src/dd/Operations.cpp @@ -43,16 +43,17 @@ MatrixDD getStandardOperationDD(Package& dd, const qc::OpType type, throw std::invalid_argument( "Expected exactly one target qubit for single-qubit gate"); } - return dd.makeGateDD(opToSingleQubitGateMatrix(type, params), controls, - targets[0U]); + return dd.matrices().makeGateDD(opToSingleQubitGateMatrix(type, params), + controls, targets[0U]); } if (qc::isTwoQubitGate(type)) { if (targets.size() != 2) { throw std::invalid_argument( "Expected two target qubits for two-qubit gate"); } - return dd.makeTwoQubitGateDD(opToTwoQubitGateMatrix(type, params), controls, - targets[0U], targets[1U]); + return dd.matrices().makeTwoQubitGateDD( + opToTwoQubitGateMatrix(type, params), controls, targets[0U], + targets[1U]); } throw std::runtime_error("Unexpected operation type"); } @@ -145,7 +146,7 @@ MatrixDD getDD(const qc::Operation& op, Package& dd, const auto type = op.getType(); if (type == qc::Barrier) { - return Package::makeIdent(); + return MatrixDDContainer::makeIdent(); } if (type == qc::GPhase) { @@ -153,7 +154,7 @@ MatrixDD getDD(const qc::Operation& op, Package& dd, if (inverse) { phase = -phase; } - auto id = Package::makeIdent(); + auto id = MatrixDDContainer::makeIdent(); id.w = dd.cn.lookup(std::cos(phase), std::sin(phase)); return id; } @@ -168,7 +169,7 @@ MatrixDD getDD(const qc::Operation& op, Package& dd, if (op.isCompoundOperation()) { const auto& compoundOp = dynamic_cast(op); - auto e = Package::makeIdent(); + auto e = MatrixDDContainer::makeIdent(); if (inverse) { for (const auto& operation : compoundOp) { e = dd.multiply(e, getInverseDD(*operation, dd, permutation)); diff --git a/src/dd/Package.cpp b/src/dd/Package.cpp index 79283409f..0700b6905 100644 --- a/src/dd/Package.cpp +++ b/src/dd/Package.cpp @@ -65,65 +65,47 @@ void Package::resize(const std::size_t nq) { "package with a wider Qubit type!"); } nqubits = nq; - vUniqueTable.resize(nqubits); - mUniqueTable.resize(nqubits); - dUniqueTable.resize(nqubits); + vectors().resize(nqubits); + matrices().resize(nqubits); + densities().resize(nqubits); stochasticNoiseOperationCache.resize(nqubits); } void Package::reset() { - clearUniqueTables(); - resetMemoryManagers(); - clearComputeTables(); -} - -void Package::resetMemoryManagers(const bool resizeToTotal) { - vMemoryManager.reset(resizeToTotal); - mMemoryManager.reset(resizeToTotal); - dMemoryManager.reset(resizeToTotal); - cMemoryManager.reset(resizeToTotal); -} + cUt.clear(); -void Package::clearUniqueTables() { - vUniqueTable.clear(); - mUniqueTable.clear(); - dUniqueTable.clear(); - cUniqueTable.clear(); + vectors().reset(); + matrices().reset(); + densities().reset(); } bool Package::garbageCollect(bool force) { // return immediately if no table needs collection - if (!force && !vUniqueTable.possiblyNeedsCollection() && - !mUniqueTable.possiblyNeedsCollection() && - !dUniqueTable.possiblyNeedsCollection() && - !cUniqueTable.possiblyNeedsCollection()) { + if (!force && !cUt.possiblyNeedsCollection() && + !matrices().possiblyNeedsCollection() && + !densities().possiblyNeedsCollection() && + !vectors().possiblyNeedsCollection()) { return false; } - const auto cCollect = cUniqueTable.garbageCollect(force); + const auto cCollect = cUt.garbageCollect(force); if (cCollect > 0) { // Collecting garbage in the complex numbers table requires collecting the // node tables as well force = true; } - const auto vCollect = vUniqueTable.garbageCollect(force); - const auto mCollect = mUniqueTable.garbageCollect(force); - const auto dCollect = dUniqueTable.garbageCollect(force); + const auto vCollect = vectors().garbageCollect(force); + const auto mCollect = matrices().garbageCollect(force); + const auto dCollect = densities().garbageCollect(force); // invalidate all compute tables involving vectors if any vector node has // been collected if (vCollect > 0) { - vectorAdd.clear(); - vectorInnerProduct.clear(); - vectorKronecker.clear(); matrixVectorMultiplication.clear(); } // invalidate all compute tables involving matrices if any matrix node has // been collected if (mCollect > 0) { - matrixAdd.clear(); - conjugateMatrixTranspose.clear(); - matrixKronecker.clear(); matrixTrace.clear(); matrixVectorMultiplication.clear(); matrixMatrixMultiplication.clear(); @@ -132,7 +114,6 @@ bool Package::garbageCollect(bool force) { // invalidate all compute tables involving density matrices if any density // matrix node has been collected if (dCollect > 0) { - densityAdd.clear(); densityDensityMultiplication.clear(); densityNoise.clear(); densityTrace.clear(); @@ -142,13 +123,8 @@ bool Package::garbageCollect(bool force) { if (cCollect > 0) { matrixVectorMultiplication.clear(); matrixMatrixMultiplication.clear(); - conjugateMatrixTranspose.clear(); - vectorInnerProduct.clear(); - vectorKronecker.clear(); - matrixKronecker.clear(); matrixTrace.clear(); stochasticNoiseOperationCache.clear(); - densityAdd.clear(); densityDensityMultiplication.clear(); densityNoise.clear(); densityTrace.clear(); @@ -156,616 +132,27 @@ bool Package::garbageCollect(bool force) { return vCollect > 0 || mCollect > 0 || cCollect > 0; } -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; -} - -vEdge Package::makeZeroState(const std::size_t n, const std::size_t start) { - if (n + start > nqubits) { - throw std::runtime_error{ - "Requested state with " + std::to_string(n + start) + - " qubits, but current package configuration only supports up to " + - std::to_string(nqubits) + - " qubits. Please allocate a larger package instance."}; - } - auto f = vEdge::one(); - for (std::size_t p = start; p < n + start; p++) { - f = makeDDNode(static_cast(p), std::array{f, vEdge::zero()}); - } - incRef(f); - return f; -} - -vEdge Package::makeBasisState(const std::size_t n, - const std::vector& state, - const std::size_t start) { - if (n + start > nqubits) { - throw std::runtime_error{ - "Requested state with " + std::to_string(n + start) + - " qubits, but current package configuration only supports up to " + - std::to_string(nqubits) + - " qubits. Please allocate a larger package instance."}; - } - auto f = vEdge::one(); - for (std::size_t p = start; p < n + start; ++p) { - if (!state[p]) { - f = makeDDNode(static_cast(p), std::array{f, vEdge::zero()}); - } else { - f = makeDDNode(static_cast(p), std::array{vEdge::zero(), f}); - } - } - incRef(f); - return f; -} -vEdge Package::makeBasisState(const std::size_t n, - const std::vector& state, - const std::size_t start) { - if (n + start > nqubits) { - throw std::runtime_error{ - "Requested state with " + std::to_string(n + start) + - " qubits, but current package configuration only supports up to " + - std::to_string(nqubits) + - " qubits. Please allocate a larger package instance."}; - } - if (state.size() < n) { - throw std::runtime_error("Insufficient qubit states provided. Requested " + - std::to_string(n) + ", but received " + - std::to_string(state.size())); - } - - auto f = vCachedEdge::one(); - for (std::size_t p = start; p < n + start; ++p) { - switch (state[p]) { - case BasisStates::zero: - f = makeDDNode(static_cast(p), std::array{f, vCachedEdge::zero()}); - break; - case BasisStates::one: - f = makeDDNode(static_cast(p), std::array{vCachedEdge::zero(), f}); - break; - case BasisStates::plus: - f = makeDDNode(static_cast(p), - std::array{ - {{f.p, dd::SQRT2_2}, {f.p, dd::SQRT2_2}}}); - break; - case BasisStates::minus: - f = makeDDNode(static_cast(p), - std::array{ - {{f.p, dd::SQRT2_2}, {f.p, -dd::SQRT2_2}}}); - break; - case BasisStates::right: - f = makeDDNode(static_cast(p), - std::array{ - {{f.p, dd::SQRT2_2}, {f.p, {0, dd::SQRT2_2}}}}); - break; - case BasisStates::left: - f = makeDDNode(static_cast(p), - std::array{ - {{f.p, dd::SQRT2_2}, {f.p, {0, -dd::SQRT2_2}}}}); - break; - } - } - const vEdge e{f.p, cn.lookup(f.w)}; - incRef(e); - return e; -} -vEdge Package::makeGHZState(const std::size_t n) { - if (n > nqubits) { - throw std::runtime_error{ - "Requested state with " + std::to_string(n) + - " qubits, but current package configuration only supports up to " + - std::to_string(nqubits) + - " qubits. Please allocate a larger package instance."}; - } - - if (n == 0U) { - return vEdge::one(); - } - - auto leftSubtree = vEdge::one(); - auto rightSubtree = vEdge::one(); - - for (std::size_t p = 0; p < n - 1; ++p) { - leftSubtree = makeDDNode(static_cast(p), - std::array{leftSubtree, vEdge::zero()}); - rightSubtree = makeDDNode(static_cast(p), - std::array{vEdge::zero(), rightSubtree}); - } - - const vEdge e = makeDDNode( - static_cast(n - 1), - std::array{ - {{leftSubtree.p, {&constants::sqrt2over2, &constants::zero}}, - {rightSubtree.p, {&constants::sqrt2over2, &constants::zero}}}}); - incRef(e); - return e; -} -vEdge Package::makeWState(const std::size_t n) { - if (n > nqubits) { - throw std::runtime_error{ - "Requested state with " + std::to_string(n) + - " qubits, but current package configuration only supports up to " + - std::to_string(nqubits) + - " qubits. Please allocate a larger package instance."}; - } - - if (n == 0U) { - return vEdge::one(); - } - - auto leftSubtree = vEdge::zero(); - if ((1. / sqrt(static_cast(n))) < RealNumber::eps) { - throw std::runtime_error( - "Requested qubit size for generating W-state would lead to an " - "underflow due to 1 / sqrt(n) being smaller than the currently set " - "tolerance " + - std::to_string(RealNumber::eps) + - ". If you still wanna run the computation, please lower " - "the tolerance accordingly."); - } +void Package::performCollapsingMeasurement(vEdge& rootEdge, const Qubit index, + const fp probability, + const bool measureZero) { + const GateMatrix measurementMatrix = + measureZero ? MEAS_ZERO_MAT : MEAS_ONE_MAT; - auto rightSubtree = vEdge::terminal(cn.lookup(1. / std::sqrt(n))); - for (size_t p = 0; p < n; ++p) { - leftSubtree = makeDDNode(static_cast(p), - std::array{leftSubtree, rightSubtree}); - if (p != n - 1U) { - rightSubtree = makeDDNode(static_cast(p), - std::array{rightSubtree, vEdge::zero()}); - } - } - incRef(leftSubtree); - return leftSubtree; -} -vEdge Package::makeStateFromVector(const CVec& stateVector) { - if (stateVector.empty()) { - return vEdge::one(); - } - const auto& length = stateVector.size(); - if ((length & (length - 1)) != 0) { - throw std::invalid_argument( - "State vector must have a length of a power of two."); - } + const auto measurementGate = matrices().makeGateDD(measurementMatrix, index); - if (length == 1) { - return vEdge::terminal(cn.lookup(stateVector[0])); - } + vEdge e = multiply(measurementGate, rootEdge); - const auto level = static_cast(std::log2(length) - 1); - const auto state = - makeStateFromVector(stateVector.begin(), stateVector.end(), level); - const vEdge e{state.p, cn.lookup(state.w)}; + assert(probability > 0.); + e.w = cn.lookup(e.w / std::sqrt(probability)); incRef(e); - return e; -} -mEdge Package::makeGateDD(const GateMatrix& mat, const qc::Qubit target) { - return makeGateDD(mat, qc::Controls{}, target); -} -mEdge Package::makeGateDD(const GateMatrix& mat, const qc::Control& control, - const qc::Qubit target) { - return makeGateDD(mat, qc::Controls{control}, target); -} -mEdge Package::makeGateDD(const GateMatrix& mat, const qc::Controls& controls, - const qc::Qubit target) { - if (std::any_of(controls.begin(), controls.end(), - [this](const auto& c) { - return c.qubit > static_cast(nqubits - 1U); - }) || - target > static_cast(nqubits - 1U)) { - throw std::runtime_error{ - "Requested gate acting on qubit(s) with index larger than " + - std::to_string(nqubits - 1U) + - " while the package configuration only supports up to " + - std::to_string(nqubits) + - " qubits. Please allocate a larger package instance."}; - } - std::array em{}; - for (auto i = 0U; i < NEDGE; ++i) { - em[i] = mCachedEdge::terminal(mat[i]); - } - - if (controls.empty()) { - // Single qubit operation - const auto e = makeDDNode(static_cast(target), em); - return {e.p, cn.lookup(e.w)}; - } - - auto it = controls.begin(); - auto edges = std::array{mCachedEdge::zero(), mCachedEdge::zero(), - mCachedEdge::zero(), mCachedEdge::zero()}; - - // process lines below target - for (; it != controls.end() && it->qubit < target; ++it) { - for (auto i1 = 0U; i1 < RADIX; ++i1) { - for (auto i2 = 0U; i2 < RADIX; ++i2) { - const auto i = (i1 * RADIX) + i2; - if (it->type == qc::Control::Type::Neg) { // neg. control - edges[0] = em[i]; - edges[3] = (i1 == i2) ? mCachedEdge::one() : mCachedEdge::zero(); - } else { // pos. control - edges[0] = (i1 == i2) ? mCachedEdge::one() : mCachedEdge::zero(); - edges[3] = em[i]; - } - em[i] = makeDDNode(static_cast(it->qubit), edges); - } - } - } - - // target line - auto e = makeDDNode(static_cast(target), em); - - // process lines above target - for (; it != controls.end(); ++it) { - if (it->type == qc::Control::Type::Neg) { // neg. control - edges[0] = e; - edges[3] = mCachedEdge::one(); - e = makeDDNode(static_cast(it->qubit), edges); - } else { // pos. control - edges[0] = mCachedEdge::one(); - edges[3] = e; - e = makeDDNode(static_cast(it->qubit), edges); - } - } - return {e.p, cn.lookup(e.w)}; -} -mEdge Package::makeTwoQubitGateDD(const TwoQubitGateMatrix& mat, - const qc::Qubit target0, - const qc::Qubit target1) { - return makeTwoQubitGateDD(mat, qc::Controls{}, target0, target1); -} -mEdge Package::makeTwoQubitGateDD(const TwoQubitGateMatrix& mat, - const qc::Control& control, - const qc::Qubit target0, - const qc::Qubit target1) { - return makeTwoQubitGateDD(mat, qc::Controls{control}, target0, target1); -} -mEdge Package::makeTwoQubitGateDD(const TwoQubitGateMatrix& mat, - const qc::Controls& controls, - const qc::Qubit target0, - const qc::Qubit target1) { - // sanity check - if (std::any_of(controls.begin(), controls.end(), - [this](const auto& c) { - return c.qubit > static_cast(nqubits - 1U); - }) || - target0 > static_cast(nqubits - 1U) || - target1 > static_cast(nqubits - 1U)) { - throw std::runtime_error{ - "Requested gate acting on qubit(s) with index larger than " + - std::to_string(nqubits - 1U) + - " while the package configuration only supports up to " + - std::to_string(nqubits) + - " qubits. Please allocate a larger package instance."}; - } - - // create terminal edge matrix - std::array, NEDGE> em{}; - for (auto i1 = 0U; i1 < NEDGE; i1++) { - const auto& matRow = mat.at(i1); - auto& emRow = em.at(i1); - for (auto i2 = 0U; i2 < NEDGE; i2++) { - emRow.at(i2) = mCachedEdge::terminal(matRow.at(i2)); - } - } - - // process lines below smaller target - auto it = controls.begin(); - const auto smallerTarget = std::min(target0, target1); - - auto edges = std::array{mCachedEdge::zero(), mCachedEdge::zero(), - mCachedEdge::zero(), mCachedEdge::zero()}; - - for (; it != controls.end() && it->qubit < smallerTarget; ++it) { - for (auto row = 0U; row < NEDGE; ++row) { - for (auto col = 0U; col < NEDGE; ++col) { - if (it->type == qc::Control::Type::Neg) { // negative control - edges[0] = em[row][col]; - edges[3] = (row == col) ? mCachedEdge::one() : mCachedEdge::zero(); - } else { // positive control - edges[0] = (row == col) ? mCachedEdge::one() : mCachedEdge::zero(); - edges[3] = em[row][col]; - } - em[row][col] = makeDDNode(static_cast(it->qubit), edges); - } - } - } - - // process the smaller target by taking the 16 submatrices and appropriately - // combining them into four DDs. - std::array em0{}; - for (std::size_t row = 0; row < RADIX; ++row) { - for (std::size_t col = 0; col < RADIX; ++col) { - std::array local{}; - if (target0 > target1) { - for (std::size_t i = 0; i < RADIX; ++i) { - for (std::size_t j = 0; j < RADIX; ++j) { - local.at((i * RADIX) + j) = - em.at((row * RADIX) + i).at((col * RADIX) + j); - } - } - } else { - for (std::size_t i = 0; i < RADIX; ++i) { - for (std::size_t j = 0; j < RADIX; ++j) { - local.at((i * RADIX) + j) = - em.at((i * RADIX) + row).at((j * RADIX) + col); - } - } - } - em0.at((row * RADIX) + col) = - makeDDNode(static_cast(smallerTarget), local); - } - } - - // process lines between the two targets - const auto largerTarget = std::max(target0, target1); - for (; it != controls.end() && it->qubit < largerTarget; ++it) { - for (auto i = 0U; i < NEDGE; ++i) { - if (it->type == qc::Control::Type::Neg) { // negative control - edges[0] = em0[i]; - edges[3] = - (i == 0 || i == 3) ? mCachedEdge::one() : mCachedEdge::zero(); - } else { // positive control - edges[0] = - (i == 0 || i == 3) ? mCachedEdge::one() : mCachedEdge::zero(); - edges[3] = em0[i]; - } - em0[i] = makeDDNode(static_cast(it->qubit), edges); - } - } - - // process the larger target by combining the four DDs from the smaller - // target - auto e = makeDDNode(static_cast(largerTarget), em0); - - // process lines above the larger target - for (; it != controls.end(); ++it) { - if (it->type == qc::Control::Type::Neg) { // negative control - edges[0] = e; - edges[3] = mCachedEdge::one(); - } else { // positive control - edges[0] = mCachedEdge::one(); - edges[3] = e; - } - e = makeDDNode(static_cast(it->qubit), edges); - } - - return {e.p, cn.lookup(e.w)}; -} -mEdge Package::makeDDFromMatrix(const CMat& matrix) { - if (matrix.empty()) { - return mEdge::one(); - } - - const auto& length = matrix.size(); - if ((length & (length - 1)) != 0) { - throw std::invalid_argument("Matrix must have a length of a power of two."); - } - - const auto& width = matrix[0].size(); - if (length != width) { - throw std::invalid_argument("Matrix must be square."); - } - - if (length == 1) { - return mEdge::terminal(cn.lookup(matrix[0][0])); - } - - const auto level = static_cast(std::log2(length) - 1); - const auto matrixDD = makeDDFromMatrix(matrix, level, 0, length, 0, width); - return {matrixDD.p, cn.lookup(matrixDD.w)}; -} -vCachedEdge Package::makeStateFromVector(const CVec::const_iterator& begin, - const CVec::const_iterator& end, - const Qubit level) { - if (level == 0U) { - assert(std::distance(begin, end) == 2); - const auto zeroSuccessor = vCachedEdge::terminal(*begin); - const auto oneSuccessor = vCachedEdge::terminal(*(begin + 1)); - return makeDDNode(0, {zeroSuccessor, oneSuccessor}); - } - - const auto half = std::distance(begin, end) / 2; - const auto zeroSuccessor = - makeStateFromVector(begin, begin + half, level - 1); - const auto oneSuccessor = makeStateFromVector(begin + half, end, level - 1); - return makeDDNode(level, {zeroSuccessor, oneSuccessor}); -} -mCachedEdge Package::makeDDFromMatrix(const CMat& matrix, const Qubit level, - const std::size_t rowStart, - const std::size_t rowEnd, - const std::size_t colStart, - const std::size_t colEnd) { - // base case - if (level == 0U) { - assert(rowEnd - rowStart == 2); - assert(colEnd - colStart == 2); - return makeDDNode( - 0U, {mCachedEdge::terminal(matrix[rowStart][colStart]), - mCachedEdge::terminal(matrix[rowStart][colStart + 1]), - mCachedEdge::terminal(matrix[rowStart + 1][colStart]), - mCachedEdge::terminal(matrix[rowStart + 1][colStart + 1])}); - } - - // recursively call the function on all quadrants - const auto rowMid = (rowStart + rowEnd) / 2; - const auto colMid = (colStart + colEnd) / 2; - const auto l = static_cast(level - 1U); - - return makeDDNode( - level, {makeDDFromMatrix(matrix, l, rowStart, rowMid, colStart, colMid), - makeDDFromMatrix(matrix, l, rowStart, rowMid, colMid, colEnd), - makeDDFromMatrix(matrix, l, rowMid, rowEnd, colStart, colMid), - makeDDFromMatrix(matrix, l, rowMid, rowEnd, colMid, colEnd)}); -} -void Package::clearComputeTables() { - vectorAdd.clear(); - matrixAdd.clear(); - vectorAddMagnitudes.clear(); - matrixAddMagnitudes.clear(); - conjugateVector.clear(); - conjugateMatrixTranspose.clear(); - matrixMatrixMultiplication.clear(); - matrixVectorMultiplication.clear(); - vectorInnerProduct.clear(); - 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) { - if (std::abs(ComplexNumbers::mag2(rootEdge.w) - 1.0) > epsilon) { - if (rootEdge.w.approximatelyZero()) { - throw std::runtime_error( - "Numerical instabilities led to a 0-vector! Abort simulation!"); - } - std::cerr << "WARNING in MAll: numerical instability occurred during " - "simulation: |alpha|^2 + |beta|^2 = " - << ComplexNumbers::mag2(rootEdge.w) << ", but should be 1!\n"; - } - - if (rootEdge.isTerminal()) { - return ""; - } - - vEdge cur = rootEdge; - const auto numberOfQubits = static_cast(rootEdge.p->v) + 1U; - - std::string result(numberOfQubits, '0'); - - std::uniform_real_distribution dist(0.0, 1.0L); - - for (auto i = numberOfQubits; i > 0; --i) { - fp p0 = ComplexNumbers::mag2(cur.p->e.at(0).w); - const fp p1 = ComplexNumbers::mag2(cur.p->e.at(1).w); - const fp tmp = p0 + p1; - - if (std::abs(tmp - 1.0) > epsilon) { - throw std::runtime_error("Added probabilities differ from 1 by " + - std::to_string(std::abs(tmp - 1.0))); - } - p0 /= tmp; - - const fp threshold = dist(mt); - if (threshold < p0) { - cur = cur.p->e.at(0); - } else { - result[cur.p->v] = '1'; - cur = cur.p->e.at(1); - } - } - - if (collapse) { - vEdge e = vEdge::one(); - std::array edges{}; - for (std::size_t p = 0U; p < numberOfQubits; ++p) { - if (result[p] == '0') { - edges[0] = e; - edges[1] = vEdge::zero(); - } else { - edges[0] = vEdge::zero(); - edges[1] = e; - } - e = makeDDNode(static_cast(p), edges); - } - incRef(e); - decRef(rootEdge); - rootEdge = e; - } - - return std::string{result.rbegin(), result.rend()}; -} -fp Package::assignProbabilities(const vEdge& edge, - std::unordered_map& probs) { - auto it = probs.find(edge.p); - if (it != probs.end()) { - return ComplexNumbers::mag2(edge.w) * it->second; - } - double sum{1}; - if (!edge.isTerminal()) { - sum = assignProbabilities(edge.p->e[0], probs) + - assignProbabilities(edge.p->e[1], probs); - } - - probs.insert({edge.p, sum}); - - return ComplexNumbers::mag2(edge.w) * sum; + decRef(rootEdge); + rootEdge = e; } -std::pair -Package::determineMeasurementProbabilities(const vEdge& rootEdge, - const Qubit index) { - std::map measurementProbabilities; - std::set visited; - std::queue q; - - measurementProbabilities[rootEdge.p] = ComplexNumbers::mag2(rootEdge.w); - visited.insert(rootEdge.p); - q.push(rootEdge.p); - - while (q.front()->v != index) { - const auto* ptr = q.front(); - q.pop(); - const fp prob = measurementProbabilities[ptr]; - - const auto& s0 = ptr->e[0]; - if (const auto s0w = static_cast(s0.w); - !s0w.approximatelyZero()) { - const fp tmp1 = prob * s0w.mag2(); - if (visited.find(s0.p) != visited.end()) { - measurementProbabilities[s0.p] = measurementProbabilities[s0.p] + tmp1; - } else { - measurementProbabilities[s0.p] = tmp1; - visited.insert(s0.p); - q.push(s0.p); - } - } - - const auto& s1 = ptr->e[1]; - if (const auto s1w = static_cast(s1.w); - !s1w.approximatelyZero()) { - const fp tmp1 = prob * s1w.mag2(); - if (visited.find(s1.p) != visited.end()) { - measurementProbabilities[s1.p] = measurementProbabilities[s1.p] + tmp1; - } else { - measurementProbabilities[s1.p] = tmp1; - visited.insert(s1.p); - q.push(s1.p); - } - } - } - - fp pzero{0}; - fp pone{0}; - while (!q.empty()) { - const auto* ptr = q.front(); - q.pop(); - const auto& s0 = ptr->e[0]; - if (const auto s0w = static_cast(s0.w); - !s0w.approximatelyZero()) { - pzero += measurementProbabilities[ptr] * s0w.mag2(); - } - const auto& s1 = ptr->e[1]; - if (const auto s1w = static_cast(s1.w); - !s1w.approximatelyZero()) { - pone += measurementProbabilities[ptr] * s1w.mag2(); - } - } - return {pzero, pone}; -} char Package::measureOneCollapsing(vEdge& rootEdge, const Qubit index, std::mt19937_64& mt, const fp epsilon) { const auto& [pzero, pone] = - determineMeasurementProbabilities(rootEdge, index); + vectors().determineMeasurementProbabilities(rootEdge, index); const fp sum = pzero + pone; if (std::abs(sum - 1) > epsilon) { throw std::runtime_error( @@ -782,6 +169,7 @@ 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'; @@ -789,17 +177,17 @@ char Package::measureOneCollapsing(dEdge& e, const Qubit index, const auto nrQubits = e.p->v + 1U; dEdge::setDensityMatrixTrue(e); - auto const measZeroDd = makeGateDD(MEAS_ZERO_MAT, index); + auto const measZeroDd = matrices().makeGateDD(MEAS_ZERO_MAT, index); - auto tmp0 = conjugateTranspose(measZeroDd); + auto tmp0 = matrices().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); + auto const measOneDd = matrices().makeGateDD(MEAS_ONE_MAT, index); + tmp0 = matrices().conjugateTranspose(measOneDd); tmp1 = multiply(e, densityFromMatrixEdge(tmp0), false); tmp2 = multiply(densityFromMatrixEdge(measOneDd), tmp1, true); measuredResult = '1'; @@ -819,95 +207,27 @@ char Package::measureOneCollapsing(dEdge& e, const Qubit index, cn.incRef(e.w); return measuredResult; } -void Package::performCollapsingMeasurement(vEdge& rootEdge, const Qubit index, - const fp probability, - const bool measureZero) { - const GateMatrix measurementMatrix = - measureZero ? MEAS_ZERO_MAT : MEAS_ONE_MAT; - - const auto measurementGate = makeGateDD(measurementMatrix, index); - - vEdge e = multiply(measurementGate, rootEdge); - - assert(probability > 0.); - e.w = cn.lookup(e.w / std::sqrt(probability)); - incRef(e); - decRef(rootEdge); - rootEdge = e; -} -vEdge Package::conjugate(const vEdge& a) { - const auto r = conjugateRec(a); - return {r.p, cn.lookup(r.w)}; -} -vCachedEdge Package::conjugateRec(const vEdge& a) { - if (a.isZeroTerminal()) { - return vCachedEdge::zero(); - } - - if (a.isTerminal()) { - return {a.p, ComplexNumbers::conj(a.w)}; - } - - if (const auto* r = conjugateVector.lookup(a.p); r != nullptr) { - return {r->p, r->w * ComplexNumbers::conj(a.w)}; - } - - std::array e{}; - e[0] = conjugateRec(a.p->e[0]); - e[1] = conjugateRec(a.p->e[1]); - auto res = makeDDNode(a.p->v, e); - conjugateVector.insert(a.p, res); - res.w = res.w * ComplexNumbers::conj(a.w); - return res; -} -mEdge Package::conjugateTranspose(const mEdge& a) { - const auto r = conjugateTransposeRec(a); - return {r.p, cn.lookup(r.w)}; -} -mCachedEdge Package::conjugateTransposeRec(const mEdge& a) { - if (a.isTerminal()) { // terminal case - return {a.p, ComplexNumbers::conj(a.w)}; - } - - // check if in compute table - if (const auto* r = conjugateMatrixTranspose.lookup(a.p); r != nullptr) { - return {r->p, r->w * ComplexNumbers::conj(a.w)}; - } - std::array e{}; - // conjugate transpose submatrices and rearrange as required - for (auto i = 0U; i < RADIX; ++i) { - for (auto j = 0U; j < RADIX; ++j) { - e[(RADIX * i) + j] = conjugateTransposeRec(a.p->e[(RADIX * j) + i]); - } - } - // create new top node - auto res = makeDDNode(a.p->v, e); - - // put it in the compute table - conjugateMatrixTranspose.insert(a.p, res); - - // adjust top weight including conjugate - return {res.p, res.w * ComplexNumbers::conj(a.w)}; -} -VectorDD Package::applyOperation(const MatrixDD& operation, const VectorDD& e) { - const auto tmp = multiply(operation, e); +mEdge Package::applyOperation(const mEdge& operation, const mEdge& e, + const bool applyFromLeft) { + const mEdge tmp = + applyFromLeft ? multiply(operation, e) : multiply(e, operation); incRef(tmp); decRef(e); garbageCollect(); return tmp; } -MatrixDD Package::applyOperation(const MatrixDD& operation, const MatrixDD& e, - const bool applyFromLeft) { - const MatrixDD tmp = - applyFromLeft ? multiply(operation, e) : multiply(e, operation); + +vEdge Package::applyOperation(const mEdge& operation, const vEdge& e) { + const auto tmp = multiply(operation, e); incRef(tmp); decRef(e); garbageCollect(); return tmp; } + dEdge Package::applyOperationToDensity(dEdge& e, const mEdge& operation) { - const auto tmp0 = conjugateTranspose(operation); + const auto tmp0 = matrices().conjugateTranspose(operation); const auto tmp1 = multiply(e, densityFromMatrixEdge(tmp0), false); const auto tmp2 = multiply(densityFromMatrixEdge(operation), tmp1, true); incRef(tmp2); @@ -917,108 +237,7 @@ dEdge Package::applyOperationToDensity(dEdge& e, const mEdge& operation) { 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 - return 0; - } - const auto w = std::max(x.p->v, y.p->v); - // Overall normalization factor needs to be conjugated - // before input into recursive private function - auto xCopy = vEdge{x.p, ComplexNumbers::conj(x.w)}; - return innerProduct(xCopy, y, w + 1U); -} -fp Package::fidelity(const vEdge& x, const vEdge& y) { - return innerProduct(x, y).mag2(); -} -fp Package::fidelityOfMeasurementOutcomes(const vEdge& e, - const SparsePVec& probs, - const qc::Permutation& permutation) { - if (e.w.approximatelyZero()) { - return 0.; - } - return fidelityOfMeasurementOutcomesRecursive(e, probs, 0, permutation, - e.p->v + 1U); -} -ComplexValue Package::innerProduct(const vEdge& x, const vEdge& y, - const Qubit var) { - const auto xWeight = static_cast(x.w); - if (xWeight.approximatelyZero()) { - return 0; - } - const auto yWeight = static_cast(y.w); - if (yWeight.approximatelyZero()) { - return 0; - } - - const auto rWeight = xWeight * yWeight; - if (var == 0) { // Multiplies terminal weights - return rWeight; - } - - if (const auto* r = vectorInnerProduct.lookup(x.p, y.p); r != nullptr) { - return r->w * rWeight; - } - - const auto w = static_cast(var - 1U); - ComplexValue sum = 0; - // Iterates through edge weights recursively until terminal - for (auto i = 0U; i < RADIX; i++) { - vEdge e1{}; - if (!x.isTerminal() && x.p->v == w) { - e1 = x.p->e[i]; - e1.w = ComplexNumbers::conj(e1.w); - } else { - e1 = {x.p, Complex::one()}; - } - vEdge e2{}; - if (!y.isTerminal() && y.p->v == w) { - e2 = y.p->e[i]; - } else { - e2 = {y.p, Complex::one()}; - } - sum += innerProduct(e1, e2, w); - } - vectorInnerProduct.insert(x.p, y.p, vCachedEdge::terminal(sum)); - return sum * rWeight; -} -fp Package::fidelityOfMeasurementOutcomesRecursive( - const vEdge& e, const SparsePVec& probs, const std::size_t i, - const qc::Permutation& permutation, const std::size_t nQubits) { - const auto top = ComplexNumbers::mag(e.w); - if (e.isTerminal()) { - auto idx = i; - if (!permutation.empty()) { - const auto binaryString = intToBinaryString(i, nQubits); - std::string filteredString(permutation.size(), '0'); - for (const auto& [physical, logical] : permutation) { - filteredString[logical] = binaryString[physical]; - } - idx = std::stoull(filteredString, nullptr, 2); - } - if (auto it = probs.find(idx); it != probs.end()) { - return top * std::sqrt(it->second); - } - return 0.; - } - - const std::size_t leftIdx = i; - fp leftContribution = 0.; - if (!e.p->e[0].w.approximatelyZero()) { - leftContribution = fidelityOfMeasurementOutcomesRecursive( - e.p->e[0], probs, leftIdx, permutation, nQubits); - } - - const std::size_t rightIdx = i | (1ULL << e.p->v); - auto rightContribution = 0.; - if (!e.p->e[1].w.approximatelyZero()) { - rightContribution = fidelityOfMeasurementOutcomesRecursive( - e.p->e[1], probs, rightIdx, permutation, nQubits); - } - - return top * (leftContribution + rightContribution); -} fp Package::expectationValue(const mEdge& x, const vEdge& y) { assert(!x.isZeroTerminal() && !y.isTerminal()); if (!x.isTerminal() && x.p->v > y.p->v) { @@ -1028,7 +247,7 @@ fp Package::expectationValue(const mEdge& x, const vEdge& y) { } const auto yPrime = multiply(x, y); - const ComplexValue expValue = innerProduct(y, yPrime); + const ComplexValue expValue = vectors().innerProduct(y, yPrime); assert(RealNumber::approximatelyZero(expValue.i)); return expValue.r; @@ -1038,440 +257,5 @@ mEdge Package::partialTrace(const mEdge& a, auto r = trace(a, eliminate, eliminate.size()); return {r.p, cn.lookup(r.w)}; } -bool Package::isCloseToIdentity(const mEdge& m, const fp tol, - const std::vector& garbage, - const bool checkCloseToOne) const { - std::unordered_set visited{}; - visited.reserve(mUniqueTable.getNumActiveEntries()); - return isCloseToIdentityRecursive(m, visited, tol, garbage, checkCloseToOne); -} -bool Package::isCloseToIdentityRecursive( - const mEdge& m, std::unordered_set& visited, const fp tol, - const std::vector& garbage, const bool checkCloseToOne) { - // immediately return if this node is identical to the identity or zero - if (m.isTerminal()) { - return true; - } - - // immediately return if this node has already been visited - if (visited.find(m.p) != visited.end()) { - return true; - } - - const auto n = m.p->v; - - if (garbage.size() > n && garbage[n]) { - return isCloseToIdentityRecursive(m.p->e[0U], visited, tol, garbage, - checkCloseToOne) && - isCloseToIdentityRecursive(m.p->e[1U], visited, tol, garbage, - checkCloseToOne) && - isCloseToIdentityRecursive(m.p->e[2U], visited, tol, garbage, - checkCloseToOne) && - isCloseToIdentityRecursive(m.p->e[3U], visited, tol, garbage, - checkCloseToOne); - } - - // check whether any of the middle successors is non-zero, i.e., m = [ x 0 0 - // y ] - const auto mag1 = dd::ComplexNumbers::mag2(m.p->e[1U].w); - const auto mag2 = dd::ComplexNumbers::mag2(m.p->e[2U].w); - if (mag1 > tol || mag2 > tol) { - return false; - } - if (checkCloseToOne) { - // check whether m = [ ~1 0 0 y ] - const auto mag0 = dd::ComplexNumbers::mag2(m.p->e[0U].w); - if (std::abs(mag0 - 1.0) > tol) { - return false; - } - const auto arg0 = dd::ComplexNumbers::arg(m.p->e[0U].w); - if (std::abs(arg0) > tol) { - return false; - } - - // check whether m = [ x 0 0 ~1 ] or m = [ x 0 0 ~0 ] (the last case is - // true for an ancillary qubit) - const auto mag3 = dd::ComplexNumbers::mag2(m.p->e[3U].w); - if (mag3 > tol) { - if (std::abs(mag3 - 1.0) > tol) { - return false; - } - const auto arg3 = dd::ComplexNumbers::arg(m.p->e[3U].w); - if (std::abs(arg3) > tol) { - return false; - } - } - } - // m either has the form [ ~1 0 0 ~1 ] or [ ~1 0 0 ~0 ] - const auto ident0 = isCloseToIdentityRecursive(m.p->e[0U], visited, tol, - garbage, checkCloseToOne); - - if (!ident0) { - return false; - } - // m either has the form [ I 0 0 ~1 ] or [ I 0 0 ~0 ] - const auto ident3 = isCloseToIdentityRecursive(m.p->e[3U], visited, tol, - garbage, checkCloseToOne); - - visited.insert(m.p); - return ident3; -} -mEdge Package::makeIdent() { return mEdge::one(); } -mEdge Package::createInitialMatrix(const std::vector& ancillary) { - return reduceAncillae(makeIdent(), ancillary); -} -mEdge Package::reduceAncillae(mEdge e, const std::vector& ancillary, - const bool regular) { - // return if no more ancillaries left - if (std::none_of(ancillary.begin(), ancillary.end(), - [](const bool v) { return v; }) || - e.isZeroTerminal()) { - return e; - } - - // if we have only identities and no other nodes - if (e.isIdentity()) { - auto g = e; - for (auto i = 0U; i < ancillary.size(); ++i) { - if (ancillary[i]) { - g = makeDDNode( - static_cast(i), - std::array{g, mEdge::zero(), mEdge::zero(), mEdge::zero()}); - } - } - incRef(g); - return g; - } - - Qubit lowerbound = 0; - for (auto i = 0U; i < ancillary.size(); ++i) { - if (ancillary[i]) { - lowerbound = static_cast(i); - break; - } - } - - auto g = CachedEdge{e.p, 1.}; - if (e.p->v >= lowerbound) { - g = reduceAncillaeRecursion(e.p, ancillary, lowerbound, regular); - } - - for (std::size_t i = e.p->v + 1; i < ancillary.size(); ++i) { - if (ancillary[i]) { - g = makeDDNode(static_cast(i), - std::array{g, mCachedEdge::zero(), mCachedEdge::zero(), - mCachedEdge::zero()}); - } - } - const auto res = mEdge{g.p, cn.lookup(g.w * e.w)}; - incRef(res); - decRef(e); - return res; -} -vEdge Package::reduceGarbage(vEdge& e, const std::vector& garbage, - const bool normalizeWeights) { - // return if no more garbage left - if (!normalizeWeights && - (std::none_of(garbage.begin(), garbage.end(), [](bool v) { return v; }) || - e.isTerminal())) { - return e; - } - Qubit lowerbound = 0; - for (std::size_t i = 0U; i < garbage.size(); ++i) { - if (garbage[i]) { - lowerbound = static_cast(i); - break; - } - } - if (!normalizeWeights && e.p->v < lowerbound) { - return e; - } - const auto f = - reduceGarbageRecursion(e.p, garbage, lowerbound, normalizeWeights); - auto weight = e.w * f.w; - if (normalizeWeights) { - weight = weight.mag(); - } - const auto res = vEdge{f.p, cn.lookup(weight)}; - incRef(res); - decRef(e); - return res; -} -mEdge Package::reduceGarbage(const mEdge& e, const std::vector& garbage, - const bool regular, const bool normalizeWeights) { - // return if no more garbage left - if (!normalizeWeights && - (std::none_of(garbage.begin(), garbage.end(), [](bool v) { return v; }) || - e.isZeroTerminal())) { - return e; - } - - // if we have only identities and no other nodes - if (e.isIdentity()) { - auto g = e; - for (auto i = 0U; i < garbage.size(); ++i) { - if (garbage[i]) { - if (regular) { - g = makeDDNode(static_cast(i), - std::array{g, g, mEdge::zero(), mEdge::zero()}); - } else { - g = makeDDNode(static_cast(i), - std::array{g, mEdge::zero(), g, mEdge::zero()}); - } - } - } - incRef(g); - return g; - } - - Qubit lowerbound = 0; - for (auto i = 0U; i < garbage.size(); ++i) { - if (garbage[i]) { - lowerbound = static_cast(i); - break; - } - } - - auto g = CachedEdge{e.p, 1.}; - if (e.p->v >= lowerbound || normalizeWeights) { - g = reduceGarbageRecursion(e.p, garbage, lowerbound, regular, - normalizeWeights); - } - - for (std::size_t i = e.p->v + 1; i < garbage.size(); ++i) { - if (garbage[i]) { - if (regular) { - g = makeDDNode( - static_cast(i), - std::array{g, g, mCachedEdge::zero(), mCachedEdge::zero()}); - } else { - g = makeDDNode( - static_cast(i), - std::array{g, mCachedEdge::zero(), g, mCachedEdge::zero()}); - } - } - } - - auto weight = g.w * e.w; - if (normalizeWeights) { - weight = weight.mag(); - } - const auto res = mEdge{g.p, cn.lookup(weight)}; - - incRef(res); - decRef(e); - return res; -} -mCachedEdge Package::reduceAncillaeRecursion(mNode* p, - const std::vector& ancillary, - const Qubit lowerbound, - const bool regular) { - if (p->v < lowerbound) { - return {p, 1.}; - } - - std::array edges{}; - std::bitset handled{}; - for (auto i = 0U; i < NEDGE; ++i) { - if (ancillary[p->v]) { - // no need to reduce ancillaries for entries that will be zeroed anyway - if ((i == 3) || (i == 1 && regular) || (i == 2 && !regular)) { - continue; - } - } - if (handled.test(i)) { - continue; - } - - if (p->e[i].isZeroTerminal()) { - edges[i] = {p->e[i].p, p->e[i].w}; - handled.set(i); - continue; - } - - if (p->e[i].isIdentity()) { - auto g = mCachedEdge::one(); - for (auto j = lowerbound; j < p->v; ++j) { - if (ancillary[j]) { - g = makeDDNode(j, - std::array{g, mCachedEdge::zero(), mCachedEdge::zero(), - mCachedEdge::zero()}); - } - } - edges[i] = {g.p, p->e[i].w}; - handled.set(i); - continue; - } - - edges[i] = - reduceAncillaeRecursion(p->e[i].p, ancillary, lowerbound, regular); - for (Qubit j = p->e[i].p->v + 1U; j < p->v; ++j) { - if (ancillary[j]) { - edges[i] = - makeDDNode(j, std::array{edges[i], mCachedEdge::zero(), - mCachedEdge::zero(), mCachedEdge::zero()}); - } - } - - for (auto j = i + 1U; j < NEDGE; ++j) { - if (p->e[i].p == p->e[j].p) { - edges[j] = edges[i]; - edges[j].w = edges[j].w * p->e[j].w; - handled.set(j); - } - } - edges[i].w = edges[i].w * p->e[i].w; - handled.set(i); - } - if (!ancillary[p->v]) { - return makeDDNode(p->v, edges); - } - - // something to reduce for this qubit - if (regular) { - return makeDDNode(p->v, std::array{edges[0], mCachedEdge::zero(), edges[2], - mCachedEdge::zero()}); - } - return makeDDNode(p->v, std::array{edges[0], edges[1], mCachedEdge::zero(), - mCachedEdge::zero()}); -} -vCachedEdge Package::reduceGarbageRecursion(vNode* p, - const std::vector& garbage, - const Qubit lowerbound, - const bool normalizeWeights) { - if (!normalizeWeights && p->v < lowerbound) { - return {p, 1.}; - } - - std::array edges{}; - std::bitset handled{}; - for (auto i = 0U; i < RADIX; ++i) { - if (!handled.test(i)) { - if (p->e[i].isTerminal()) { - const auto weight = normalizeWeights - ? ComplexNumbers::mag(p->e[i].w) - : static_cast(p->e[i].w); - edges[i] = {p->e[i].p, weight}; - } else { - edges[i] = reduceGarbageRecursion(p->e[i].p, garbage, lowerbound, - normalizeWeights); - for (auto j = i + 1; j < RADIX; ++j) { - if (p->e[i].p == p->e[j].p) { - edges[j] = edges[i]; - edges[j].w = edges[j].w * p->e[j].w; - if (normalizeWeights) { - edges[j].w = edges[j].w.mag(); - } - handled.set(j); - } - } - edges[i].w = edges[i].w * p->e[i].w; - if (normalizeWeights) { - edges[i].w = edges[i].w.mag(); - } - } - handled.set(i); - } - } - if (!garbage[p->v]) { - return makeDDNode(p->v, edges); - } - // something to reduce for this qubit - return makeDDNode(p->v, - std::array{addMagnitudes(edges[0], edges[1], p->v - 1), - vCachedEdge ::zero()}); -} -mCachedEdge Package::reduceGarbageRecursion(mNode* p, - const std::vector& garbage, - const Qubit lowerbound, - const bool regular, - const bool normalizeWeights) { - if (!normalizeWeights && p->v < lowerbound) { - return {p, 1.}; - } - - std::array edges{}; - std::bitset handled{}; - for (auto i = 0U; i < NEDGE; ++i) { - if (handled.test(i)) { - continue; - } - - if (p->e[i].isZeroTerminal()) { - edges[i] = mCachedEdge::zero(); - handled.set(i); - continue; - } - - if (p->e[i].isIdentity()) { - edges[i] = mCachedEdge::one(); - for (auto j = lowerbound; j < p->v; ++j) { - if (garbage[j]) { - if (regular) { - edges[i] = makeDDNode(j, std::array{edges[i], edges[i], - mCachedEdge::zero(), - mCachedEdge::zero()}); - } else { - edges[i] = makeDDNode(j, std::array{edges[i], mCachedEdge::zero(), - edges[i], mCachedEdge::zero()}); - } - } - } - if (normalizeWeights) { - edges[i].w = edges[i].w * ComplexNumbers::mag(p->e[i].w); - } else { - edges[i].w = edges[i].w * p->e[i].w; - } - handled.set(i); - continue; - } - - edges[i] = reduceGarbageRecursion(p->e[i].p, garbage, lowerbound, regular, - normalizeWeights); - for (Qubit j = p->e[i].p->v + 1U; j < p->v; ++j) { - if (garbage[j]) { - if (regular) { - edges[i] = - makeDDNode(j, std::array{edges[i], edges[i], mCachedEdge::zero(), - mCachedEdge::zero()}); - } else { - edges[i] = makeDDNode(j, std::array{edges[i], mCachedEdge::zero(), - edges[i], mCachedEdge::zero()}); - } - } - } - - for (auto j = i + 1; j < NEDGE; ++j) { - if (p->e[i].p == p->e[j].p) { - edges[j] = edges[i]; - edges[j].w = edges[j].w * p->e[j].w; - if (normalizeWeights) { - edges[j].w = edges[j].w.mag(); - } - handled.set(j); - } - } - edges[i].w = edges[i].w * p->e[i].w; - if (normalizeWeights) { - edges[i].w = edges[i].w.mag(); - } - handled.set(i); - } - if (!garbage[p->v]) { - return makeDDNode(p->v, edges); - } - - if (regular) { - return makeDDNode(p->v, - std::array{addMagnitudes(edges[0], edges[2], p->v - 1), - addMagnitudes(edges[1], edges[3], p->v - 1), - mCachedEdge::zero(), mCachedEdge::zero()}); - } - return makeDDNode(p->v, - std::array{addMagnitudes(edges[0], edges[1], p->v - 1), - mCachedEdge::zero(), - addMagnitudes(edges[2], edges[3], p->v - 1), - mCachedEdge::zero()}); -} } // namespace dd diff --git a/src/dd/Simulation.cpp b/src/dd/Simulation.cpp index 5e1d39639..a95ff059e 100644 --- a/src/dd/Simulation.cpp +++ b/src/dd/Simulation.cpp @@ -118,7 +118,7 @@ std::map sample(const qc::QuantumComputation& qc, std::map counts{}; for (std::size_t i = 0U; i < shots; ++i) { // measure all returns a string of the form "q(n-1) ... q(0)" - auto measurement = dd.measureAll(e, false, mt); + auto measurement = dd.vectors().measureAll(e, false, mt); counts.operator[](measurement) += 1U; } // reduce reference count of measured state @@ -246,6 +246,6 @@ std::map sample(const qc::QuantumComputation& qc, const std::size_t seed) { const auto nqubits = qc.getNqubits(); const auto dd = std::make_unique(nqubits); - return sample(qc, dd->makeZeroState(nqubits), *dd, shots, seed); + return sample(qc, dd->vectors().makeZeroState(nqubits), *dd, shots, seed); } } // namespace dd diff --git a/src/dd/VectorDDContainer.cpp b/src/dd/VectorDDContainer.cpp new file mode 100644 index 000000000..7fe6ccb1b --- /dev/null +++ b/src/dd/VectorDDContainer.cpp @@ -0,0 +1,572 @@ +/* + * Copyright (c) 2025 Chair for Design Automation, TUM + * All rights reserved. + * + * SPDX-License-Identifier: MIT + * + * Licensed under the MIT License + */ + +#include "dd/VectorDDContainer.hpp" + +#include "dd/GateMatrixDefinitions.hpp" + +#include + +namespace dd { + +vEdge VectorDDContainer::makeZeroState(const std::size_t n, + const std::size_t start) { + if (n + start > nqubits) { + throw std::runtime_error{ + "Requested state with " + std::to_string(n + start) + + " qubits, but current package configuration only supports up to " + + std::to_string(nqubits) + + " qubits. Please allocate a larger package instance."}; + } + auto f = vEdge::one(); + for (std::size_t p = start; p < n + start; p++) { + f = makeDDNode(static_cast(p), std::array{f, vEdge::zero()}); + } + incRef(f); + return f; +} + +vEdge VectorDDContainer::makeBasisState(const std::size_t n, + const std::vector& state, + const std::size_t start) { + if (n + start > nqubits) { + throw std::runtime_error{ + "Requested state with " + std::to_string(n + start) + + " qubits, but current package configuration only supports up to " + + std::to_string(nqubits) + + " qubits. Please allocate a larger package instance."}; + } + auto f = vEdge::one(); + for (std::size_t p = start; p < n + start; ++p) { + if (!state[p]) { + f = makeDDNode(static_cast(p), std::array{f, vEdge::zero()}); + } else { + f = makeDDNode(static_cast(p), std::array{vEdge::zero(), f}); + } + } + incRef(f); + return f; +} + +vEdge VectorDDContainer::makeBasisState(const std::size_t n, + const std::vector& state, + const std::size_t start) { + if (n + start > nqubits) { + throw std::runtime_error{ + "Requested state with " + std::to_string(n + start) + + " qubits, but current package configuration only supports up to " + + std::to_string(nqubits) + + " qubits. Please allocate a larger package instance."}; + } + if (state.size() < n) { + throw std::runtime_error("Insufficient qubit states provided. Requested " + + std::to_string(n) + ", but received " + + std::to_string(state.size())); + } + + auto f = vCachedEdge::one(); + for (std::size_t p = start; p < n + start; ++p) { + switch (state[p]) { + case BasisStates::zero: + f = makeDDNode(static_cast(p), std::array{f, vCachedEdge::zero()}); + break; + case BasisStates::one: + f = makeDDNode(static_cast(p), std::array{vCachedEdge::zero(), f}); + break; + case BasisStates::plus: + f = makeDDNode(static_cast(p), + std::array{ + {{f.p, dd::SQRT2_2}, {f.p, dd::SQRT2_2}}}); + break; + case BasisStates::minus: + f = makeDDNode(static_cast(p), + std::array{ + {{f.p, dd::SQRT2_2}, {f.p, -dd::SQRT2_2}}}); + break; + case BasisStates::right: + f = makeDDNode(static_cast(p), + std::array{ + {{f.p, dd::SQRT2_2}, {f.p, {0, dd::SQRT2_2}}}}); + break; + case BasisStates::left: + f = makeDDNode(static_cast(p), + std::array{ + {{f.p, dd::SQRT2_2}, {f.p, {0, -dd::SQRT2_2}}}}); + break; + } + } + const vEdge e{f.p, getCn().lookup(f.w)}; + incRef(e); + return e; +} + +vEdge VectorDDContainer::makeGHZState(const std::size_t n) { + if (n > nqubits) { + throw std::runtime_error{ + "Requested state with " + std::to_string(n) + + " qubits, but current package configuration only supports up to " + + std::to_string(nqubits) + + " qubits. Please allocate a larger package instance."}; + } + + if (n == 0U) { + return vEdge::one(); + } + + auto leftSubtree = vEdge::one(); + auto rightSubtree = vEdge::one(); + + for (std::size_t p = 0; p < n - 1; ++p) { + leftSubtree = makeDDNode(static_cast(p), + std::array{leftSubtree, vEdge::zero()}); + rightSubtree = makeDDNode(static_cast(p), + std::array{vEdge::zero(), rightSubtree}); + } + + const vEdge e = makeDDNode( + static_cast(n - 1), + std::array{ + {{leftSubtree.p, {&constants::sqrt2over2, &constants::zero}}, + {rightSubtree.p, {&constants::sqrt2over2, &constants::zero}}}}); + incRef(e); + return e; +} + +vEdge VectorDDContainer::makeWState(const std::size_t n) { + if (n > nqubits) { + throw std::runtime_error{ + "Requested state with " + std::to_string(n) + + " qubits, but current package configuration only supports up to " + + std::to_string(nqubits) + + " qubits. Please allocate a larger package instance."}; + } + + if (n == 0U) { + return vEdge::one(); + } + + auto leftSubtree = vEdge::zero(); + if ((1. / sqrt(static_cast(n))) < RealNumber::eps) { + throw std::runtime_error( + "Requested qubit size for generating W-state would lead to an " + "underflow due to 1 / sqrt(n) being smaller than the currently set " + "tolerance " + + std::to_string(RealNumber::eps) + + ". If you still wanna run the computation, please lower " + "the tolerance accordingly."); + } + + auto rightSubtree = vEdge::terminal(getCn().lookup(1. / std::sqrt(n))); + for (size_t p = 0; p < n; ++p) { + leftSubtree = makeDDNode(static_cast(p), + std::array{leftSubtree, rightSubtree}); + if (p != n - 1U) { + rightSubtree = makeDDNode(static_cast(p), + std::array{rightSubtree, vEdge::zero()}); + } + } + incRef(leftSubtree); + return leftSubtree; +} + +vEdge VectorDDContainer::makeStateFromVector(const CVec& stateVector) { + if (stateVector.empty()) { + return vEdge::one(); + } + const auto& length = stateVector.size(); + if ((length & (length - 1)) != 0) { + throw std::invalid_argument( + "State vector must have a length of a power of two."); + } + + if (length == 1) { + return vEdge::terminal(getCn().lookup(stateVector[0])); + } + + const auto level = static_cast(std::log2(length) - 1); + const auto state = + makeStateFromVector(stateVector.begin(), stateVector.end(), level); + const vEdge e{state.p, getCn().lookup(state.w)}; + incRef(e); + return e; +} + +vCachedEdge +VectorDDContainer::makeStateFromVector(const CVec::const_iterator& begin, + const CVec::const_iterator& end, + const Qubit level) { + if (level == 0U) { + assert(std::distance(begin, end) == 2); + const auto zeroSuccessor = vCachedEdge::terminal(*begin); + const auto oneSuccessor = vCachedEdge::terminal(*(begin + 1)); + return makeDDNode(0, {zeroSuccessor, oneSuccessor}); + } + + const auto half = std::distance(begin, end) / 2; + const auto zeroSuccessor = + makeStateFromVector(begin, begin + half, level - 1); + const auto oneSuccessor = makeStateFromVector(begin + half, end, level - 1); + return makeDDNode(level, {zeroSuccessor, oneSuccessor}); +} + +std::string VectorDDContainer::measureAll(vEdge& rootEdge, const bool collapse, + std::mt19937_64& mt, + const fp epsilon) { + if (std::abs(ComplexNumbers::mag2(rootEdge.w) - 1.0) > epsilon) { + if (rootEdge.w.approximatelyZero()) { + throw std::runtime_error( + "Numerical instabilities led to a 0-vector! Abort simulation!"); + } + std::cerr << "WARNING in MAll: numerical instability occurred during " + "simulation: |alpha|^2 + |beta|^2 = " + << ComplexNumbers::mag2(rootEdge.w) << ", but should be 1!\n"; + } + + if (rootEdge.isTerminal()) { + return ""; + } + + vEdge cur = rootEdge; + const auto numberOfQubits = static_cast(rootEdge.p->v) + 1U; + + std::string result(numberOfQubits, '0'); + + std::uniform_real_distribution dist(0.0, 1.0L); + + for (auto i = numberOfQubits; i > 0; --i) { + fp p0 = ComplexNumbers::mag2(cur.p->e.at(0).w); + const fp p1 = ComplexNumbers::mag2(cur.p->e.at(1).w); + const fp tmp = p0 + p1; + + if (std::abs(tmp - 1.0) > epsilon) { + throw std::runtime_error("Added probabilities differ from 1 by " + + std::to_string(std::abs(tmp - 1.0))); + } + p0 /= tmp; + + const fp threshold = dist(mt); + if (threshold < p0) { + cur = cur.p->e.at(0); + } else { + result[cur.p->v] = '1'; + cur = cur.p->e.at(1); + } + } + + if (collapse) { + vEdge e = vEdge::one(); + std::array edges{}; + for (std::size_t p = 0U; p < numberOfQubits; ++p) { + if (result[p] == '0') { + edges[0] = e; + edges[1] = vEdge::zero(); + } else { + edges[0] = vEdge::zero(); + edges[1] = e; + } + e = makeDDNode(static_cast(p), edges); + } + incRef(e); + decRef(rootEdge); + rootEdge = e; + } + + return std::string{result.rbegin(), result.rend()}; +} + +fp VectorDDContainer::assignProbabilities( + const vEdge& edge, std::unordered_map& probs) { + auto it = probs.find(edge.p); + if (it != probs.end()) { + return ComplexNumbers::mag2(edge.w) * it->second; + } + double sum{1}; + if (!edge.isTerminal()) { + sum = assignProbabilities(edge.p->e[0], probs) + + assignProbabilities(edge.p->e[1], probs); + } + + probs.insert({edge.p, sum}); + + return ComplexNumbers::mag2(edge.w) * sum; +} +std::pair +VectorDDContainer::determineMeasurementProbabilities(const vEdge& rootEdge, + const Qubit index) { + std::map measurementProbabilities; + std::set visited; + std::queue q; + + measurementProbabilities[rootEdge.p] = ComplexNumbers::mag2(rootEdge.w); + visited.insert(rootEdge.p); + q.push(rootEdge.p); + + while (q.front()->v != index) { + const auto* ptr = q.front(); + q.pop(); + const fp prob = measurementProbabilities[ptr]; + + const auto& s0 = ptr->e[0]; + if (const auto s0w = static_cast(s0.w); + !s0w.approximatelyZero()) { + const fp tmp1 = prob * s0w.mag2(); + if (visited.find(s0.p) != visited.end()) { + measurementProbabilities[s0.p] = measurementProbabilities[s0.p] + tmp1; + } else { + measurementProbabilities[s0.p] = tmp1; + visited.insert(s0.p); + q.push(s0.p); + } + } + + const auto& s1 = ptr->e[1]; + if (const auto s1w = static_cast(s1.w); + !s1w.approximatelyZero()) { + const fp tmp1 = prob * s1w.mag2(); + if (visited.find(s1.p) != visited.end()) { + measurementProbabilities[s1.p] = measurementProbabilities[s1.p] + tmp1; + } else { + measurementProbabilities[s1.p] = tmp1; + visited.insert(s1.p); + q.push(s1.p); + } + } + } + + fp pzero{0}; + fp pone{0}; + while (!q.empty()) { + const auto* ptr = q.front(); + q.pop(); + const auto& s0 = ptr->e[0]; + if (const auto s0w = static_cast(s0.w); + !s0w.approximatelyZero()) { + pzero += measurementProbabilities[ptr] * s0w.mag2(); + } + const auto& s1 = ptr->e[1]; + if (const auto s1w = static_cast(s1.w); + !s1w.approximatelyZero()) { + pone += measurementProbabilities[ptr] * s1w.mag2(); + } + } + + return {pzero, pone}; +} + +vEdge VectorDDContainer::conjugate(const vEdge& a) { + const auto r = conjugateRec(a); + return {r.p, getCn().lookup(r.w)}; +} + +vCachedEdge VectorDDContainer::conjugateRec(const vEdge& a) { + if (a.isZeroTerminal()) { + return vCachedEdge::zero(); + } + + if (a.isTerminal()) { + return {a.p, ComplexNumbers::conj(a.w)}; + } + + if (const auto* r = conjugateVector.lookup(a.p); r != nullptr) { + return {r->p, r->w * ComplexNumbers::conj(a.w)}; + } + + std::array e{}; + e[0] = conjugateRec(a.p->e[0]); + e[1] = conjugateRec(a.p->e[1]); + auto res = makeDDNode(a.p->v, e); + conjugateVector.insert(a.p, res); + res.w = res.w * ComplexNumbers::conj(a.w); + return res; +} + +ComplexValue VectorDDContainer::innerProduct(const vEdge& x, const vEdge& y) { + if (x.isTerminal() || y.isTerminal() || x.w.approximatelyZero() || + y.w.approximatelyZero()) { // the 0 case + return 0; + } + + const auto w = std::max(x.p->v, y.p->v); + // Overall normalization factor needs to be conjugated + // before input into recursive private function + auto xCopy = vEdge{x.p, ComplexNumbers::conj(x.w)}; + return innerProduct(xCopy, y, w + 1U); +} + +fp VectorDDContainer::fidelity(const vEdge& x, const vEdge& y) { + return innerProduct(x, y).mag2(); +} + +fp VectorDDContainer::fidelityOfMeasurementOutcomes( + const vEdge& e, const SparsePVec& probs, + const qc::Permutation& permutation) { + if (e.w.approximatelyZero()) { + return 0.; + } + return fidelityOfMeasurementOutcomesRecursive(e, probs, 0, permutation, + e.p->v + 1U); +} + +ComplexValue VectorDDContainer::innerProduct(const vEdge& x, const vEdge& y, + const Qubit var) { + const auto xWeight = static_cast(x.w); + if (xWeight.approximatelyZero()) { + return 0; + } + const auto yWeight = static_cast(y.w); + if (yWeight.approximatelyZero()) { + return 0; + } + + const auto rWeight = xWeight * yWeight; + if (var == 0) { // Multiplies terminal weights + return rWeight; + } + + if (const auto* r = vectorInnerProduct.lookup(x.p, y.p); r != nullptr) { + return r->w * rWeight; + } + + const auto w = static_cast(var - 1U); + ComplexValue sum = 0; + // Iterates through edge weights recursively until terminal + for (auto i = 0U; i < RADIX; i++) { + vEdge e1{}; + if (!x.isTerminal() && x.p->v == w) { + e1 = x.p->e[i]; + e1.w = ComplexNumbers::conj(e1.w); + } else { + e1 = {x.p, Complex::one()}; + } + vEdge e2{}; + if (!y.isTerminal() && y.p->v == w) { + e2 = y.p->e[i]; + } else { + e2 = {y.p, Complex::one()}; + } + sum += innerProduct(e1, e2, w); + } + vectorInnerProduct.insert(x.p, y.p, vCachedEdge::terminal(sum)); + return sum * rWeight; +} + +fp VectorDDContainer::fidelityOfMeasurementOutcomesRecursive( + const vEdge& e, const SparsePVec& probs, const std::size_t i, + const qc::Permutation& permutation, const std::size_t nQubits) { + const auto top = ComplexNumbers::mag(e.w); + if (e.isTerminal()) { + auto idx = i; + if (!permutation.empty()) { + const auto binaryString = intToBinaryString(i, nQubits); + std::string filteredString(permutation.size(), '0'); + for (const auto& [physical, logical] : permutation) { + filteredString[logical] = binaryString[physical]; + } + idx = std::stoull(filteredString, nullptr, 2); + } + if (auto it = probs.find(idx); it != probs.end()) { + return top * std::sqrt(it->second); + } + return 0.; + } + + const std::size_t leftIdx = i; + fp leftContribution = 0.; + if (!e.p->e[0].w.approximatelyZero()) { + leftContribution = fidelityOfMeasurementOutcomesRecursive( + e.p->e[0], probs, leftIdx, permutation, nQubits); + } + + const std::size_t rightIdx = i | (1ULL << e.p->v); + auto rightContribution = 0.; + if (!e.p->e[1].w.approximatelyZero()) { + rightContribution = fidelityOfMeasurementOutcomesRecursive( + e.p->e[1], probs, rightIdx, permutation, nQubits); + } + + return top * (leftContribution + rightContribution); +} + +vEdge VectorDDContainer::reduceGarbage(vEdge& e, + const std::vector& garbage, + const bool normalizeWeights) { + // return if no more garbage left + if (!normalizeWeights && + (std::none_of(garbage.begin(), garbage.end(), [](bool v) { return v; }) || + e.isTerminal())) { + return e; + } + Qubit lowerbound = 0; + for (std::size_t i = 0U; i < garbage.size(); ++i) { + if (garbage[i]) { + lowerbound = static_cast(i); + break; + } + } + if (!normalizeWeights && e.p->v < lowerbound) { + return e; + } + const auto f = + reduceGarbageRecursion(e.p, garbage, lowerbound, normalizeWeights); + auto weight = e.w * f.w; + if (normalizeWeights) { + weight = weight.mag(); + } + const auto res = vEdge{f.p, getCn().lookup(weight)}; + incRef(res); + decRef(e); + return res; +} + +vCachedEdge VectorDDContainer::reduceGarbageRecursion( + vNode* p, const std::vector& garbage, const Qubit lowerbound, + const bool normalizeWeights) { + if (!normalizeWeights && p->v < lowerbound) { + return {p, 1.}; + } + + std::array edges{}; + std::bitset handled{}; + for (auto i = 0U; i < RADIX; ++i) { + if (!handled.test(i)) { + if (p->e[i].isTerminal()) { + const auto weight = normalizeWeights + ? ComplexNumbers::mag(p->e[i].w) + : static_cast(p->e[i].w); + edges[i] = {p->e[i].p, weight}; + } else { + edges[i] = reduceGarbageRecursion(p->e[i].p, garbage, lowerbound, + normalizeWeights); + for (auto j = i + 1; j < RADIX; ++j) { + if (p->e[i].p == p->e[j].p) { + edges[j] = edges[i]; + edges[j].w = edges[j].w * p->e[j].w; + if (normalizeWeights) { + edges[j].w = edges[j].w.mag(); + } + handled.set(j); + } + } + edges[i].w = edges[i].w * p->e[i].w; + if (normalizeWeights) { + edges[i].w = edges[i].w.mag(); + } + } + handled.set(i); + } + } + if (!garbage[p->v]) { + return makeDDNode(p->v, edges); + } + // something to reduce for this qubit + return makeDDNode(p->v, + std::array{addMagnitudes(edges[0], edges[1], p->v - 1), + vCachedEdge ::zero()}); +} +} // namespace dd diff --git a/src/dd/statistics/PackageStatistics.cpp b/src/dd/statistics/PackageStatistics.cpp index 275ab70ef..a91a08a54 100644 --- a/src/dd/statistics/PackageStatistics.cpp +++ b/src/dd/statistics/PackageStatistics.cpp @@ -23,74 +23,31 @@ namespace dd { -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); -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(const Package& package) { - const auto vActiveEntries = - static_cast(package.vUniqueTable.getNumActiveEntries()); - const auto mActiveEntries = - static_cast(package.mUniqueTable.getNumActiveEntries()); - const auto dActiveEntries = - static_cast(package.dUniqueTable.getNumActiveEntries()); - - 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 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 vectorMem = package.vectors().computeActiveMemoryMiB(); + const auto matrixMem = package.matrices().computeActiveMemoryMiB(); + const auto densityMem = package.densities().computeActiveMemoryMiB(); const auto activeRealNumbers = - static_cast(package.cUniqueTable.getStats().numActiveEntries); + static_cast(package.cUt.getStats().numActiveEntries); const auto memoryForRealNumbers = activeRealNumbers * REAL_NUMBER_MEMORY_MIB; - return memoryForNodes + memoryForEdges + memoryForRealNumbers; + return vectorMem + matrixMem + densityMem + memoryForRealNumbers; } double computePeakMemoryMiB(const Package& package) { - const auto vPeakUsedEntries = - 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 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 vectorMem = package.vectors().computePeakMemoryMiB(); + const auto matrixMem = package.matrices().computePeakMemoryMiB(); + const auto densityMem = package.densities().computePeakMemoryMiB(); const auto peakRealNumbers = - static_cast(package.cMemoryManager.getStats().peakNumUsed); + static_cast(package.cMm.getStats().peakNumUsed); const auto memoryForRealNumbers = peakRealNumbers * REAL_NUMBER_MEMORY_MIB; - return memoryForNodes + memoryForEdges + memoryForRealNumbers; + return vectorMem + matrixMem + densityMem + memoryForRealNumbers; } nlohmann::basic_json<> getStatistics(const Package& package, @@ -100,44 +57,41 @@ nlohmann::basic_json<> getStatistics(const Package& package, j["data_structure"] = getDataStructureStatistics(); auto& vector = j["vector"]; - vector["unique_table"] = - package.vUniqueTable.getStatsJson(includeIndividualTables); - vector["memory_manager"] = package.vMemoryManager.getStats().json(); + package.vectors().addStatsJson(vector, includeIndividualTables); auto& matrix = j["matrix"]; - matrix["unique_table"] = - package.mUniqueTable.getStatsJson(includeIndividualTables); - matrix["memory_manager"] = package.mMemoryManager.getStats().json(); + package.matrices().addStatsJson(matrix, includeIndividualTables); auto& densityMatrix = j["density_matrix"]; - densityMatrix["unique_table"] = - package.dUniqueTable.getStatsJson(includeIndividualTables); - densityMatrix["memory_manager"] = package.dMemoryManager.getStats().json(); + package.densities().addStatsJson(densityMatrix, includeIndividualTables); auto& realNumbers = j["real_numbers"]; - realNumbers["unique_table"] = package.cUniqueTable.getStats().json(); - realNumbers["memory_manager"] = package.cMemoryManager.getStats().json(); - - 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(); + realNumbers["unique_table"] = package.cUt.getStats().json(); + realNumbers["memory_manager"] = package.cMm.getStats().json(); + + // TODO: auto& computeTables = j["compute_tables"]; + // TODO: computeTables["vector_add"] = package.vectorAdd.getStats().json(); + // TODO: computeTables["matrix_add"] = package.matrixAdd.getStats().json(); + // TODO: computeTables["density_matrix_add"] = + // package.densityAdd.getStats().json(); + // TODO: computeTables["matrix_conjugate_transpose"] = + // TODO: package.conjugateMatrixTranspose.getStats().json(); + // TODO: computeTables["matrix_vector_mult"] = + // TODO: package.matrixVectorMultiplication.getStats().json(); + // TODO: computeTables["matrix_matrix_mult"] = + // TODO: package.matrixMatrixMultiplication.getStats().json(); + // TODO: computeTables["density_density_mult"] = + // TODO: package.densityDensityMultiplication.getStats().json(); + // TODO: computeTables["vector_kronecker"] = + // package.vectorKronecker.getStats().json(); + // TODO: computeTables["matrix_kronecker"] = + // package.matrixKronecker.getStats().json(); + // TODO: computeTables["vector_inner_product"] = + // TODO: package.vectorInnerProduct.getStats().json(); + // TODO: computeTables["stochastic_noise_operations"] = + // TODO: package.stochasticNoiseOperationCache.getStats().json(); + // TODO: computeTables["density_noise_operations"] = + // TODO: package.densityNoise.getStats().json(); j["active_memory_mib"] = computeActiveMemoryMiB(package); j["peak_memory_mib"] = computePeakMemoryMiB(package); @@ -197,63 +151,71 @@ nlohmann::basic_json<> getDataStructureStatistics() { // Information about all the compute table entries // For every entry, we store the size in bytes and the alignment in bytes auto& ctEntries = j["ComplexTableEntries"]; - auto& vectorAdd = ctEntries["vector_add"]; - vectorAdd["size_B"] = sizeof(typename decltype(Package::vectorAdd)::Entry); - vectorAdd["alignment_B"] = - alignof(typename decltype(Package::vectorAdd)::Entry); - - auto& matrixAdd = ctEntries["matrix_add"]; - matrixAdd["size_B"] = sizeof(typename decltype(Package::matrixAdd)::Entry); - 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); - conjugateMatrixTranspose["alignment_B"] = - alignof(typename decltype(Package::conjugateMatrixTranspose)::Entry); - - auto& matrixVectorMult = ctEntries["matrix_vector_mult"]; - matrixVectorMult["size_B"] = - sizeof(typename decltype(Package::matrixVectorMultiplication)::Entry); - matrixVectorMult["alignment_B"] = - alignof(typename decltype(Package::matrixVectorMultiplication)::Entry); - - auto& matrixMatrixMult = ctEntries["matrix_matrix_mult"]; - matrixMatrixMult["size_B"] = - sizeof(typename decltype(Package::matrixMatrixMultiplication)::Entry); - 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); - vectorKronecker["alignment_B"] = - alignof(typename decltype(Package::vectorKronecker)::Entry); - - auto& matrixKronecker = ctEntries["matrix_kronecker"]; - matrixKronecker["size_B"] = - sizeof(typename decltype(Package::matrixKronecker)::Entry); - matrixKronecker["alignment_B"] = - alignof(typename decltype(Package::matrixKronecker)::Entry); - - auto& vectorInnerProduct = ctEntries["vector_inner_product"]; - vectorInnerProduct["size_B"] = - sizeof(typename decltype(Package::vectorInnerProduct)::Entry); - vectorInnerProduct["alignment_B"] = - alignof(typename decltype(Package::vectorInnerProduct)::Entry); - + // TODO: auto& vectorAdd = ctEntries["vector_add"]; + // TODO: vectorAdd["size_B"] = sizeof(typename + // decltype(Package::vectorAdd)::Entry); + // TODO: vectorAdd["alignment_B"] = + // TODO: alignof(typename decltype(Package::vectorAdd)::Entry); + // TODO: + // TODO: auto& matrixAdd = ctEntries["matrix_add"]; + // TODO: matrixAdd["size_B"] = sizeof(typename + // decltype(Package::matrixAdd)::Entry); + // TODO: matrixAdd["alignment_B"] = + // TODO: alignof(typename decltype(Package::matrixAdd)::Entry); + // TODO: + // TODO: auto& densityAdd = ctEntries["density_add"]; + // TODO: densityAdd["size_B"] = sizeof(typename + // decltype(Package::densityAdd)::Entry); + // TODO: densityAdd["alignment_B"] = + // TODO: alignof(typename decltype(Package::densityAdd)::Entry); + + // TODO: auto& conjugateMatrixTranspose = + // ctEntries["conjugate_matrix_transpose"]; + // TODO: conjugateMatrixTranspose["size_B"] = + // TODO: sizeof(typename + // decltype(Package::conjugateMatrixTranspose)::Entry); + // TODO: conjugateMatrixTranspose["alignment_B"] = + // TODO: alignof(typename + // decltype(Package::conjugateMatrixTranspose)::Entry); + + // auto& matrixVectorMult = ctEntries["matrix_vector_mult"]; + // matrixVectorMult["size_B"] = + // sizeof(typename decltype(Package::matrixVectorMultiplication)::Entry); + // matrixVectorMult["alignment_B"] = + // alignof(typename decltype(Package::matrixVectorMultiplication)::Entry); + // + // auto& matrixMatrixMult = ctEntries["matrix_matrix_mult"]; + // matrixMatrixMult["size_B"] = + // sizeof(typename decltype(Package::matrixMatrixMultiplication)::Entry); + // 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); + // vectorKronecker["alignment_B"] = + // alignof(typename decltype(Package::vectorKronecker)::Entry); + // + // auto& matrixKronecker = ctEntries["matrix_kronecker"]; + // matrixKronecker["size_B"] = + // sizeof(typename decltype(Package::matrixKronecker)::Entry); + // matrixKronecker["alignment_B"] = + // alignof(typename decltype(Package::matrixKronecker)::Entry); + // + // auto& vectorInnerProduct = ctEntries["vector_inner_product"]; + // vectorInnerProduct["size_B"] = + // sizeof(typename decltype(Package::vectorInnerProduct)::Entry); + // vectorInnerProduct["alignment_B"] = + // alignof(typename decltype(Package::vectorInnerProduct)::Entry); + // return j; } diff --git a/src/python/dd/register_dd.cpp b/src/python/dd/register_dd.cpp index fb22e7c8e..207093bc5 100644 --- a/src/python/dd/register_dd.cpp +++ b/src/python/dd/register_dd.cpp @@ -63,7 +63,7 @@ PYBIND11_MODULE(dd, mod, py::mod_gil_not_used()) { "simulate_statevector", [](const qc::QuantumComputation& qc) { auto dd = std::make_unique(qc.getNqubits()); - auto in = dd->makeZeroState(qc.getNqubits()); + auto in = dd->vectors().makeZeroState(qc.getNqubits()); const auto sim = dd::simulate(qc, in, *dd); return getVector(sim); }, diff --git a/src/python/dd/register_dd_package.cpp b/src/python/dd/register_dd_package.cpp index a37fff816..b2dcf937d 100644 --- a/src/python/dd/register_dd_package.cpp +++ b/src/python/dd/register_dd_package.cpp @@ -238,7 +238,9 @@ void registerDDPackage(const py::module& mod) { }, "vec"_a, "collapse"_a = false); - dd.def_static("identity", &dd::Package::makeIdent); + dd.def("identity", &dd::MatrixDDContainer::makeIdent, + // keep the DD package alive while the returned matrix DD is alive. + py::keep_alive<0, 1>()); using NumPyMatrix = py::array_t, py::array::c_style | py::array::forcecast>; diff --git a/test/algorithms/eval_dynamic_circuits.cpp b/test/algorithms/eval_dynamic_circuits.cpp index 4683a7044..4f0456097 100644 --- a/test/algorithms/eval_dynamic_circuits.cpp +++ b/test/algorithms/eval_dynamic_circuits.cpp @@ -126,7 +126,7 @@ TEST_P(DynamicCircuitEvalExactQPE, UnitaryTransformation) { iqpe.reorderOperations(); const auto finishedTransformation = std::chrono::steady_clock::now(); - dd::MatrixDD e = dd::Package::makeIdent(); + dd::MatrixDD e = dd::MatrixDDContainer::makeIdent(); dd->incRef(e); auto leftIt = qpe.begin(); @@ -295,7 +295,7 @@ TEST_P(DynamicCircuitEvalInexactQPE, UnitaryTransformation) { iqpe.reorderOperations(); const auto finishedTransformation = std::chrono::steady_clock::now(); - dd::MatrixDD e = dd::Package::makeIdent(); + dd::MatrixDD e = dd::MatrixDDContainer::makeIdent(); dd->incRef(e); auto leftIt = qpe.begin(); @@ -408,7 +408,7 @@ TEST_P(DynamicCircuitEvalBV, UnitaryTransformation) { dbv.reorderOperations(); const auto finishedTransformation = std::chrono::steady_clock::now(); - dd::MatrixDD e = dd::Package::makeIdent(); + dd::MatrixDD e = dd::MatrixDDContainer::makeIdent(); dd->incRef(e); auto leftIt = bv.begin(); @@ -517,7 +517,7 @@ TEST_P(DynamicCircuitEvalQFT, UnitaryTransformation) { dqft.reorderOperations(); const auto finishedTransformation = std::chrono::steady_clock::now(); - dd::MatrixDD e = dd::Package::makeIdent(); + dd::MatrixDD e = dd::MatrixDDContainer::makeIdent(); dd->incRef(e); auto leftIt = qft.begin(); diff --git a/test/algorithms/test_bernsteinvazirani.cpp b/test/algorithms/test_bernsteinvazirani.cpp index 0f0462ec8..e9804f133 100644 --- a/test/algorithms/test_bernsteinvazirani.cpp +++ b/test/algorithms/test_bernsteinvazirani.cpp @@ -122,7 +122,7 @@ TEST_P(BernsteinVazirani, DynamicEquivalenceSimulation) { qc::CircuitOptimizer::removeFinalMeasurements(bv); // simulate circuit - auto e = dd::simulate(bv, dd->makeZeroState(bv.getNqubits()), *dd); + auto e = dd::simulate(bv, dd->vectors().makeZeroState(bv.getNqubits()), *dd); // create dynamic BV circuit auto dbv = qc::createIterativeBernsteinVazirani(s); @@ -137,9 +137,10 @@ TEST_P(BernsteinVazirani, DynamicEquivalenceSimulation) { qc::CircuitOptimizer::removeFinalMeasurements(dbv); // simulate circuit - auto f = dd::simulate(dbv, dd->makeZeroState(dbv.getNqubits()), *dd); + auto f = + dd::simulate(dbv, dd->vectors().makeZeroState(dbv.getNqubits()), *dd); // calculate fidelity between both results - auto fidelity = dd->fidelity(e, f); + auto fidelity = dd->vectors().fidelity(e, f); EXPECT_NEAR(fidelity, 1.0, 1e-4); } diff --git a/test/algorithms/test_entanglement.cpp b/test/algorithms/test_entanglement.cpp index ca48ab671..9b05bc92e 100644 --- a/test/algorithms/test_entanglement.cpp +++ b/test/algorithms/test_entanglement.cpp @@ -45,23 +45,24 @@ TEST_P(Entanglement, FunctionTest) { const auto qc = qc::createGHZState(nq); const auto e = dd::buildFunctionality(qc, *dd); ASSERT_EQ(qc.getNops(), nq); - const auto r = dd->multiply(e, dd->makeZeroState(nq)); + const auto r = dd->multiply(e, dd->vectors().makeZeroState(nq)); ASSERT_EQ(r.getValueByPath(nq, std::string(nq, '0')), dd::SQRT2_2); ASSERT_EQ(r.getValueByPath(nq, std::string(nq, '1')), dd::SQRT2_2); } TEST_P(Entanglement, GHZRoutineFunctionTest) { const auto qc = qc::createGHZState(nq); - const auto e = dd::simulate(qc, dd->makeZeroState(nq), *dd); - const auto f = dd->makeGHZState(nq); + const auto e = dd::simulate(qc, dd->vectors().makeZeroState(nq), *dd); + const auto f = dd->vectors().makeGHZState(nq); EXPECT_EQ(e, f); } TEST(Entanglement, GHZStateEdgeCasesTest) { auto dd = std::make_unique(3); - EXPECT_EQ(dd->makeGHZState(0), - dd->makeBasisState(0, {dd::BasisStates::zero})); - EXPECT_EQ(dd->makeGHZState(0), dd->makeBasisState(0, {dd::BasisStates::one})); - ASSERT_THROW(dd->makeGHZState(6), std::runtime_error); + EXPECT_EQ(dd->vectors().makeGHZState(0), + dd->vectors().makeBasisState(0, {dd::BasisStates::zero})); + EXPECT_EQ(dd->vectors().makeGHZState(0), + dd->vectors().makeBasisState(0, {dd::BasisStates::one})); + ASSERT_THROW(dd->vectors().makeGHZState(6), std::runtime_error); } diff --git a/test/algorithms/test_qft.cpp b/test/algorithms/test_qft.cpp index 07dce1a51..a81b1215a 100644 --- a/test/algorithms/test_qft.cpp +++ b/test/algorithms/test_qft.cpp @@ -169,7 +169,7 @@ TEST_P(QFT, Simulation) { // there should be no error simulating the circuit ASSERT_NO_THROW({ - auto in = dd->makeZeroState(nqubits); + auto in = dd->vectors().makeZeroState(nqubits); sim = simulate(qc, in, *dd); }); qc.printStatistics(std::cout); diff --git a/test/algorithms/test_qpe.cpp b/test/algorithms/test_qpe.cpp index 42f735292..5ac3a8b73 100644 --- a/test/algorithms/test_qpe.cpp +++ b/test/algorithms/test_qpe.cpp @@ -135,8 +135,9 @@ TEST_P(QPE, QPETest) { ASSERT_NO_THROW({ qc::CircuitOptimizer::removeFinalMeasurements(qc); }); dd::VectorDD e{}; - ASSERT_NO_THROW( - { e = dd::simulate(qc, dd->makeZeroState(qc.getNqubits()), *dd); }); + ASSERT_NO_THROW({ + e = dd::simulate(qc, dd->vectors().makeZeroState(qc.getNqubits()), *dd); + }); // account for the eigenstate qubit by adding an offset const auto offset = 1ULL << (e.p->v + 1); @@ -218,7 +219,8 @@ TEST_P(QPE, DynamicEquivalenceSimulation) { qc::CircuitOptimizer::removeFinalMeasurements(qpe); // simulate circuit - auto e = dd::simulate(qpe, dd->makeZeroState(qpe.getNqubits()), *dd); + auto e = + dd::simulate(qpe, dd->vectors().makeZeroState(qpe.getNqubits()), *dd); // create standard IQPE circuit auto iqpe = qc::createIterativeQPE(lambda, precision); @@ -232,10 +234,11 @@ TEST_P(QPE, DynamicEquivalenceSimulation) { qc::CircuitOptimizer::removeFinalMeasurements(iqpe); // simulate circuit - auto f = dd::simulate(iqpe, dd->makeZeroState(iqpe.getNqubits()), *dd); + auto f = + dd::simulate(iqpe, dd->vectors().makeZeroState(iqpe.getNqubits()), *dd); // calculate fidelity between both results - auto fidelity = dd->fidelity(e, f); + auto fidelity = dd->vectors().fidelity(e, f); std::cout << "Fidelity of both circuits: " << fidelity << "\n"; EXPECT_NEAR(fidelity, 1.0, 1e-4); diff --git a/test/algorithms/test_random_clifford.cpp b/test/algorithms/test_random_clifford.cpp index f63c9aa23..90b965f53 100644 --- a/test/algorithms/test_random_clifford.cpp +++ b/test/algorithms/test_random_clifford.cpp @@ -43,7 +43,7 @@ TEST_P(RandomClifford, simulate) { for (size_t i = 0; i < numReps; ++i) { auto qc = qc::createRandomCliffordCircuit(nq, static_cast(nq) * nq); - auto in = dd->makeZeroState(nq); + auto in = dd->vectors().makeZeroState(nq); ASSERT_NO_THROW({ dd::simulate(qc, in, *dd); }); qc.printStatistics(std::cout); } diff --git a/test/algorithms/test_statepreparation.cpp b/test/algorithms/test_statepreparation.cpp index d3f994cf8..38a107f09 100644 --- a/test/algorithms/test_statepreparation.cpp +++ b/test/algorithms/test_statepreparation.cpp @@ -64,8 +64,9 @@ TEST_P(StatePreparation, StatePreparationCircuitSimulation) { ASSERT_NO_THROW({ qc = qc::createStatePreparationCircuit(amplitudes); }); auto dd = std::make_unique(qc.getNqubits()); dd::VectorDD e{}; - ASSERT_NO_THROW( - { e = dd::simulate(qc, dd->makeZeroState(qc.getNqubits()), *dd); }); + ASSERT_NO_THROW({ + e = dd::simulate(qc, dd->vectors().makeZeroState(qc.getNqubits()), *dd); + }); auto result = e.getVector(); for (size_t i = 0; i < expectedAmplitudes.size(); ++i) { ASSERT_NEAR(expectedAmplitudes[i].real(), result[i].real(), EPS); diff --git a/test/algorithms/test_wstate.cpp b/test/algorithms/test_wstate.cpp index 40396184a..36658b0d9 100644 --- a/test/algorithms/test_wstate.cpp +++ b/test/algorithms/test_wstate.cpp @@ -65,8 +65,8 @@ TEST_P(WState, RoutineFunctionTest) { const auto qc = qc::createWState(nq); const auto dd = std::make_unique(qc.getNqubits()); const dd::VectorDD e = - dd::simulate(qc, dd->makeZeroState(qc.getNqubits()), *dd); - const auto f = dd->makeWState(nq); + dd::simulate(qc, dd->vectors().makeZeroState(qc.getNqubits()), *dd); + const auto f = dd->vectors().makeWState(nq); EXPECT_EQ(e, f); } @@ -76,10 +76,13 @@ TEST(WState, WStateEdgeCasesTest) { const auto tolerance = dd::RealNumber::eps; dd::ComplexNumbers::setTolerance(0.1); - ASSERT_THROW(dd->makeWState(101), std::runtime_error); - EXPECT_EQ(dd->makeWState(0), dd->makeBasisState(0, {dd::BasisStates::zero})); - EXPECT_EQ(dd->makeWState(0), dd->makeBasisState(0, {dd::BasisStates::one})); - EXPECT_EQ(dd->makeWState(1), dd->makeBasisState(1, {dd::BasisStates::one})); - ASSERT_THROW(dd->makeWState(127), std::runtime_error); + ASSERT_THROW(dd->vectors().makeWState(101), std::runtime_error); + EXPECT_EQ(dd->vectors().makeWState(0), + dd->vectors().makeBasisState(0, {dd::BasisStates::zero})); + EXPECT_EQ(dd->vectors().makeWState(0), + dd->vectors().makeBasisState(0, {dd::BasisStates::one})); + EXPECT_EQ(dd->vectors().makeWState(1), + dd->vectors().makeBasisState(1, {dd::BasisStates::one})); + ASSERT_THROW(dd->vectors().makeWState(127), std::runtime_error); dd::ComplexNumbers::setTolerance(tolerance); } diff --git a/test/dd/test_dd_functionality.cpp b/test/dd/test_dd_functionality.cpp index 6d10635bc..f8958be0f 100644 --- a/test/dd/test_dd_functionality.cpp +++ b/test/dd/test_dd_functionality.cpp @@ -55,7 +55,7 @@ class DDFunctionality : public testing::TestWithParam { initialComplexCount = dd->cn.realCount(); // initial state preparation - e = ident = dd::Package::makeIdent(); + e = ident = dd::MatrixDDContainer::makeIdent(); dd->incRef(ident); std::array @@ -371,7 +371,8 @@ TEST_F(DDFunctionality, changePermutation) { "qreg q[2];" "x q[0];\n"; const auto qc = qasm3::Importer::imports(testfile); - const auto sim = simulate(qc, dd->makeZeroState(qc.getNqubits()), *dd); + const auto sim = + simulate(qc, dd->vectors().makeZeroState(qc.getNqubits()), *dd); EXPECT_TRUE(sim.p->e[0].isZeroTerminal()); EXPECT_TRUE(sim.p->e[1].w.exactlyOne()); EXPECT_TRUE(sim.p->e[1].p->e[1].isZeroTerminal()); @@ -485,8 +486,8 @@ TEST_F(DDFunctionality, classicControlledOperationConditions) { TEST_F(DDFunctionality, vectorKroneckerWithTerminal) { constexpr auto root = dd::vEdge::one(); - const auto zeroState = dd->makeZeroState(1); - const auto extendedRoot = dd->kronecker(zeroState, root, 0); + const auto zeroState = dd->vectors().makeZeroState(1); + const auto extendedRoot = dd->vectors().kronecker(zeroState, root, 0); EXPECT_EQ(zeroState, extendedRoot); } diff --git a/test/dd/test_dd_noise_functionality.cpp b/test/dd/test_dd_noise_functionality.cpp index 451a2ba6f..48dc3bf54 100644 --- a/test/dd/test_dd_noise_functionality.cpp +++ b/test/dd/test_dd_noise_functionality.cpp @@ -79,7 +79,7 @@ TEST_F(DDNoiseFunctionalityTest, DetSimulateAdder4TrackAPD) { auto dd = std::make_unique( qc.getNqubits(), dd::DENSITY_MATRIX_SIMULATOR_DD_PACKAGE_CONFIG); - auto rootEdge = dd->makeZeroDensityOperator(qc.getNqubits()); + auto rootEdge = dd->densities().makeZeroDensityOperator(qc.getNqubits()); const auto* const noiseEffects = "APDI"; @@ -111,7 +111,7 @@ TEST_F(DDNoiseFunctionalityTest, DetSimulateAdder4TrackD) { auto dd = std::make_unique( qc.getNqubits(), dd::DENSITY_MATRIX_SIMULATOR_DD_PACKAGE_CONFIG); - auto rootEdge = dd->makeZeroDensityOperator(qc.getNqubits()); + auto rootEdge = dd->densities().makeZeroDensityOperator(qc.getNqubits()); const auto* const noiseEffects = "D"; @@ -143,7 +143,7 @@ TEST_F(DDNoiseFunctionalityTest, testingMeasure) { auto dd = std::make_unique( qcOp.getNqubits(), dd::DENSITY_MATRIX_SIMULATOR_DD_PACKAGE_CONFIG); - auto rootEdge = dd->makeZeroDensityOperator(qcOp.getNqubits()); + auto rootEdge = dd->densities().makeZeroDensityOperator(qcOp.getNqubits()); auto deterministicNoiseFunctionality = dd::DeterministicNoiseFunctionality( *dd, qcOp.getNqubits(), 0.01, 0.02, 0.02, 0.04, {}); @@ -209,7 +209,7 @@ TEST_F(DDNoiseFunctionalityTest, StochSimulateAdder4TrackAPD) { *dd, qc.getNqubits(), 0.01, 0.02, 2., noiseEffects); for (size_t i = 0U; i < stochRuns; i++) { - auto rootEdge = dd->makeZeroState(qc.getNqubits()); + auto rootEdge = dd->vectors().makeZeroState(qc.getNqubits()); dd->incRef(rootEdge); for (auto const& op : qc) { @@ -262,7 +262,7 @@ TEST_F(DDNoiseFunctionalityTest, StochSimulateAdder4IdentityError) { *dd, qc.getNqubits(), 0.01, 0.02, 2., noiseEffects); for (size_t i = 0U; i < stochRuns; i++) { - auto rootEdge = dd->makeZeroState(qc.getNqubits()); + auto rootEdge = dd->vectors().makeZeroState(qc.getNqubits()); dd->incRef(rootEdge); for (auto const& op : qc) { diff --git a/test/dd/test_edge_functionality.cpp b/test/dd/test_edge_functionality.cpp index 8111f7370..34432e383 100644 --- a/test/dd/test_edge_functionality.cpp +++ b/test/dd/test_edge_functionality.cpp @@ -42,7 +42,7 @@ TEST(VectorFunctionality, GetValueByIndexEndianness) { auto dd = std::make_unique(2); const CVec state = {std::sqrt(0.1), std::sqrt(0.2), std::sqrt(0.3), std::sqrt(0.4)}; - const auto stateDD = dd->makeStateFromVector(state); + const auto stateDD = dd->vectors().makeStateFromVector(state); for (std::size_t i = 0U; i < state.size(); ++i) { EXPECT_EQ(state[i], stateDD.getValueByIndex(i)); @@ -58,7 +58,7 @@ TEST(VectorFunctionality, GetVectorRoundtrip) { auto dd = std::make_unique(2); const CVec state = {std::sqrt(0.1), std::sqrt(0.2), std::sqrt(0.3), std::sqrt(0.4)}; - const auto stateDD = dd->makeStateFromVector(state); + const auto stateDD = dd->vectors().makeStateFromVector(state); const auto stateVec = stateDD.getVector(); EXPECT_EQ(stateVec, state); } @@ -67,7 +67,7 @@ TEST(VectorFunctionality, GetVectorTolerance) { auto dd = std::make_unique(2); const CVec state = {std::sqrt(0.1), std::sqrt(0.2), std::sqrt(0.3), std::sqrt(0.4)}; - const auto stateDD = dd->makeStateFromVector(state); + const auto stateDD = dd->vectors().makeStateFromVector(state); const auto stateVec = stateDD.getVector(std::sqrt(0.1)); EXPECT_EQ(stateVec, state); const auto stateVec2 = stateDD.getVector(std::sqrt(0.1) + RealNumber::eps); @@ -86,7 +86,7 @@ TEST(VectorFunctionality, GetSparseVectorConsistency) { auto dd = std::make_unique(2); const CVec state = {std::sqrt(0.1), std::sqrt(0.2), std::sqrt(0.3), std::sqrt(0.4)}; - const auto stateDD = dd->makeStateFromVector(state); + const auto stateDD = dd->vectors().makeStateFromVector(state); const auto stateSparseVec = stateDD.getSparseVector(); const auto stateVec = stateDD.getVector(); for (const auto& [index, value] : stateSparseVec) { @@ -98,7 +98,7 @@ TEST(VectorFunctionality, GetSparseVectorTolerance) { auto dd = std::make_unique(2); const CVec state = {std::sqrt(0.1), std::sqrt(0.2), std::sqrt(0.3), std::sqrt(0.4)}; - const auto stateDD = dd->makeStateFromVector(state); + const auto stateDD = dd->vectors().makeStateFromVector(state); const auto stateSparseVec = stateDD.getSparseVector(std::sqrt(0.1)); for (const auto& [index, value] : stateSparseVec) { EXPECT_EQ(value, state[index]); @@ -124,7 +124,7 @@ TEST(VectorFunctionality, PrintVector) { auto dd = std::make_unique(2); const CVec state = {std::sqrt(0.1), std::sqrt(0.2), std::sqrt(0.3), std::sqrt(0.4)}; - const auto stateDD = dd->makeStateFromVector(state); + const auto stateDD = dd->vectors().makeStateFromVector(state); testing::internal::CaptureStdout(); stateDD.printVector(); const auto stateStr = testing::internal::GetCapturedStdout(); @@ -144,7 +144,7 @@ TEST(VectorFunctionality, AddToVector) { auto dd = std::make_unique(2); const CVec state = {std::sqrt(0.1), std::sqrt(0.2), std::sqrt(0.3), std::sqrt(0.4)}; - const auto stateDD = dd->makeStateFromVector(state); + const auto stateDD = dd->vectors().makeStateFromVector(state); stateDD.addToVector(vec); EXPECT_EQ(vec, state); } @@ -157,7 +157,7 @@ TEST(VectorFunctionality, SizeTerminal) { TEST(VectorFunctionality, SizeBellState) { auto dd = std::make_unique(2); const CVec state = {SQRT2_2, 0., 0., SQRT2_2}; - const auto bell = dd->makeStateFromVector(state); + const auto bell = dd->vectors().makeStateFromVector(state); EXPECT_EQ(bell.size(), 4); } @@ -185,7 +185,7 @@ TEST(MatrixFunctionality, GetValueByIndexEndianness) { {-std::sqrt(0.4), -std::sqrt(0.1), -std::sqrt(0.2), std::sqrt(0.3)}}; // clang-format on - const auto matDD = dd->makeDDFromMatrix(mat); + const auto matDD = dd->matrices().makeDDFromMatrix(mat); for (std::size_t i = 0U; i < mat.size(); ++i) { for (std::size_t j = 0U; j < mat.size(); ++j) { @@ -212,7 +212,7 @@ TEST(MatrixFunctionality, GetMatrixRoundtrip) { {-std::sqrt(0.4), -std::sqrt(0.1), -std::sqrt(0.2), std::sqrt(0.3)}}; // clang-format on - const auto matDD = dd->makeDDFromMatrix(mat); + const auto matDD = dd->matrices().makeDDFromMatrix(mat); const auto matVec = matDD.getMatrix(dd->qubits()); for (std::size_t i = 0U; i < mat.size(); ++i) { for (std::size_t j = 0U; j < mat.size(); ++j) { @@ -234,7 +234,7 @@ TEST(MatrixFunctionality, GetMatrixTolerance) { {-std::sqrt(0.4), -std::sqrt(0.1), -std::sqrt(0.2), std::sqrt(0.3)}}; // clang-format on - const auto matDD = dd->makeDDFromMatrix(mat); + const auto matDD = dd->matrices().makeDDFromMatrix(mat); const auto matVec = matDD.getMatrix(dd->qubits(), std::sqrt(0.1)); for (std::size_t i = 0U; i < mat.size(); ++i) { for (std::size_t j = 0U; j < mat.size(); ++j) { @@ -270,7 +270,7 @@ TEST(MatrixFunctionality, GetSparseMatrixConsistency) { {-std::sqrt(0.4), -std::sqrt(0.1), -std::sqrt(0.2), std::sqrt(0.3)}}; // clang-format on - const auto matDD = dd->makeDDFromMatrix(mat); + const auto matDD = dd->matrices().makeDDFromMatrix(mat); const auto matSparse = matDD.getSparseMatrix(dd->qubits()); const auto matDense = matDD.getMatrix(dd->qubits()); for (const auto& [index, value] : matSparse) { @@ -290,7 +290,7 @@ TEST(MatrixFunctionality, GetSparseMatrixTolerance) { {-std::sqrt(0.4), -std::sqrt(0.1), -std::sqrt(0.2), std::sqrt(0.3)}}; // clang-format on - const auto matDD = dd->makeDDFromMatrix(mat); + const auto matDD = dd->matrices().makeDDFromMatrix(mat); const auto matSparse = matDD.getSparseMatrix(dd->qubits(), std::sqrt(0.1)); const auto matDense = matDD.getMatrix(dd->qubits()); for (const auto& [index, value] : matSparse) { @@ -328,7 +328,7 @@ TEST(MatrixFunctionality, PrintMatrix) { {-std::sqrt(0.4), -std::sqrt(0.1), -std::sqrt(0.2), std::sqrt(0.3)}}; // clang-format on - const auto matDD = dd->makeDDFromMatrix(mat); + const auto matDD = dd->matrices().makeDDFromMatrix(mat); testing::internal::CaptureStdout(); matDD.printMatrix(dd->qubits()); const auto matStr = testing::internal::GetCapturedStdout(); @@ -353,7 +353,7 @@ TEST(MatrixFunctionality, SizeBellState) { {SQRT2_2, 0., 0., -SQRT2_2}}; // clang-format on - const auto bell = dd->makeDDFromMatrix(mat); + const auto bell = dd->matrices().makeDDFromMatrix(mat); EXPECT_EQ(bell.size(), 3); } @@ -374,7 +374,7 @@ TEST(DensityMatrixFunctionality, GetValueByIndexTerminal) { TEST(DensityMatrixFunctionality, GetValueByIndexProperDensityMatrix) { const auto nqubits = 1U; auto dd = std::make_unique(nqubits); - auto zero = dd->makeZeroDensityOperator(nqubits); + auto zero = dd->densities().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); @@ -408,7 +408,7 @@ TEST(DensityMatrixFunctionality, GetSparseMatrixTerminal) { TEST(DensityMatrixFunctionality, GetSparseMatrixConsistency) { const auto nqubits = 1U; auto dd = std::make_unique(nqubits); - auto zero = dd->makeZeroDensityOperator(nqubits); + auto zero = dd->densities().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); @@ -438,7 +438,7 @@ TEST(DensityMatrixFunctionality, PrintMatrixTerminal) { TEST(DensityMatrixFunctionality, PrintMatrix) { const auto nqubits = 1U; auto dd = std::make_unique(nqubits); - auto zero = dd->makeZeroDensityOperator(nqubits); + auto zero = dd->densities().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); diff --git a/test/dd/test_package.cpp b/test/dd/test_package.cpp index b36ad1a27..627ad252e 100644 --- a/test/dd/test_package.cpp +++ b/test/dd/test_package.cpp @@ -60,15 +60,15 @@ TEST(DDPackageTest, TrivialTest) { ASSERT_EQ(hGate.getValueByPath(1, "0"), SQRT2_2); - auto zeroState = dd->makeZeroState(1); + auto zeroState = dd->vectors().makeZeroState(1); auto hState = dd->multiply(hGate, zeroState); auto oneState = dd->multiply(xGate, zeroState); - ASSERT_EQ(dd->fidelity(zeroState, oneState), 0.0); + ASSERT_EQ(dd->vectors().fidelity(zeroState, oneState), 0.0); // repeat the same calculation — triggering compute table hit - ASSERT_EQ(dd->fidelity(zeroState, oneState), 0.0); - ASSERT_NEAR(dd->fidelity(zeroState, hState), 0.5, RealNumber::eps); - ASSERT_NEAR(dd->fidelity(oneState, hState), 0.5, RealNumber::eps); + ASSERT_EQ(dd->vectors().fidelity(zeroState, oneState), 0.0); + ASSERT_NEAR(dd->vectors().fidelity(zeroState, hState), 0.5, RealNumber::eps); + ASSERT_NEAR(dd->vectors().fidelity(oneState, hState), 0.5, RealNumber::eps); } TEST(DDPackageTest, BellState) { @@ -76,7 +76,7 @@ TEST(DDPackageTest, BellState) { auto hGate = getDD(qc::StandardOperation(1, qc::H), *dd); auto cxGate = getDD(qc::StandardOperation(1_pc, 0, qc::X), *dd); - auto zeroState = dd->makeZeroState(2); + auto zeroState = dd->vectors().makeZeroState(2); auto bellState = dd->multiply(dd->multiply(cxGate, hGate), zeroState); bellState.printVector(); @@ -98,7 +98,7 @@ TEST(DDPackageTest, BellState) { auto goalState = CVec{{SQRT2_2, 0.}, {0., 0.}, {0., 0.}, {SQRT2_2, 0.}}; ASSERT_EQ(bellState.getVector(), goalState); - ASSERT_DOUBLE_EQ(dd->fidelity(zeroState, bellState), 0.5); + ASSERT_DOUBLE_EQ(dd->vectors().fidelity(zeroState, bellState), 0.5); export2Dot(bellState, "bell_state_colored_labels.dot", true, true, false, false, false); @@ -159,7 +159,7 @@ TEST(DDPackageTest, QFTState) { qftOp = dd->multiply(h2Gate, qftOp); qftOp = dd->multiply(swapGate, qftOp); - auto qftState = dd->multiply(qftOp, dd->makeZeroState(3)); + auto qftState = dd->multiply(qftOp, dd->vectors().makeZeroState(3)); qftState.printVector(); @@ -270,18 +270,19 @@ TEST(DDPackageTest, CorruptedBellState) { auto hGate = getDD(qc::StandardOperation(1, qc::H), *dd); auto cxGate = getDD(qc::StandardOperation(1_pc, 0, qc::X), *dd); - auto zeroState = dd->makeZeroState(2); + auto zeroState = dd->vectors().makeZeroState(2); auto bellState = dd->multiply(dd->multiply(cxGate, hGate), zeroState); bellState.w = dd->cn.lookup(0.5, 0); // prints a warning std::mt19937_64 mt; // NOLINT(cert-msc51-cpp) - std::cout << dd->measureAll(bellState, false, mt) << "\n"; + std::cout << dd->vectors().measureAll(bellState, false, mt) << "\n"; bellState.w = Complex::zero(); - ASSERT_THROW(dd->measureAll(bellState, false, mt), std::runtime_error); + ASSERT_THROW(dd->vectors().measureAll(bellState, false, mt), + std::runtime_error); ASSERT_THROW(dd->measureOneCollapsing(bellState, 0, mt), std::runtime_error); } @@ -315,14 +316,14 @@ TEST(DDPackageTest, NegativeControl) { auto dd = std::make_unique(2); auto xGate = getDD(qc::StandardOperation(1_nc, 0, qc::X), *dd); - auto zeroState = dd->makeZeroState(2); + auto zeroState = dd->vectors().makeZeroState(2); auto state01 = dd->multiply(xGate, zeroState); EXPECT_EQ(state01.getValueByIndex(0b01).real(), 1.); } TEST(DDPackageTest, IdentityTrace) { auto dd = std::make_unique(4); - auto fullTrace = dd->trace(Package::makeIdent(), 4); + auto fullTrace = dd->trace(MatrixDDContainer::makeIdent(), 4); ASSERT_EQ(fullTrace.r, 1.); } @@ -330,14 +331,14 @@ TEST(DDPackageTest, IdentityTrace) { TEST(DDPackageTest, CNotKronTrace) { auto dd = std::make_unique(4); auto cxGate = getDD(qc::StandardOperation(1_pc, 0, qc::X), *dd); - auto cxGateKron = dd->kronecker(cxGate, cxGate, 2); + auto cxGateKron = dd->matrices().kronecker(cxGate, cxGate, 2); auto fullTrace = dd->trace(cxGateKron, 4); ASSERT_EQ(fullTrace, 0.25); } TEST(DDPackageTest, PartialIdentityTrace) { auto dd = std::make_unique(2); - auto tr = dd->partialTrace(Package::makeIdent(), {false, true}); + auto tr = dd->partialTrace(MatrixDDContainer::makeIdent(), {false, true}); auto mul = dd->multiply(tr, tr); EXPECT_EQ(RealNumber::val(mul.w.r), 1.); } @@ -367,7 +368,7 @@ TEST(DDPackageTest, PartialTraceKeepInnerQubits) { getDD(qc::StandardOperation(qc::Targets{0, 1}, qc::SWAP), *dd); auto swapKron = swapGate; for (std::size_t i = 0; i < 3; ++i) { - swapKron = dd->kronecker(swapKron, swapGate, 2); + swapKron = dd->matrices().kronecker(swapKron, swapGate, 2); } auto fullTraceOriginal = dd->trace(swapKron, numQubits); auto ptr = dd->partialTrace( @@ -389,7 +390,7 @@ TEST(DDPackageTest, TraceComplexity) { const auto hGate = getDD(qc::StandardOperation(0, qc::H), *dd); auto hKron = hGate; for (std::size_t i = 0; i < numQubits - 1; ++i) { - hKron = dd->kronecker(hKron, hGate, 1); + hKron = dd->matrices().kronecker(hKron, hGate, 1); } dd->trace(hKron, numQubits); const auto& stats = computeTable.getStats(); @@ -404,11 +405,11 @@ TEST(DDPackageTest, KeepBottomQubitsPartialTraceComplexity) { // recurse further but immediately returns the current CachedEdge. constexpr std::size_t numQubits = 8; auto dd = std::make_unique(numQubits); - auto& uniqueTable = dd->getUniqueTable(); + auto& uniqueTable = dd->matrices().getUniqueTable(); const auto hGate = getDD(qc::StandardOperation(0, qc::H), *dd); auto hKron = hGate; for (std::size_t i = 0; i < numQubits - 1; ++i) { - hKron = dd->kronecker(hKron, hGate, 1); + hKron = dd->matrices().kronecker(hKron, hGate, 1); } constexpr std::size_t maxNodeVal = 6; @@ -433,13 +434,13 @@ TEST(DDPackageTest, PartialTraceComplexity) { // bottom qubits. constexpr std::size_t numQubits = 9; auto dd = std::make_unique(numQubits); - auto& uniqueTable = dd->getUniqueTable(); + const auto& uniqueTable = dd->matrices().getUniqueTable(); const auto hGate = getDD(qc::StandardOperation(0, qc::H), *dd); auto hKron = hGate; for (std::size_t i = 0; i < numQubits - 2; ++i) { - hKron = dd->kronecker(hKron, hGate, 1); + hKron = dd->matrices().kronecker(hKron, hGate, 1); } - hKron = dd->kronecker(hKron, Package::makeIdent(), 1); + hKron = dd->matrices().kronecker(hKron, MatrixDDContainer::makeIdent(), 1); constexpr std::size_t maxNodeVal = 6; std::array lookupValues{}; @@ -464,11 +465,11 @@ TEST(DDPackageTest, StateGenerationManipulation) { auto dd = std::make_unique(nqubits); auto b = std::vector(nqubits, false); b[0] = b[1] = true; - auto e = dd->makeBasisState(nqubits, b); - auto f = dd->makeBasisState(nqubits, {BasisStates::zero, BasisStates::one, - BasisStates::plus, BasisStates::minus, - BasisStates::left, BasisStates::right}); - dd->vUniqueTable.print(); + auto e = dd->vectors().makeBasisState(nqubits, b); + auto f = dd->vectors().makeBasisState( + nqubits, {BasisStates::zero, BasisStates::one, BasisStates::plus, + BasisStates::minus, BasisStates::left, BasisStates::right}); + dd->vectors().getUniqueTable().print(); dd->decRef(e); dd->decRef(f); } @@ -478,17 +479,19 @@ TEST(DDPackageTest, VectorSerializationTest) { auto hGate = getDD(qc::StandardOperation(1, qc::H), *dd); auto cxGate = getDD(qc::StandardOperation(1_pc, 0, qc::X), *dd); - auto zeroState = dd->makeZeroState(2); + auto zeroState = dd->vectors().makeZeroState(2); auto bellState = dd->multiply(dd->multiply(cxGate, hGate), zeroState); serialize(bellState, "bell_state.dd", false); - auto deserializedBellState = dd->deserialize("bell_state.dd", false); + auto deserializedBellState = + dd->vectors().deserialize("bell_state.dd", false); EXPECT_EQ(bellState, deserializedBellState); std::filesystem::remove("bell_state.dd"); serialize(bellState, "bell_state_binary.dd", true); - deserializedBellState = dd->deserialize("bell_state_binary.dd", true); + deserializedBellState = + dd->vectors().deserialize("bell_state_binary.dd", true); EXPECT_EQ(bellState, deserializedBellState); std::filesystem::remove("bell_state_binary.dd"); } @@ -584,13 +587,14 @@ TEST(DDPackageTest, MatrixSerializationTest) { auto bellMatrix = dd->multiply(cxGate, hGate); serialize(bellMatrix, "bell_matrix.dd", false); - auto deserializedBellMatrix = dd->deserialize("bell_matrix.dd", false); + auto deserializedBellMatrix = + dd->matrices().deserialize("bell_matrix.dd", false); EXPECT_EQ(bellMatrix, deserializedBellMatrix); std::filesystem::remove("bell_matrix.dd"); serialize(bellMatrix, "bell_matrix_binary.dd", true); deserializedBellMatrix = - dd->deserialize("bell_matrix_binary.dd", true); + dd->matrices().deserialize("bell_matrix_binary.dd", true); EXPECT_EQ(bellMatrix, deserializedBellMatrix); std::filesystem::remove("bell_matrix_binary.dd"); } @@ -600,50 +604,75 @@ TEST(DDPackageTest, SerializationErrors) { auto hGate = getDD(qc::StandardOperation(1, qc::H), *dd); auto cxGate = getDD(qc::StandardOperation(1_pc, 0, qc::X), *dd); - auto zeroState = dd->makeZeroState(2); + auto zeroState = dd->vectors().makeZeroState(2); auto bellState = dd->multiply(dd->multiply(cxGate, hGate), zeroState); // test non-existing file EXPECT_THROW(serialize(bellState, "./path/that/does/not/exist/filename.dd"), std::invalid_argument); EXPECT_THROW( - dd->deserialize("./path/that/does/not/exist/filename.dd", true), + { + [[maybe_unused]] const auto& v = dd->vectors().deserialize( + "./path/that/does/not/exist/filename.dd", true); + }, std::invalid_argument); // test wrong version number std::stringstream ss{}; ss << 2 << "\n"; - EXPECT_THROW(dd->deserialize(ss, false), std::runtime_error); + EXPECT_THROW( + { + [[maybe_unused]] const auto& v = dd->vectors().deserialize(ss, false); + }, + std::runtime_error); ss << 2 << "\n"; - EXPECT_THROW(dd->deserialize(ss, false), std::runtime_error); + EXPECT_THROW( + { + [[maybe_unused]] const auto& v = dd->matrices().deserialize(ss, false); + }, + std::runtime_error); ss.str(""); std::remove_const_t version = 2; ss.write(reinterpret_cast(&version), sizeof(decltype(SERIALIZATION_VERSION))); - EXPECT_THROW(dd->deserialize(ss, true), std::runtime_error); + EXPECT_THROW( + { [[maybe_unused]] const auto& v = dd->vectors().deserialize(ss, true); }, + std::runtime_error); ss.write(reinterpret_cast(&version), sizeof(decltype(SERIALIZATION_VERSION))); - EXPECT_THROW(dd->deserialize(ss, true), std::runtime_error); + EXPECT_THROW( + { + [[maybe_unused]] const auto& v = dd->matrices().deserialize(ss, true); + }, + std::runtime_error); // test wrong format ss.str(""); ss << "1\n"; ss << "not_complex\n"; - EXPECT_THROW(dd->deserialize(ss), std::runtime_error); + EXPECT_THROW( + { [[maybe_unused]] const auto& v = dd->vectors().deserialize(ss); }, + std::runtime_error); ss << "1\n"; ss << "not_complex\n"; - EXPECT_THROW(dd->deserialize(ss), std::runtime_error); + EXPECT_THROW( + { [[maybe_unused]] const auto& v = dd->matrices().deserialize(ss); }, + std::runtime_error); ss.str(""); ss << "1\n"; ss << "1.0\n"; ss << "no_node_here\n"; - EXPECT_THROW(dd->deserialize(ss), std::runtime_error); + EXPECT_THROW( + { [[maybe_unused]] const auto& v = dd->vectors().deserialize(ss); }, + std::runtime_error); ss << "1\n"; ss << "1.0\n"; ss << "no_node_here\n"; - EXPECT_THROW(dd->deserialize(ss), std::runtime_error); + EXPECT_THROW( + { [[maybe_unused]] const auto& v = dd->matrices().deserialize(ss); }, + std::runtime_error); } TEST(DDPackageTest, Ancillaries) { @@ -654,12 +683,12 @@ TEST(DDPackageTest, Ancillaries) { dd->incRef(bellMatrix); auto reducedBellMatrix = - dd->reduceAncillae(bellMatrix, {false, false, false, false}); + dd->matrices().reduceAncillae(bellMatrix, {false, false, false, false}); EXPECT_EQ(bellMatrix, reducedBellMatrix); dd->incRef(bellMatrix); reducedBellMatrix = - dd->reduceAncillae(bellMatrix, {false, false, true, true}); + dd->matrices().reduceAncillae(bellMatrix, {false, false, true, true}); EXPECT_TRUE(reducedBellMatrix.p->e[1].isZeroTerminal()); EXPECT_TRUE(reducedBellMatrix.p->e[2].isZeroTerminal()); EXPECT_TRUE(reducedBellMatrix.p->e[3].isZeroTerminal()); @@ -670,8 +699,8 @@ TEST(DDPackageTest, Ancillaries) { EXPECT_TRUE(reducedBellMatrix.p->e[0].p->e[3].isZeroTerminal()); dd->incRef(bellMatrix); - reducedBellMatrix = - dd->reduceAncillae(bellMatrix, {false, false, true, true}, false); + reducedBellMatrix = dd->matrices().reduceAncillae( + bellMatrix, {false, false, true, true}, false); EXPECT_TRUE(reducedBellMatrix.p->e[1].isZeroTerminal()); EXPECT_TRUE(reducedBellMatrix.p->e[2].isZeroTerminal()); EXPECT_TRUE(reducedBellMatrix.p->e[3].isZeroTerminal()); @@ -686,7 +715,7 @@ TEST(DDPackageTest, GarbageVector) { auto dd = std::make_unique(4); auto hGate = getDD(qc::StandardOperation(0, qc::H), *dd); auto cxGate = getDD(qc::StandardOperation(0_pc, 1, qc::X), *dd); - auto zeroState = dd->makeZeroState(2); + auto zeroState = dd->vectors().makeZeroState(2); auto bellState = dd->multiply(dd->multiply(cxGate, hGate), zeroState); std::cout << "Bell State:\n"; bellState.printVector(); @@ -724,16 +753,16 @@ TEST(DDPackageTest, GarbageMatrix) { dd->incRef(bellMatrix); auto reducedBellMatrix = - dd->reduceGarbage(bellMatrix, {false, false, false, false}); + dd->matrices().reduceGarbage(bellMatrix, {false, false, false, false}); EXPECT_EQ(bellMatrix, reducedBellMatrix); dd->incRef(bellMatrix); reducedBellMatrix = - dd->reduceGarbage(bellMatrix, {false, false, true, false}); + dd->matrices().reduceGarbage(bellMatrix, {false, false, true, false}); EXPECT_NE(bellMatrix, reducedBellMatrix); dd->incRef(bellMatrix); reducedBellMatrix = - dd->reduceGarbage(bellMatrix, {false, true, false, false}); + dd->matrices().reduceGarbage(bellMatrix, {false, true, false, false}); auto mat = reducedBellMatrix.getMatrix(2); auto zero = CVec{{0., 0.}, {0., 0.}, {0., 0.}, {0., 0.}}; EXPECT_EQ(mat[2], zero); @@ -741,14 +770,14 @@ TEST(DDPackageTest, GarbageMatrix) { dd->incRef(bellMatrix); reducedBellMatrix = - dd->reduceGarbage(bellMatrix, {true, false, false, false}); + dd->matrices().reduceGarbage(bellMatrix, {true, false, false, false}); mat = reducedBellMatrix.getMatrix(2); EXPECT_EQ(mat[1], zero); EXPECT_EQ(mat[3], zero); dd->incRef(bellMatrix); - reducedBellMatrix = - dd->reduceGarbage(bellMatrix, {false, true, false, false}, false); + reducedBellMatrix = dd->matrices().reduceGarbage( + bellMatrix, {false, true, false, false}, false); EXPECT_TRUE(reducedBellMatrix.p->e[1].isZeroTerminal()); EXPECT_TRUE(reducedBellMatrix.p->e[3].isZeroTerminal()); } @@ -757,7 +786,7 @@ TEST(DDPackageTest, ReduceGarbageVector) { auto dd = std::make_unique(3); auto xGate = getDD(qc::StandardOperation(2, qc::X), *dd); auto hGate = getDD(qc::StandardOperation(2, qc::H), *dd); - auto zeroState = dd->makeZeroState(3); + auto zeroState = dd->vectors().makeZeroState(3); auto initialState = dd->multiply(dd->multiply(hGate, xGate), zeroState); std::cout << "Initial State:\n"; initialState.printVector(); @@ -766,13 +795,13 @@ TEST(DDPackageTest, ReduceGarbageVector) { auto reducedState = dd->reduceGarbage(initialState, {false, true, true}); std::cout << "After reduceGarbage():\n"; reducedState.printVector(); - EXPECT_EQ(reducedState, dd->makeZeroState(3)); + EXPECT_EQ(reducedState, dd->vectors().makeZeroState(3)); dd->incRef(initialState); auto reducedState2 = dd->reduceGarbage(initialState, {false, true, true}, true); - EXPECT_EQ(reducedState2, dd->makeZeroState(3)); + EXPECT_EQ(reducedState2, dd->vectors().makeZeroState(3)); } TEST(DDPackageTest, ReduceGarbageVectorTGate) { @@ -782,7 +811,7 @@ TEST(DDPackageTest, ReduceGarbageVectorTGate) { const auto xGate1 = getDD(qc::StandardOperation(1, qc::X), *dd); const auto tdgGate0 = getDD(qc::StandardOperation(0, qc::Tdg), *dd); - auto zeroState = dd->makeZeroState(nqubits); + auto zeroState = dd->vectors().makeZeroState(nqubits); auto initialState = dd->multiply( dd->multiply(tdgGate0, dd->multiply(xGate0, xGate1)), zeroState); std::cout << "Initial State:\n"; @@ -807,8 +836,8 @@ TEST(DDPackageTest, ReduceGarbageMatrix) { initialState.printMatrix(dd->qubits()); dd->incRef(initialState); - auto reducedState1 = - dd->reduceGarbage(initialState, {false, true, true}, true, true); + auto reducedState1 = dd->matrices().reduceGarbage( + initialState, {false, true, true}, true, true); std::cout << "After reduceGarbage(q1 and q2 are garbage):\n"; reducedState1.printMatrix(dd->qubits()); @@ -824,8 +853,8 @@ TEST(DDPackageTest, ReduceGarbageMatrix) { EXPECT_EQ(reducedState1.getMatrix(dd->qubits()), expectedMatrix1); dd->incRef(initialState); - auto reducedState2 = - dd->reduceGarbage(initialState, {true, false, false}, true, true); + auto reducedState2 = dd->matrices().reduceGarbage( + initialState, {true, false, false}, true, true); std::cout << "After reduceGarbage(q0 is garbage):\n"; reducedState2.printMatrix(dd->qubits()); @@ -876,7 +905,7 @@ TEST(DDPackageTest, ReduceGarbageMatrixNoGarbage) { const auto tdgGate0 = getDD(qc::StandardOperation(0, qc::Tdg), *dd); const auto tdgGate1 = getDD(qc::StandardOperation(1, qc::Tdg), *dd); - auto c1 = Package::makeIdent(); + auto c1 = MatrixDDContainer::makeIdent(); auto c2 = dd->multiply(tdgGate0, tdgGate1); std::cout << "c2:\n"; @@ -895,7 +924,7 @@ TEST(DDPackageTest, ReduceGarbageMatrixTGate) { const auto tdgGate0 = getDD(qc::StandardOperation(0, qc::Tdg), *dd); const auto tdgGate1 = getDD(qc::StandardOperation(1, qc::Tdg), *dd); - auto c1 = Package::makeIdent(); + auto c1 = MatrixDDContainer::makeIdent(); auto c2 = dd->multiply(tdgGate0, tdgGate1); std::cout << "c1:\n"; @@ -919,10 +948,12 @@ TEST(DDPackageTest, InvalidMakeBasisStateAndGate) { auto nqubits = 2U; auto dd = std::make_unique(nqubits); auto basisState = std::vector{BasisStates::zero}; - EXPECT_THROW(dd->makeBasisState(nqubits, basisState), std::runtime_error); - EXPECT_THROW(dd->makeZeroState(3), std::runtime_error); - EXPECT_THROW(dd->makeBasisState(3, {true, true, true}), std::runtime_error); - EXPECT_THROW(dd->makeBasisState( + EXPECT_THROW(dd->vectors().makeBasisState(nqubits, basisState), + std::runtime_error); + EXPECT_THROW(dd->vectors().makeZeroState(3), std::runtime_error); + EXPECT_THROW(dd->vectors().makeBasisState(3, {true, true, true}), + std::runtime_error); + EXPECT_THROW(dd->vectors().makeBasisState( 3, {BasisStates::one, BasisStates::one, BasisStates::one}), std::runtime_error); EXPECT_THROW(getDD(qc::StandardOperation(3, qc::X), *dd), std::runtime_error); @@ -942,9 +973,9 @@ TEST(DDPackageTest, PackageReset) { // one node in unique table of variable 0 auto xGate = getDD(qc::StandardOperation(0, qc::X), *dd); - const auto& unique = dd->mUniqueTable.getTables(); + const auto& unique = dd->matrices().getUniqueTable().getTables(); const auto& table = unique[0]; - auto ihash = dd->mUniqueTable.hash(*xGate.p); + auto ihash = dd->matrices().getUniqueTable().hash(*xGate.p); const auto* node = table[ihash]; std::cout << ihash << ": " << reinterpret_cast(xGate.p) << "\n"; // node should be the first in this unique table bucket @@ -970,43 +1001,43 @@ TEST(DDPackageTest, MaxRefCount) { TEST(DDPackageTest, Inverse) { auto dd = std::make_unique(1); auto x = getDD(qc::StandardOperation(0, qc::X), *dd); - auto xdag = dd->conjugateTranspose(x); + auto xdag = dd->matrices().conjugateTranspose(x); EXPECT_EQ(x, xdag); dd->garbageCollect(); // nothing should have been collected since the threshold is not reached - EXPECT_EQ(dd->mUniqueTable.getNumEntries(), 1); + EXPECT_EQ(dd->matrices().getUniqueTable().getNumEntries(), 1); dd->incRef(x); dd->garbageCollect(true); // nothing should have been collected since the lone node has a non-zero ref // count - EXPECT_EQ(dd->mUniqueTable.getNumEntries(), 1); + EXPECT_EQ(dd->matrices().getUniqueTable().getNumEntries(), 1); dd->decRef(x); dd->garbageCollect(true); // now the node should have been collected - EXPECT_EQ(dd->mUniqueTable.getNumEntries(), 0); + EXPECT_EQ(dd->matrices().getUniqueTable().getNumEntries(), 0); } TEST(DDPackageTest, UniqueTableAllocation) { auto dd = std::make_unique(1); - auto allocs = dd->vMemoryManager.getStats().numAllocated; + auto allocs = dd->vectors().getMemoryManager().getStats().numAllocated; std::cout << allocs << "\n"; std::vector nodes{allocs}; // get all the nodes that are pre-allocated for (auto i = 0U; i < allocs; ++i) { - nodes[i] = dd->vMemoryManager.get(); + nodes[i] = dd->vectors().getMemoryManager().get(); } // trigger new allocation - const auto* node = dd->vMemoryManager.get(); + const auto* node = dd->vectors().getMemoryManager().get(); ASSERT_NE(node, nullptr); - EXPECT_EQ(dd->vMemoryManager.getStats().numAllocated, + EXPECT_EQ(dd->vectors().getMemoryManager().getStats().numAllocated, (1. + MemoryManager::GROWTH_FACTOR) * static_cast(allocs)); // clearing the unique table should reduce the allocated size to the original // size - dd->vMemoryManager.reset(); - EXPECT_EQ(dd->vMemoryManager.getStats().numAllocated, allocs); + dd->vectors().getMemoryManager().reset(); + EXPECT_EQ(dd->vectors().getMemoryManager().getStats().numAllocated, allocs); } TEST(DDPackageTest, SpecialCaseTerminal) { @@ -1030,23 +1061,23 @@ TEST(DDPackageTest, SpecialCaseTerminal) { std::filesystem::remove(filename); } - EXPECT_EQ(dd->vUniqueTable.lookup(one.p), one.p); + EXPECT_EQ(dd->vectors().getUniqueTable().lookup(one.p), one.p); auto zero = vEdge::zero(); - EXPECT_TRUE(dd->kronecker(zero, one, 0).isZeroTerminal()); - EXPECT_TRUE(dd->kronecker(one, one, 0).isOneTerminal()); + EXPECT_TRUE(dd->vectors().kronecker(zero, one, 0).isZeroTerminal()); + EXPECT_TRUE(dd->vectors().kronecker(one, one, 0).isOneTerminal()); EXPECT_EQ(one.getValueByPath(0, ""), 1.); EXPECT_EQ(one.getValueByIndex(0), 1.); EXPECT_EQ(mEdge::one().getValueByIndex(0, 0, 0), 1.); - EXPECT_EQ(dd->innerProduct(zero, zero), ComplexValue(0.)); + EXPECT_EQ(dd->vectors().innerProduct(zero, zero), ComplexValue(0.)); } TEST(DDPackageTest, KroneckerProduct) { auto dd = std::make_unique(2); auto x = getDD(qc::StandardOperation(0, qc::X), *dd); - auto kronecker = dd->kronecker(x, x, 1); + auto kronecker = dd->matrices().kronecker(x, x, 1); EXPECT_EQ(kronecker.p->v, 1); EXPECT_TRUE(kronecker.p->e[0].isZeroTerminal()); EXPECT_EQ(kronecker.p->e[0], kronecker.p->e[3]); @@ -1057,16 +1088,16 @@ TEST(DDPackageTest, KroneckerProduct) { EXPECT_TRUE(kronecker.p->e[1].p->e[1].isOneTerminal()); EXPECT_EQ(kronecker.p->e[1].p->e[1], kronecker.p->e[1].p->e[2]); - auto kronecker2 = dd->kronecker(x, x, 1); + auto kronecker2 = dd->matrices().kronecker(x, x, 1); EXPECT_EQ(kronecker, kronecker2); } TEST(DDPackageTest, KroneckerProductVectors) { auto dd = std::make_unique(2); - auto zeroState = dd->makeZeroState(1); - auto kronecker = dd->kronecker(zeroState, zeroState, 1); + auto zeroState = dd->vectors().makeZeroState(1); + auto kronecker = dd->vectors().kronecker(zeroState, zeroState, 1); - auto expected = dd->makeZeroState(2); + auto expected = dd->vectors().makeZeroState(2); EXPECT_EQ(kronecker, expected); } @@ -1075,9 +1106,9 @@ TEST(DDPackageTest, KroneckerIdentityHandling) { // create a Hadamard gate on the middle qubit auto h = getDD(qc::StandardOperation(1U, qc::H), *dd); // create a single qubit identity - auto id = Package::makeIdent(); + auto id = MatrixDDContainer::makeIdent(); // kronecker both DDs - const auto combined = dd->kronecker(h, id, 1); + const auto combined = dd->matrices().kronecker(h, id, 1); const auto matrix = combined.getMatrix(dd->qubits()); const auto expectedMatrix = CMat{ {SQRT2_2, 0, 0, 0, SQRT2_2, 0, 0, 0}, @@ -1096,56 +1127,56 @@ TEST(DDPackageTest, NearZeroNormalize) { auto dd = std::make_unique(2); const fp nearZero = RealNumber::eps / 10; vEdge ve{}; - ve.p = dd->vMemoryManager.get(); + ve.p = dd->vectors().getMemoryManager().get(); ve.p->v = 1; ve.w = Complex::one(); std::array edges{}; for (auto& edge : edges) { - edge.p = dd->vMemoryManager.get(); + edge.p = dd->vectors().getMemoryManager().get(); edge.p->v = 0; edge.w = nearZero; edge.p->e = {vEdge::one(), vEdge::one()}; } - auto veNormalizedCached = - vCachedEdge::normalize(ve.p, edges, dd->vMemoryManager, dd->cn); + auto veNormalizedCached = vCachedEdge::normalize( + ve.p, edges, dd->vectors().getMemoryManager(), dd->cn); EXPECT_EQ(veNormalizedCached, vCachedEdge::zero()); std::array edges2{}; for (auto& edge : edges2) { - edge.p = dd->vMemoryManager.get(); + edge.p = dd->vectors().getMemoryManager().get(); edge.p->v = 0; edge.w = dd->cn.lookup(nearZero); edge.p->e = {vEdge::one(), vEdge::one()}; } auto veNormalized = - vEdge::normalize(ve.p, edges2, dd->vMemoryManager, dd->cn); + vEdge::normalize(ve.p, edges2, dd->vectors().getMemoryManager(), dd->cn); EXPECT_TRUE(veNormalized.isZeroTerminal()); mEdge me{}; - me.p = dd->mMemoryManager.get(); + me.p = dd->matrices().getMemoryManager().get(); me.p->v = 1; me.w = Complex::one(); std::array edges3{}; for (auto& edge : edges3) { - edge.p = dd->mMemoryManager.get(); + edge.p = dd->matrices().getMemoryManager().get(); edge.p->v = 0; edge.w = nearZero; edge.p->e = {mEdge::one(), mEdge::one(), mEdge::one(), mEdge::one()}; } - auto meNormalizedCached = - mCachedEdge::normalize(me.p, edges3, dd->mMemoryManager, dd->cn); + auto meNormalizedCached = mCachedEdge::normalize( + me.p, edges3, dd->matrices().getMemoryManager(), dd->cn); EXPECT_EQ(meNormalizedCached, mCachedEdge::zero()); - me.p = dd->mMemoryManager.get(); + me.p = dd->matrices().getMemoryManager().get(); std::array edges4{}; for (auto& edge : edges4) { - edge.p = dd->mMemoryManager.get(); + edge.p = dd->matrices().getMemoryManager().get(); edge.p->v = 0; edge.w = dd->cn.lookup(nearZero, 0.); edge.p->e = {mEdge::one(), mEdge::one(), mEdge::one(), mEdge::one()}; } auto meNormalized = - mEdge::normalize(me.p, edges4, dd->mMemoryManager, dd->cn); + mEdge::normalize(me.p, edges4, dd->matrices().getMemoryManager(), dd->cn); EXPECT_TRUE(meNormalized.isZeroTerminal()); } @@ -1154,7 +1185,7 @@ TEST(DDPackageTest, DestructiveMeasurementAll) { auto hGate0 = getDD(qc::StandardOperation(0, qc::H), *dd); auto hGate1 = getDD(qc::StandardOperation(1, qc::H), *dd); auto plusMatrix = dd->multiply(hGate0, hGate1); - auto zeroState = dd->makeZeroState(2); + auto zeroState = dd->vectors().makeZeroState(2); auto plusState = dd->multiply(plusMatrix, zeroState); dd->incRef(plusState); @@ -1179,7 +1210,7 @@ TEST(DDPackageTest, DestructiveMeasurementOne) { auto hGate0 = getDD(qc::StandardOperation(0, qc::H), *dd); auto hGate1 = getDD(qc::StandardOperation(1, qc::H), *dd); auto plusMatrix = dd->multiply(hGate0, hGate1); - auto zeroState = dd->makeZeroState(2); + auto zeroState = dd->vectors().makeZeroState(2); auto plusState = dd->multiply(plusMatrix, zeroState); dd->incRef(plusState); @@ -1323,7 +1354,7 @@ TEST(DDPackageTest, BasicNumericStabilityTest) { auto dd = std::make_unique(1); auto tol = RealNumber::eps; ComplexNumbers::setTolerance(limits::epsilon()); - auto state = dd->makeZeroState(1); + auto state = dd->vectors().makeZeroState(1); auto h = getDD(qc::StandardOperation(0, qc::H), *dd); auto state1 = dd->multiply(h, state); auto z = getDD(qc::StandardOperation(0, qc::Z), *dd); @@ -1357,8 +1388,8 @@ TEST(DDPackageTest, NormalizationNumericStabilityTest) { auto pdag = getDD(qc::StandardOperation(0, qc::P, {-lambda}), *dd); auto result = dd->multiply(p, pdag); EXPECT_TRUE(result.isIdentity()); - dd->cUniqueTable.clear(); - dd->cMemoryManager.reset(); + dd->cUt.clear(); + dd->cMm.reset(); } } @@ -1368,7 +1399,7 @@ TEST(DDPackageTest, FidelityOfMeasurementOutcomes) { const auto hGate = getDD(qc::StandardOperation(2, qc::H), *dd); const auto cxGate1 = getDD(qc::StandardOperation(2_pc, 1, qc::X), *dd); const auto cxGate2 = getDD(qc::StandardOperation(1_pc, 0, qc::X), *dd); - const auto zeroState = dd->makeZeroState(3); + const auto zeroState = dd->vectors().makeZeroState(3); const auto ghzState = dd->multiply( cxGate2, dd->multiply(cxGate1, dd->multiply(hGate, zeroState))); @@ -1376,41 +1407,42 @@ TEST(DDPackageTest, FidelityOfMeasurementOutcomes) { SparsePVec probs{}; probs[0] = 0.5; probs[7] = 0.5; - const auto fidelity = Package::fidelityOfMeasurementOutcomes(ghzState, probs); + const auto fidelity = + VectorDDContainer::fidelityOfMeasurementOutcomes(ghzState, probs); EXPECT_NEAR(fidelity, 1.0, RealNumber::eps); } TEST(DDPackageTest, CloseToIdentity) { auto dd = std::make_unique(3); - auto id = Package::makeIdent(); - EXPECT_TRUE(dd->isCloseToIdentity(id)); + auto id = MatrixDDContainer::makeIdent(); + EXPECT_TRUE(dd->matrices().isCloseToIdentity(id)); mEdge close{}; close.p = id.p; close.w = dd->cn.lookup(1e-11, 0); auto id2 = dd->makeDDNode(1, std::array{id, mEdge::zero(), mEdge::zero(), close}); - EXPECT_TRUE(dd->isCloseToIdentity(id2)); + EXPECT_TRUE(dd->matrices().isCloseToIdentity(id2)); auto noId = dd->makeDDNode(1, std::array{mEdge::zero(), id, mEdge::zero(), close}); - EXPECT_FALSE(dd->isCloseToIdentity(noId)); + EXPECT_FALSE(dd->matrices().isCloseToIdentity(noId)); mEdge notClose{}; notClose.p = id.p; notClose.w = dd->cn.lookup(1e-9, 0); auto noId2 = dd->makeDDNode( 1, std::array{notClose, mEdge::zero(), mEdge::zero(), close}); - EXPECT_FALSE(dd->isCloseToIdentity(noId2)); + EXPECT_FALSE(dd->matrices().isCloseToIdentity(noId2)); auto noId3 = dd->makeDDNode( 1, std::array{close, mEdge::zero(), mEdge::zero(), notClose}); - EXPECT_FALSE(dd->isCloseToIdentity(noId3)); + EXPECT_FALSE(dd->matrices().isCloseToIdentity(noId3)); auto notClose2 = dd->makeDDNode( 0, std::array{mEdge::zero(), mEdge::one(), mEdge::one(), mEdge::zero()}); auto notClose3 = dd->makeDDNode( 1, std::array{notClose2, mEdge::zero(), mEdge::zero(), notClose2}); - EXPECT_FALSE(dd->isCloseToIdentity(notClose3)); + EXPECT_FALSE(dd->matrices().isCloseToIdentity(notClose3)); } TEST(DDPackageTest, CloseToIdentityWithGarbageAtTheBeginning) { @@ -1432,10 +1464,10 @@ TEST(DDPackageTest, CloseToIdentityWithGarbageAtTheBeginning) { auto c1MultipliedWithC2 = dd->multiply(c1, dd->conjugateTranspose(c2)); - EXPECT_TRUE(dd->isCloseToIdentity(c1MultipliedWithC2, tol, - {false, true, true}, false)); - EXPECT_FALSE(dd->isCloseToIdentity(c1MultipliedWithC2, tol, - {false, false, true}, false)); + EXPECT_TRUE(dd->matrices().isCloseToIdentity(c1MultipliedWithC2, tol, + {false, true, true}, false)); + EXPECT_FALSE(dd->matrices().isCloseToIdentity(c1MultipliedWithC2, tol, + {false, false, true}, false)); } TEST(DDPackageTest, CloseToIdentityWithGarbageAtTheEnd) { @@ -1459,12 +1491,12 @@ TEST(DDPackageTest, CloseToIdentityWithGarbageAtTheEnd) { const auto c3MultipliedWithC4 = dd->multiply(c3, dd->conjugateTranspose(c4)); - EXPECT_FALSE(dd->isCloseToIdentity(c3MultipliedWithC4, tol, - {false, true, true}, false)); - EXPECT_FALSE(dd->isCloseToIdentity(c3MultipliedWithC4, tol, - {true, false, true}, false)); - EXPECT_TRUE(dd->isCloseToIdentity(c3MultipliedWithC4, tol, - {true, true, false}, false)); + EXPECT_FALSE(dd->matrices().isCloseToIdentity(c3MultipliedWithC4, tol, + {false, true, true}, false)); + EXPECT_FALSE(dd->matrices().isCloseToIdentity(c3MultipliedWithC4, tol, + {true, false, true}, false)); + EXPECT_TRUE(dd->matrices().isCloseToIdentity(c3MultipliedWithC4, tol, + {true, true, false}, false)); } TEST(DDPackageTest, CloseToIdentityWithGarbageInTheMiddle) { @@ -1488,12 +1520,12 @@ TEST(DDPackageTest, CloseToIdentityWithGarbageInTheMiddle) { const auto c5MultipliedWithC6 = dd->multiply(c5, dd->conjugateTranspose(c6)); - EXPECT_FALSE(dd->isCloseToIdentity(c5MultipliedWithC6, tol, - {false, true, true}, false)); - EXPECT_FALSE(dd->isCloseToIdentity(c5MultipliedWithC6, tol, - {true, true, false}, false)); - EXPECT_TRUE(dd->isCloseToIdentity(c5MultipliedWithC6, tol, - {true, false, true}, false)); + EXPECT_FALSE(dd->matrices().isCloseToIdentity(c5MultipliedWithC6, tol, + {false, true, true}, false)); + EXPECT_FALSE(dd->matrices().isCloseToIdentity(c5MultipliedWithC6, tol, + {true, true, false}, false)); + EXPECT_TRUE(dd->matrices().isCloseToIdentity(c5MultipliedWithC6, tol, + {true, false, true}, false)); } TEST(DDPackageTest, dNodeMultiply) { @@ -1502,7 +1534,7 @@ TEST(DDPackageTest, dNodeMultiply) { const auto dd = std::make_unique( nrQubits, DENSITY_MATRIX_SIMULATOR_DD_PACKAGE_CONFIG); // Make zero density matrix - auto state = dd->makeZeroDensityOperator(dd->qubits()); + auto state = dd->densities().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)); @@ -1549,7 +1581,7 @@ TEST(DDPackageTest, dNodeMultiply2) { const auto dd = std::make_unique( nrQubits, DENSITY_MATRIX_SIMULATOR_DD_PACKAGE_CONFIG); // Make zero density matrix - auto state = dd->makeZeroDensityOperator(dd->qubits()); + auto state = dd->densities().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)); @@ -1589,12 +1621,12 @@ TEST(DDPackageTest, dNodeMulCache1) { const auto dd = std::make_unique( nrQubits, DENSITY_MATRIX_SIMULATOR_DD_PACKAGE_CONFIG); // Make zero density matrix - auto state = dd->makeZeroDensityOperator(nrQubits); + auto state = dd->densities().makeZeroDensityOperator(nrQubits); const auto operation = getDD(qc::StandardOperation(0, qc::H), *dd); dd->applyOperationToDensity(state, operation); - state = dd->makeZeroDensityOperator(nrQubits); + state = dd->densities().makeZeroDensityOperator(nrQubits); auto& computeTable = dd->getMultiplicationComputeTable(); const auto& densityMatrix0 = @@ -1635,7 +1667,7 @@ TEST(DDPackageTest, dNoiseCache) { constexpr auto nrQubits = 1U; const auto dd = std::make_unique(nrQubits); // Make zero density matrix - const auto initialState = dd->makeZeroDensityOperator(nrQubits); + const auto initialState = dd->densities().makeZeroDensityOperator(nrQubits); // nothing pre-cached const std::vector target = {0}; @@ -1710,7 +1742,7 @@ TEST(DDPackageTest, dStochCache) { TEST(DDPackageTest, stateFromVectorBell) { const auto dd = std::make_unique(2); const auto v = std::vector>{SQRT2_2, 0, 0, SQRT2_2}; - const auto s = dd->makeStateFromVector(v); + const auto s = dd->vectors().makeStateFromVector(v); ASSERT_NE(s.p, nullptr); EXPECT_EQ(s.p->v, 1); EXPECT_EQ(s.p->e[0].w.r->value, SQRT2_2); @@ -1732,18 +1764,18 @@ TEST(DDPackageTest, stateFromVectorBell) { TEST(DDPackageTest, stateFromVectorEmpty) { auto dd = std::make_unique(1); auto v = std::vector>{}; - EXPECT_TRUE(dd->makeStateFromVector(v).isOneTerminal()); + EXPECT_TRUE(dd->vectors().makeStateFromVector(v).isOneTerminal()); } TEST(DDPackageTest, stateFromVectorNoPowerOfTwo) { auto dd = std::make_unique(3); auto v = std::vector>{1, 2, 3, 4, 5}; - EXPECT_THROW(dd->makeStateFromVector(v), std::invalid_argument); + EXPECT_THROW(dd->vectors().makeStateFromVector(v), std::invalid_argument); } TEST(DDPackageTest, stateFromScalar) { const auto dd = std::make_unique(1); - const auto s = dd->makeStateFromVector({1}); + const auto s = dd->vectors().makeStateFromVector({1}); EXPECT_TRUE(s.isTerminal()); EXPECT_EQ(s.w.r->value, 1); EXPECT_EQ(s.w.i->value, 0); @@ -1753,7 +1785,7 @@ TEST(DDPackageTest, expectationValueGlobalOperators) { constexpr Qubit maxQubits = 3; const auto dd = std::make_unique(maxQubits); for (Qubit nrQubits = 1; nrQubits < maxQubits + 1; ++nrQubits) { - const auto zeroState = dd->makeZeroState(nrQubits); + const auto zeroState = dd->vectors().makeZeroState(nrQubits); // Definition global operators const auto singleSiteX = getDD(qc::StandardOperation(0, qc::X), *dd); @@ -1783,7 +1815,7 @@ TEST(DDPackageTest, expectationValueLocalOperators) { constexpr Qubit maxQubits = 3; const auto dd = std::make_unique(maxQubits); for (Qubit nrQubits = 1; nrQubits < maxQubits + 1; ++nrQubits) { - const auto zeroState = dd->makeZeroState(nrQubits); + const auto zeroState = dd->vectors().makeZeroState(nrQubits); // Local expectation values at each site for (Qubit site = 0; site < nrQubits - 1; ++site) { @@ -1803,7 +1835,7 @@ TEST(DDPackageTest, expectationValueExceptions) { constexpr auto nrQubits = 2U; const auto dd = std::make_unique(nrQubits); - const auto zeroState = dd->makeZeroState(nrQubits - 1); + const auto zeroState = dd->vectors().makeZeroState(nrQubits - 1); const auto xGate = getDD(qc::StandardOperation(1, qc::X), *dd); EXPECT_ANY_THROW(dd->expectationValue(xGate, zeroState)); @@ -1814,7 +1846,7 @@ TEST(DDPackageTest, DDFromSingleQubitMatrix) { constexpr auto nrQubits = 1U; const auto dd = std::make_unique(nrQubits); - const auto matDD = dd->makeDDFromMatrix(inputMatrix); + const auto matDD = dd->matrices().makeDDFromMatrix(inputMatrix); const auto outputMatrix = matDD.getMatrix(dd->qubits()); EXPECT_EQ(inputMatrix, outputMatrix); @@ -1826,7 +1858,7 @@ TEST(DDPackageTest, DDFromTwoQubitMatrix) { constexpr auto nrQubits = 2U; const auto dd = std::make_unique(nrQubits); - const auto matDD = dd->makeDDFromMatrix(inputMatrix); + const auto matDD = dd->matrices().makeDDFromMatrix(inputMatrix); const auto outputMatrix = matDD.getMatrix(dd->qubits()); EXPECT_EQ(inputMatrix, outputMatrix); @@ -1840,7 +1872,7 @@ TEST(DDPackageTest, DDFromTwoQubitAsymmetricalMatrix) { constexpr auto nrQubits = 2U; const auto dd = std::make_unique(nrQubits); - const auto matDD = dd->makeDDFromMatrix(inputMatrix); + const auto matDD = dd->matrices().makeDDFromMatrix(inputMatrix); const auto outputMatrix = matDD.getMatrix(dd->qubits()); EXPECT_EQ(inputMatrix, outputMatrix); @@ -1855,7 +1887,7 @@ TEST(DDPackageTest, DDFromThreeQubitMatrix) { constexpr auto nrQubits = 3U; const auto dd = std::make_unique(nrQubits); - const auto matDD = dd->makeDDFromMatrix(inputMatrix); + const auto matDD = dd->matrices().makeDDFromMatrix(inputMatrix); const auto outputMatrix = matDD.getMatrix(dd->qubits()); @@ -1867,7 +1899,7 @@ TEST(DDPackageTest, DDFromEmptyMatrix) { constexpr auto nrQubits = 3U; const auto dd = std::make_unique(nrQubits); - EXPECT_TRUE(dd->makeDDFromMatrix(inputMatrix).isOneTerminal()); + EXPECT_TRUE(dd->matrices().makeDDFromMatrix(inputMatrix).isOneTerminal()); } TEST(DDPackageTest, DDFromNonPowerOfTwoMatrix) { @@ -1875,7 +1907,8 @@ TEST(DDPackageTest, DDFromNonPowerOfTwoMatrix) { constexpr auto nrQubits = 3U; const auto dd = std::make_unique(nrQubits); - EXPECT_THROW(dd->makeDDFromMatrix(inputMatrix), std::invalid_argument); + EXPECT_THROW(dd->matrices().makeDDFromMatrix(inputMatrix), + std::invalid_argument); } TEST(DDPackageTest, DDFromNonSquareMatrix) { @@ -1883,7 +1916,8 @@ TEST(DDPackageTest, DDFromNonSquareMatrix) { constexpr auto nrQubits = 3U; const auto dd = std::make_unique(nrQubits); - EXPECT_THROW(dd->makeDDFromMatrix(inputMatrix), std::invalid_argument); + EXPECT_THROW(dd->matrices().makeDDFromMatrix(inputMatrix), + std::invalid_argument); } TEST(DDPackageTest, DDFromSingleElementMatrix) { @@ -1892,7 +1926,7 @@ TEST(DDPackageTest, DDFromSingleElementMatrix) { constexpr auto nrQubits = 1U; const auto dd = std::make_unique(nrQubits); - EXPECT_TRUE(dd->makeDDFromMatrix(inputMatrix).isOneTerminal()); + EXPECT_TRUE(dd->matrices().makeDDFromMatrix(inputMatrix).isOneTerminal()); } constexpr TwoQubitGateMatrix CX_MAT{ @@ -1917,9 +1951,9 @@ TEST(DDPackageTest, TwoQubitControlledGateDDConstruction) { if (control == target) { continue; } - const auto controlledGateDD = - dd->makeTwoQubitGateDD(controlledGateMatrix, control, target); - const auto gateDD = dd->makeGateDD( + const auto controlledGateDD = dd->matrices().makeTwoQubitGateDD( + controlledGateMatrix, control, target); + const auto gateDD = dd->matrices().makeGateDD( gateMatrix, qc::Control{static_cast(control)}, target); EXPECT_EQ(controlledGateDD, gateDD); } @@ -2120,7 +2154,7 @@ TEST(DDPackageTest, RZZGateDDConstruction) { } } - const auto identity = Package::makeIdent(); + const auto identity = MatrixDDContainer::makeIdent(); const auto rzzZero = getDD(qc::StandardOperation({0, 1}, qc::RZZ, {0.}), *dd); EXPECT_EQ(rzzZero, identity); @@ -2172,7 +2206,7 @@ TEST(DDPackageTest, RYYGateDDConstruction) { } } - const auto identity = Package::makeIdent(); + const auto identity = MatrixDDContainer::makeIdent(); const auto ryyZero = getDD(qc::StandardOperation({0, 1}, qc::RYY, {0.}), *dd); EXPECT_EQ(ryyZero, identity); @@ -2215,7 +2249,7 @@ TEST(DDPackageTest, RXXGateDDConstruction) { } } - const auto identity = Package::makeIdent(); + const auto identity = MatrixDDContainer::makeIdent(); const auto rxxZero = getDD(qc::StandardOperation({0, 1}, qc::RXX, {0.}), *dd); EXPECT_EQ(rxxZero, identity); @@ -2256,7 +2290,7 @@ TEST(DDPackageTest, RZXGateDDConstruction) { } } - const auto identity = Package::makeIdent(); + const auto identity = MatrixDDContainer::makeIdent(); const auto rzxZero = getDD(qc::StandardOperation({0, 1}, qc::RZX, {0.}), *dd); EXPECT_EQ(rzxZero, identity); @@ -2430,7 +2464,7 @@ TEST(DDPackageTest, InnerProductTopNodeConjugation) { // Ising model evolution up to a time T=1 constexpr auto nrQubits = 2U; const auto dd = std::make_unique(nrQubits); - const auto zeroState = dd->makeZeroState(nrQubits); + const auto zeroState = dd->vectors().makeZeroState(nrQubits); const auto rxx = getDD(qc::StandardOperation({0, 1}, qc::RXX, {-2}), *dd); const auto op = getDD(qc::StandardOperation(0, qc::Z), *dd); @@ -2453,11 +2487,11 @@ TEST(DDPackageTest, DDNodeLeakRegressionTest) { constexpr auto nqubits = 1U; auto dd = std::make_unique(nqubits); - const auto dd1 = dd->makeGateDD(MEAS_ZERO_MAT, 0U); - const auto dd2 = dd->makeGateDD(MEAS_ONE_MAT, 0U); + const auto dd1 = dd->matrices().makeGateDD(MEAS_ZERO_MAT, 0U); + const auto dd2 = dd->matrices().makeGateDD(MEAS_ONE_MAT, 0U); dd->multiply(dd1, dd2); dd->garbageCollect(true); - EXPECT_EQ(dd->mMemoryManager.getStats().numUsed, 0U); + EXPECT_EQ(dd->matrices().getMemoryManager().getStats().numUsed, 0U); } /** @@ -2482,7 +2516,7 @@ TEST(DDPackageTest, CTPerformanceRegressionTest) { // This additional check makes sure that no nodes are leaked. dd->garbageCollect(true); - EXPECT_EQ(dd->mMemoryManager.getStats().numUsed, 0U); + EXPECT_EQ(dd->matrices().getMemoryManager().getStats().numUsed, 0U); } TEST(DDPackageTest, DataStructureStatistics) { @@ -2535,9 +2569,9 @@ TEST(DDPackageTest, ReduceAncillaRegression) { const auto dd = std::make_unique(2); const auto inputMatrix = CMat{{1, 1, 1, 1}, {1, -1, 1, -1}, {1, 1, -1, -1}, {1, -1, -1, 1}}; - auto inputDD = dd->makeDDFromMatrix(inputMatrix); + auto inputDD = dd->matrices().makeDDFromMatrix(inputMatrix); dd->incRef(inputDD); - const auto outputDD = dd->reduceAncillae(inputDD, {true, false}); + const auto outputDD = dd->matrices().reduceAncillae(inputDD, {true, false}); const auto outputMatrix = outputDD.getMatrix(dd->qubits()); const auto expected = @@ -2560,10 +2594,10 @@ TEST(DDPackageTest, VectorConjugate) { {0., -0.5}, {-0.5 * SQRT2_2, -0.5 * SQRT2_2}}; - const auto vecDD = dd->makeStateFromVector(vec); + const auto vecDD = dd->vectors().makeStateFromVector(vec); std::cout << "Vector:\n"; vecDD.printVector(); - const auto conjVecDD = dd->conjugate(vecDD); + const auto conjVecDD = dd->vectors().conjugate(vecDD); std::cout << "Conjugated vector:\n"; conjVecDD.printVector(); @@ -2577,8 +2611,8 @@ TEST(DDPackageTest, VectorConjugate) { TEST(DDPackageTest, ReduceAncillaIdentity) { const auto dd = std::make_unique(2); - auto inputDD = Package::makeIdent(); - const auto outputDD = dd->reduceAncillae(inputDD, {true, true}); + auto inputDD = MatrixDDContainer::makeIdent(); + const auto outputDD = dd->matrices().reduceAncillae(inputDD, {true, true}); const auto outputMatrix = outputDD.getMatrix(dd->qubits()); const auto expected = @@ -2590,7 +2624,7 @@ TEST(DDPackageTest, ReduceAncillaIdentity) { TEST(DDPackageTest, ReduceAnicllaIdentityBeforeFirstNode) { const auto dd = std::make_unique(2); auto xGate = getDD(qc::StandardOperation(0, qc::X), *dd); - const auto outputDD = dd->reduceAncillae(xGate, {false, true}); + const auto outputDD = dd->matrices().reduceAncillae(xGate, {false, true}); const auto outputMatrix = outputDD.getMatrix(dd->qubits()); const auto expected = @@ -2602,7 +2636,7 @@ TEST(DDPackageTest, ReduceAnicllaIdentityAfterLastNode) { const auto dd = std::make_unique(2); auto xGate = getDD(qc::StandardOperation(1, qc::X), *dd); dd->incRef(xGate); - const auto outputDD = dd->reduceAncillae(xGate, {true, false}); + const auto outputDD = dd->matrices().reduceAncillae(xGate, {true, false}); const auto outputMatrix = outputDD.getMatrix(dd->qubits()); const auto expected = @@ -2617,7 +2651,8 @@ TEST(DDPackageTest, ReduceAncillaIdentityBetweenTwoNodes) { auto state = dd->multiply(xGate0, xGate2); dd->incRef(state); - const auto outputDD = dd->reduceAncillae(state, {false, true, false}); + const auto outputDD = + dd->matrices().reduceAncillae(state, {false, true, false}); const auto outputMatrix = outputDD.getMatrix(dd->qubits()); const auto expected = CMat{{0, 0, 0, 0, 0, 1, 0, 0}, {0, 0, 0, 0, 1, 0, 0, 0}, @@ -2629,7 +2664,7 @@ TEST(DDPackageTest, ReduceAncillaIdentityBetweenTwoNodes) { TEST(DDPackageTest, ReduceGarbageIdentity) { const auto dd = std::make_unique(2); - auto inputDD = Package::makeIdent(); + auto inputDD = MatrixDDContainer::makeIdent(); auto outputDD = dd->reduceGarbage(inputDD, {true, true}); auto outputMatrix = outputDD.getMatrix(dd->qubits());