diff --git a/CMakeLists.txt b/CMakeLists.txt index a76b76a..944dbf1 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -52,5 +52,5 @@ endif() configure_file(${CMAKE_CURRENT_SOURCE_DIR}/cmake/cmake_uninstall.cmake.in ${CMAKE_CURRENT_BINARY_DIR}/cmake_uninstall.cmake IMMEDIATE @ONLY) -add_custom_target(uninstall COMMAND ${CMAKE_COMMAND} -P - ${CMAKE_CURRENT_BINARY_DIR}/cmake_uninstall.cmake) +add_custom_target(uninstall-debug COMMAND ${CMAKE_COMMAND} -P + ${CMAKE_CURRENT_BINARY_DIR}/cmake_uninstall.cmake) diff --git a/cmake/ExternalDependencies.cmake b/cmake/ExternalDependencies.cmake index f2684b3..4d676f8 100644 --- a/cmake/ExternalDependencies.cmake +++ b/cmake/ExternalDependencies.cmake @@ -17,6 +17,7 @@ if(BUILD_MQT_DEBUG_BINDINGS) find_package(pybind11 2.13 CONFIG REQUIRED) endif() +# ---------------------------------------------------------------------------------Fetch MQT Core # cmake-format: off set(MQT_CORE_VERSION 2.5.2 CACHE STRING "MQT Core version") @@ -26,6 +27,7 @@ set(MQT_CORE_REPO_OWNER "cda-tum" CACHE STRING "MQT Core repository owner (change when using a fork)") # cmake-format: on if(CMAKE_VERSION VERSION_GREATER_EQUAL 3.24) + # Fetch MQT Core FetchContent_Declare( mqt-core GIT_REPOSITORY https://github.com/${MQT_CORE_REPO_OWNER}/mqt-core.git @@ -42,6 +44,43 @@ else() endif() endif() +# ---------------------------------------------------------------------------------Fetch Eigen3 +# cmake-format: off +set(EIGEN_VERSION 3.4.0 + CACHE STRING "Eigen3 version") +# cmake-format: on +if(CMAKE_VERSION VERSION_GREATER_EQUAL 3.24) + # Fetch Eigen3 + FetchContent_Declare( + Eigen3 + GIT_REPOSITORY https://gitlab.com/libeigen/eigen.git + GIT_TAG ${EIGEN_VERSION} + GIT_SHALLOW TRUE) + list(APPEND FETCH_PACKAGES Eigen3) + set(EIGEN_BUILD_TESTING + OFF + CACHE BOOL "Disable testing for Eigen") + set(EIGEN_BUILD_DOC + OFF + CACHE BOOL "Disable documentation build for Eigen") +else() + find_package(Eigen3 ${EIGEN3_VERSION} REQUIRED NO_MODULE) + if(NOT Eigen3_FOUND) + FetchContent_Declare( + Eigen3 + GIT_REPOSITORY https://gitlab.com/libeigen/eigen.git + GIT_TAG ${EIGEN3_VERSION} + GIT_SHALLOW TRUE) + list(APPEND FETCH_PACKAGES Eigen3) + set(EIGEN_BUILD_TESTING + OFF + CACHE BOOL "Disable testing for Eigen") + set(EIGEN_BUILD_DOC + OFF + CACHE BOOL "Disable documentation build for Eigen") + endif() +endif() + if(BUILD_MQT_DEBUG_TESTS) set(gtest_force_shared_crt ON diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index e460df1..abb9c7b 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -15,6 +15,7 @@ target_include_directories(${PROJECT_NAME} PUBLIC ${PROJECT_SOURCE_DIR}/include # link to the MQT::Core libraries target_link_libraries(${PROJECT_NAME} PUBLIC MQT::CoreDD MQT::CoreIR MQT::CoreCircuitOptimizer) target_link_libraries(${PROJECT_NAME} PRIVATE MQT::ProjectWarnings MQT::ProjectOptions) +target_link_libraries(${PROJECT_NAME} PRIVATE Eigen3::Eigen) # add MQT alias add_library(MQT::Debug ALIAS ${PROJECT_NAME}) diff --git a/src/backend/dd/DDSimDebug.cpp b/src/backend/dd/DDSimDebug.cpp index 283449e..cab289e 100644 --- a/src/backend/dd/DDSimDebug.cpp +++ b/src/backend/dd/DDSimDebug.cpp @@ -16,6 +16,7 @@ #include "ir/operations/ClassicControlledOperation.hpp" #include "ir/operations/OpType.hpp" +#include #include #include #include @@ -848,6 +849,20 @@ Complex complexDivision(const Complex& c1, const Complex& c2) { return {real, imaginary}; } +Complex complexMultiplication(const Complex& c1, const Complex& c2) { + const double real = c1.real * c2.real - c1.imaginary * c2.imaginary; + const double imaginary = c1.real * c2.imaginary + c1.imaginary * c2.real; + return {real, imaginary}; +} + +Complex complexConjugate(const Complex& c) { return {c.real, -c.imaginary}; } + +Complex complexAddition(const Complex& c1, const Complex& c2) { + const double real = c1.real + c2.real; + const double imaginary = c1.imaginary + c2.imaginary; + return {real, imaginary}; +} + bool areComplexEqual(const Complex& c1, const Complex& c2) { const double epsilon = 0.00000001; return std::abs(c1.real - c2.real) < epsilon && @@ -871,7 +886,167 @@ double dotProduct(const Statevector& sv1, const Statevector& sv2) { return complexMagnitude(result); } -bool areQubitsEntangled(Statevector* sv) { +size_t bitStringToNumber(const std::vector& bits) { + size_t result = 0; + for (size_t i = 0; i < bits.size(); i++) { + result |= (bits[i] ? 1ULL : 0ULL) << i; + } + return result; +} + +std::vector extractBits(const std::vector& indices, + size_t value) { + std::vector result; + for (const auto& idx : indices) { + result.push_back(((value >> idx) & 1) == 1); + } + return result; +} + +std::pair splitBitString(size_t number, size_t n, + const std::vector& bits) { + size_t lenFirst = 0; + size_t lenSecond = 0; + + size_t first = 0; + size_t second = 0; + + for (size_t index = 0; index < n; index++) { + if (std::find(bits.begin(), bits.end(), index) != bits.end()) { + first |= (number & 1) << lenFirst; + lenFirst++; + } else { + second |= (number & 1) << lenSecond; + lenSecond++; + } + number >>= 1; + } + return {first, second}; +} + +std::vector> +getPartialTrace(const std::vector>& matrix, + const std::vector& indicesToTraceOut, size_t nQubits) { + const auto traceSize = 1ULL << (nQubits - indicesToTraceOut.size()); + std::vector> traceMatrix( + traceSize, std::vector(traceSize, {0, 0})); + for (size_t i = 0; i < matrix.size(); i++) { + for (size_t j = 0; j < matrix.size(); j++) { + const auto split1 = splitBitString(i, nQubits, indicesToTraceOut); + const auto split2 = splitBitString(j, nQubits, indicesToTraceOut); + if (split1.first != split2.first) { + continue; + } + traceMatrix[split1.second][split2.second] = complexAddition( + traceMatrix[split1.second][split2.second], matrix[i][j]); + } + } + return traceMatrix; +} + +double getEntropy(const std::vector>& matrix) { + Eigen::Matrix, Eigen::Dynamic, Eigen::Dynamic> mat( + static_cast(matrix.size()), static_cast(matrix.size())); + for (size_t i = 0; i < matrix.size(); i++) { + for (size_t j = 0; j < matrix.size(); j++) { + mat(static_cast(i), static_cast(j)) = { + matrix[i][j].real, matrix[i][j].imaginary}; + } + } + + Eigen::ComplexEigenSolver< + Eigen::Matrix, Eigen::Dynamic, Eigen::Dynamic>> + solver(mat); + const Eigen::Vector, Eigen::Dynamic> eigenvalues = + solver.eigenvalues(); + double entropy = 0; + for (const auto val : eigenvalues) { + const auto value = + (val.real() > -0.00001 && val.real() < 0) ? 0 : val.real(); + if (value < 0) { + throw std::runtime_error("Negative eigenvalue"); + } + if (value != 0) { + entropy -= val.real() * log2(val.real()); + } + } + return entropy; +} + +double getSharedInformation(const std::vector>& matrix) { + const auto p0 = getPartialTrace(matrix, {1}, 2); + const auto p1 = getPartialTrace(matrix, {0}, 2); + return getEntropy(p0) + getEntropy(p1) - getEntropy(matrix); +} + +bool areQubitsEntangled(std::vector>& densityMatrix, + size_t qubit1, size_t qubit2) { + const auto numQubits = static_cast(std::log2(densityMatrix.size())); + if (numQubits == 2) { + return getSharedInformation(densityMatrix) > 0; + } + std::vector toTraceOut; + for (size_t i = 0; i < numQubits; i++) { + if (i != qubit1 && i != qubit2) { + toTraceOut.push_back(i); + } + } + const auto partialTrace = + getPartialTrace(densityMatrix, toTraceOut, numQubits); + return getSharedInformation(partialTrace) > 0; +} + +/*std::vector> getPartialTrace(const Statevector& sv, const +std::vector& traceOut) { const auto traceSize = 1ULL << (sv.numQubits - +traceOut.size()); std::vector> traceMatrix(traceSize, +std::vector(traceSize, {0, 0})); for(size_t i = 0; i < sv.numStates; +i++) { for(size_t j = 0; j < sv.numStates; j++) { const auto bitString1 = +extractBits(traceOut, i); const auto bitString2 = extractBits(traceOut, j); + const auto traced1 = bitStringToNumber(bitString1); + const auto traced2 = bitStringToNumber(bitString2); + if (traced1 != traced2) { + continue; + } + const auto product = complexMultiplication(sv.amplitudes[i], +sv.amplitudes[j]); const auto row = i - traced1; const auto col = j - traced2; + traceMatrix[row][col] = complexAddition(traceMatrix[row][col], product); + } + } + return traceMatrix; +} + +Complex getTraceOfSquare(const std::vector>& matrix) { + Complex runningSum{0, 0}; + for(size_t i = 0; i < matrix.size(); i++) { + for(size_t k = 0; k < matrix.size(); k++) { + runningSum = complexAddition(runningSum, +complexMultiplication(matrix[i][k], matrix[k][i])); + } + } + const double epsilon = 0.00000001; + runningSum; +} + +bool partialTraceIsPure(const Statevector& sv, const std::vector& +traceOut) { const auto traceMatrix = getPartialTrace(sv, traceOut); const auto +trace = getTraceOfSquare(traceMatrix); const double epsilon = 0.00000001; return +runningSum.imaginary < epsilon && runningSum.imaginary > -epsilon && +(runningSum.real - 1) < epsilon && (runningSum.real - 1) > -epsilon; +} + +bool areQubitsEntangled(Statevector* sv, size_t qubit1, size_t qubit2) { + std::vector toTraceOut; + for(size_t i = 0; i < sv->numQubits; i++) { + if(i != qubit1 && i != qubit2) { + toTraceOut.push_back(i); + } + } + const auto traceMatrix = getPartialTrace(*sv, toTraceOut); + const auto trace = getTraceOfSquare(traceMatrix); + assert(trace.imaginary == 0); + const Complex factor = trace.real >= 0 ? Complex{std::sqrt(trace.real), 0} : +Complex{0, std::sqrt(-trace.real)}; + const double epsilon = 0.0001; const Span amplitudes(sv->amplitudes, sv->numStates); const bool canBe00 = complexMagnitude(amplitudes[0]) > epsilon; @@ -892,67 +1067,82 @@ bool areQubitsEntangled(Statevector* sv) { (canBe01 && canBe10 && !(canBe00 && canBe11)); } +bool qubitIsSeparable(const Statevector& sv, size_t qubit) { + std::vector toTraceOut; + for(size_t i = 0; i < sv.numQubits; i++) { + if (i != qubit) { + toTraceOut.push_back(i); + } + } + return partialTraceIsPure(sv, toTraceOut); +}*/ + bool checkAssertionEntangled( DDSimulationState* ddsim, std::unique_ptr& assertion) { + Statevector sv; + sv.numQubits = ddsim->interface.getNumQubits(&ddsim->interface); + sv.numStates = 1 << sv.numQubits; + std::vector amplitudes(sv.numStates); + sv.amplitudes = amplitudes.data(); + ddsim->interface.getStateVectorFull(&ddsim->interface, &sv); + std::vector qubits; for (const auto& variable : assertion->getTargetQubits()) { qubits.push_back(variableToQubit(ddsim, variable)); } - Statevector sv; - sv.numQubits = 2; - sv.numStates = 4; - std::vector amplitudes(4); - sv.amplitudes = amplitudes.data(); - std::array qubitsStep = {0, 0}; + std::vector> densityMatrix( + sv.numStates, std::vector(sv.numStates, {0, 0})); + for (size_t i = 0; i < sv.numStates; i++) { + for (size_t j = 0; j < sv.numStates; j++) { + densityMatrix[i][j] = + complexMultiplication(amplitudes[i], complexConjugate(amplitudes[j])); + } + } - bool result = true; - for (size_t i = 0; i < qubits.size() && result; i++) { - for (size_t j = 0; j < qubits.size(); j++) { + for (const auto i : qubits) { + for (const auto j : qubits) { if (i == j) { continue; } - qubitsStep[0] = qubits[i]; - qubitsStep[1] = qubits[j]; - ddsim->interface.getStateVectorSub(&ddsim->interface, 2, - qubitsStep.data(), &sv); - if (!areQubitsEntangled(&sv)) { - result = false; - break; + if (!areQubitsEntangled(densityMatrix, i, j)) { + return false; } } } - - return result; + return true; } bool checkAssertionSuperposition( DDSimulationState* ddsim, std::unique_ptr& assertion) { std::vector qubits; + for (const auto& variable : assertion->getTargetQubits()) { qubits.push_back(variableToQubit(ddsim, variable)); } - Statevector sv; - sv.numQubits = qubits.size(); - sv.numStates = 1 << sv.numQubits; - std::vector amplitudes(sv.numStates); - sv.amplitudes = amplitudes.data(); - - ddsim->interface.getStateVectorSub(&ddsim->interface, sv.numQubits, - qubits.data(), &sv); - - int found = 0; - const double epsilon = 0.00000001; - for (size_t i = 0; i < sv.numStates && found < 2; i++) { - if (complexMagnitude(amplitudes[i]) > epsilon) { - found++; + Complex result; + bool firstFound = false; + std::vector bitstring; + for (size_t i = 0; + i < 1ULL << ddsim->interface.getNumQubits(&ddsim->interface); i++) { + ddsim->interface.getAmplitudeIndex(&ddsim->interface, i, &result); + if (complexMagnitude(result) > 0.00000001) { + if (!firstFound) { + firstFound = true; + bitstring = extractBits(qubits, i); + } else { + const std::vector compareBitstring = extractBits(qubits, i); + if (bitstring != compareBitstring) { + return true; + } + } } } - return found > 1; + return false; } bool checkAssertionEqualityStatevector( diff --git a/test/test_custom_code.cpp b/test/test_custom_code.cpp index 2bfb9b2..1de969c 100644 --- a/test/test_custom_code.cpp +++ b/test/test_custom_code.cpp @@ -108,3 +108,16 @@ TEST_F(CustomCodeTest, EqualityAssertion) { ASSERT_EQ(state->runAll(state, &numErrors), OK); ASSERT_EQ(numErrors, 0); } + +TEST_F(CustomCodeTest, DestructiveInterference) { + loadCode(3, 0, + "x q[0];" + "h q[0];" + "h q[1];" + "cx q[1], q[2];" + "assert-sup q[1], q[2];" + "assert-ent q[1], q[2];"); + size_t numErrors = 0; + ASSERT_EQ(state->runAll(state, &numErrors), OK); + ASSERT_EQ(numErrors, 0); +}