From 5009f032cb761a52a1954567e12526bb36bb6968 Mon Sep 17 00:00:00 2001 From: DRovara Date: Tue, 27 Aug 2024 15:10:45 +0200 Subject: [PATCH] fix: :bug: Implement correct substate vector equality assertions Now, we check if a given substate is separable, and if it is, we trace out unwanted qubits instead of performing a projection --- cmake/ExternalDependencies.cmake | 6 + include/backend/dd/DDSimDebug.hpp | 5 + src/backend/dd/DDSimDebug.cpp | 195 +++++++++++++--------- src/backend/dd/DDSimDiagnostics.cpp | 44 +++-- test/python/resources/bindings/jumps.qasm | 7 +- test/python/test_python_bindings.py | 18 +- test/test_custom_code.cpp | 46 +++++ test/test_data_retrieval.cpp | 2 +- 8 files changed, 212 insertions(+), 111 deletions(-) diff --git a/cmake/ExternalDependencies.cmake b/cmake/ExternalDependencies.cmake index 4d676f8..ec2e485 100644 --- a/cmake/ExternalDependencies.cmake +++ b/cmake/ExternalDependencies.cmake @@ -60,6 +60,9 @@ if(CMAKE_VERSION VERSION_GREATER_EQUAL 3.24) set(EIGEN_BUILD_TESTING OFF CACHE BOOL "Disable testing for Eigen") + set(BUILD_TESTING + OFF + CACHE BOOL "Disable general testing") set(EIGEN_BUILD_DOC OFF CACHE BOOL "Disable documentation build for Eigen") @@ -75,6 +78,9 @@ else() set(EIGEN_BUILD_TESTING OFF CACHE BOOL "Disable testing for Eigen") + set(BUILD_TESTING + OFF + CACHE BOOL "Disable general testing") set(EIGEN_BUILD_DOC OFF CACHE BOOL "Disable documentation build for Eigen") diff --git a/include/backend/dd/DDSimDebug.hpp b/include/backend/dd/DDSimDebug.hpp index 508982b..3696a5e 100644 --- a/include/backend/dd/DDSimDebug.hpp +++ b/include/backend/dd/DDSimDebug.hpp @@ -127,3 +127,8 @@ bool checkAssertion(DDSimulationState* ddsim, std::unique_ptr& assertion); std::string getClassicalBitName(DDSimulationState* ddsim, size_t index); size_t variableToQubit(DDSimulationState* ddsim, const std::string& variable); +bool isSubStateVectorLegal(const Statevector& full, + std::vector& targetQubits); +std::vector> +getPartialTraceFromStateVector(const Statevector& sv, + const std::vector& traceOut); diff --git a/src/backend/dd/DDSimDebug.cpp b/src/backend/dd/DDSimDebug.cpp index cab289e..65faf20 100644 --- a/src/backend/dd/DDSimDebug.cpp +++ b/src/backend/dd/DDSimDebug.cpp @@ -300,12 +300,17 @@ Result ddsimStepForward(SimulationState* self) { // - SIMULATE: run the corresponding operation on the DD backend. if (ddsim->instructionTypes[currentInstruction] == ASSERTION) { auto& assertion = ddsim->assertionInstructions[currentInstruction]; - const auto failed = !checkAssertion(ddsim, assertion); - if (failed && ddsim->lastFailedAssertion != currentInstruction) { - ddsim->lastFailedAssertion = currentInstruction; - self->stepBackward(self); + try { + const auto failed = !checkAssertion(ddsim, assertion); + if (failed && ddsim->lastFailedAssertion != currentInstruction) { + ddsim->lastFailedAssertion = currentInstruction; + self->stepBackward(self); + } + return OK; + } catch (const std::exception& e) { + std::cerr << e.what() << "\n"; + return ERROR; } - return OK; } ddsim->lastFailedAssertion = -1ULL; @@ -693,22 +698,64 @@ Result ddsimGetStateVectorSub(SimulationState* self, size_t subStateSize, std::vector amplitudes(fullState.numStates); const Span outAmplitudes(output->amplitudes, output->numStates); const Span qubitsSpan(qubits, subStateSize); + + std::vector targetQubits; + for (size_t i = 0; i < subStateSize; i++) { + targetQubits.push_back(qubitsSpan[i]); + } + fullState.amplitudes = amplitudes.data(); self->getStateVectorFull(self, &fullState); - for (size_t i = 0; i < outAmplitudes.size(); i++) { - outAmplitudes[i].real = 0; - outAmplitudes[i].imaginary = 0; + if (!isSubStateVectorLegal(fullState, targetQubits)) { + return ERROR; + } + + std::vector otherQubits; + for (size_t i = 0; i < fullState.numQubits; i++) { + if (std::find(targetQubits.begin(), targetQubits.end(), i) == + targetQubits.end()) { + otherQubits.push_back(i); + } + } + + const auto traced = getPartialTraceFromStateVector(fullState, otherQubits); + + // Create Eigen3 Matrix + Eigen::Matrix, Eigen::Dynamic, Eigen::Dynamic> mat( + static_cast(traced.size()), static_cast(traced.size())); + for (size_t i = 0; i < traced.size(); i++) { + for (size_t j = 0; j < traced.size(); j++) { + mat(static_cast(i), static_cast(j)) = { + traced[i][j].real, traced[i][j].imaginary}; + } } - for (size_t i = 0; i < fullState.numStates; i++) { - size_t outputIndex = 0; - for (size_t j = 0; j < subStateSize; j++) { - outputIndex |= ((i >> qubitsSpan[j]) & 1) << j; + const Eigen::ComplexEigenSolver< + Eigen::Matrix, Eigen::Dynamic, Eigen::Dynamic>> + solver(mat); + const auto& vectors = solver.eigenvectors(); + const auto& values = solver.eigenvalues(); + const auto epsilon = 0.000001; + int index = -1; + for (int i = 0; i < values.size(); i++) { + if (values[i].imag() < -epsilon || values[i].imag() > epsilon) { + continue; + } + if (values[i].real() - 1 < -epsilon || values[i].real() - 1 > epsilon) { + continue; } - outAmplitudes[outputIndex].real += amplitudes[i].real; - outAmplitudes[outputIndex].imaginary += amplitudes[i].imaginary; + index = i; + } + + if (index == -1) { + return ERROR; + } + + for (size_t i = 0; i < traced.size(); i++) { + outAmplitudes[i] = {vectors(static_cast(i), index).real(), + vectors(static_cast(i), index).imag()}; } return OK; @@ -996,19 +1043,23 @@ bool areQubitsEntangled(std::vector>& densityMatrix, 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) { +std::vector> +getPartialTraceFromStateVector(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 split1 = splitBitString(i, sv.numQubits, traceOut); + const auto split2 = splitBitString(j, sv.numQubits, traceOut); + if (split1.first != split2.first) { continue; } - const auto product = complexMultiplication(sv.amplitudes[i], -sv.amplitudes[j]); const auto row = i - traced1; const auto col = j - traced2; + const auto product = + complexMultiplication(sv.amplitudes[i], sv.amplitudes[j]); + const auto row = split1.second; + const auto col = split2.second; traceMatrix[row][col] = complexAddition(traceMatrix[row][col], product); } } @@ -1017,65 +1068,36 @@ sv.amplitudes[j]); const auto row = i - traced1; const auto col = j - traced2; 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])); + 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; + return runningSum; } -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); +bool partialTraceIsPure(const Statevector& sv, + const std::vector& traceOut) { + const auto traceMatrix = getPartialTraceFromStateVector(sv, traceOut); 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; - const bool canBe01 = complexMagnitude(amplitudes[1]) > epsilon; - const bool canBe10 = complexMagnitude(amplitudes[2]) > epsilon; - const bool canBe11 = complexMagnitude(amplitudes[3]) > epsilon; - - const int nonZeroCount = (canBe00 ? 1 : 0) + (canBe01 ? 1 : 0) + - (canBe10 ? 1 : 0) + (canBe11 ? 1 : 0); - - if (nonZeroCount == 4) { - const auto c1 = complexDivision(amplitudes[0], amplitudes[2]); - const auto c2 = complexDivision(amplitudes[1], amplitudes[3]); - return !areComplexEqual(c1, c2); - } - - return (canBe00 && canBe11 && !(canBe01 && canBe10)) || - (canBe01 && canBe10 && !(canBe00 && canBe11)); + const double epsilon = 0.00000001; + return trace.imaginary < epsilon && trace.imaginary > -epsilon && + (trace.real - 1) < epsilon && (trace.real - 1) > -epsilon; } -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); +bool isSubStateVectorLegal(const Statevector& full, + std::vector& targetQubits) { + const auto numQubits = full.numQubits; + std::vector ignored; + for (size_t i = 0; i < numQubits; i++) { + if (std::find(targetQubits.begin(), targetQubits.end(), i) == + targetQubits.end()) { + ignored.push_back(i); } } - return partialTraceIsPure(sv, toTraceOut); -}*/ + return partialTraceIsPure(full, ignored); +} bool checkAssertionEntangled( DDSimulationState* ddsim, @@ -1159,8 +1181,11 @@ bool checkAssertionEqualityStatevector( std::vector amplitudes(sv.numStates); sv.amplitudes = amplitudes.data(); - ddsim->interface.getStateVectorSub(&ddsim->interface, sv.numQubits, - qubits.data(), &sv); + if (ddsim->interface.getStateVectorSub(&ddsim->interface, sv.numQubits, + qubits.data(), &sv) == ERROR) { + throw std::runtime_error( + "Equality assertion on entangled sub-state is not allowed."); + } const double similarityThreshold = assertion->getSimilarityThreshold(); @@ -1172,6 +1197,11 @@ bool checkAssertionEqualityStatevector( bool checkAssertionEqualityCircuit( DDSimulationState* ddsim, std::unique_ptr& assertion) { + std::vector qubits; + for (const auto& variable : assertion->getTargetQubits()) { + qubits.push_back(variableToQubit(ddsim, variable)); + } + DDSimulationState secondSimulation; createDDSimulationState(&secondSimulation); secondSimulation.interface.loadCode(&secondSimulation.interface, @@ -1193,17 +1223,16 @@ bool checkAssertionEqualityCircuit( &sv2); destroyDDSimulationState(&secondSimulation); - 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); + if (ddsim->interface.getStateVectorSub(&ddsim->interface, sv.numQubits, + qubits.data(), &sv) == ERROR) { + throw std::runtime_error( + "Equality assertion on entangled sub-state is not allowed."); + } const double similarityThreshold = assertion->getSimilarityThreshold(); diff --git a/src/backend/dd/DDSimDiagnostics.cpp b/src/backend/dd/DDSimDiagnostics.cpp index 2daa119..7619adb 100644 --- a/src/backend/dd/DDSimDiagnostics.cpp +++ b/src/backend/dd/DDSimDiagnostics.cpp @@ -219,35 +219,49 @@ size_t tryFindZeroControls(DDDiagnostics* diagnostics, size_t instruction, return index; } +bool isAlwaysZero(const Statevector& sv, size_t qubit, bool checkOne = false) { + const auto epsilon = 1e-10; + + const Span amplitudes(sv.amplitudes, sv.numStates); + + for (size_t i = 0; i < sv.numStates; i++) { + if (((i & (1ULL << qubit)) == 0ULL && !checkOne) || + ((i & (1ULL << qubit)) != 0ULL && checkOne)) { + continue; + } + if (amplitudes[i].real > epsilon || amplitudes[i].real < -epsilon || + amplitudes[i].imaginary > epsilon || + amplitudes[i].imaginary < -epsilon) { + return false; + } + } + + return true; +} + void dddiagnosticsOnStepForward(DDDiagnostics* diagnostics, size_t instruction) { auto* ddsim = diagnostics->simulationState; if (ddsim->instructionTypes[instruction] != SIMULATE) { return; } + const auto numQubits = ddsim->interface.getNumQubits(&ddsim->interface); const auto& op = (*ddsim->iterator); const auto& controls = op->getControls(); - std::vector amplitudes(2); - Statevector sv{1, 2, amplitudes.data()}; - std::vector qubits{0}; - const double epsilon = 1e-10; + std::vector amplitudes(2ULL << numQubits); + Statevector sv{numQubits, 2ULL << numQubits, amplitudes.data()}; + ddsim->interface.getStateVectorFull(&ddsim->interface, &sv); for (const auto& control : controls) { const auto pos = control.type == qc::Control::Type::Pos; const auto qubit = control.qubit; - qubits[0] = qubit; - if (ddsim->interface.getStateVectorSub(&ddsim->interface, 1, qubits.data(), - &sv) == OK) { - const auto amplitude = amplitudes[pos ? 1 : 0]; - if (amplitude.real < epsilon && amplitude.real > -epsilon && - amplitude.imaginary < epsilon && amplitude.imaginary > -epsilon) { - if (diagnostics->zeroControls.find(instruction) == - diagnostics->zeroControls.end()) { - diagnostics->zeroControls[instruction] = std::set(); - } - diagnostics->zeroControls[instruction].insert(qubit); + if (isAlwaysZero(sv, qubit, !pos)) { + if (diagnostics->zeroControls.find(instruction) == + diagnostics->zeroControls.end()) { + diagnostics->zeroControls[instruction] = std::set(); } + diagnostics->zeroControls[instruction].insert(qubit); } } } diff --git a/test/python/resources/bindings/jumps.qasm b/test/python/resources/bindings/jumps.qasm index 1d0dbfd..26b02d3 100644 --- a/test/python/resources/bindings/jumps.qasm +++ b/test/python/resources/bindings/jumps.qasm @@ -23,10 +23,11 @@ gate create_ghz q0, q1, q2 { create_ghz q[0], q[1], q[2]; -assert-eq q[0], q[1] { - qreg q[2]; +assert-eq q[0], q[1], q[2] { + qreg q[3]; h q[0]; cx q[0], q[1]; + cx q[0], q[2]; } -assert-eq 0.9, q[0] { 0.7, 0.7 } +assert-eq 0.9, q[0], q[1], q[2] { 0.7, 0, 0, 0, 0, 0, 0, 0.7 } diff --git a/test/python/test_python_bindings.py b/test/python/test_python_bindings.py index 25bc182..289960a 100644 --- a/test/python/test_python_bindings.py +++ b/test/python/test_python_bindings.py @@ -209,15 +209,15 @@ def test_access_state(simulation_instance_jumps: SimulationInstance) -> None: assert abs(c.imaginary) < 1e-6 assert abs(c.real) < 1e-6 - sv = simulation_state.get_state_vector_sub([0, 2]) - assert sv.num_qubits == 2 - assert sv.num_states == 4 - c = sv.amplitudes[3] - assert abs(c.imaginary) < 1e-6 - assert abs(c.real - 1 / (2**0.5)) < 1e-6 - c = sv.amplitudes[2] - assert abs(c.imaginary) < 1e-6 - assert abs(c.real) < 1e-6 + +def test_get_state_vector_sub(simulation_instance_classical: SimulationInstance) -> None: + """Tests the `get_state_vector_sub()` method.""" + (simulation_state, _state_id) = simulation_instance_classical + simulation_state.set_breakpoint(170) + simulation_state.run_simulation() + assert simulation_state.get_current_instruction() == 10 + sv = simulation_state.get_state_vector_sub([0, 1]) + assert sv.amplitudes[0].real == 1 or sv.amplitudes[-1].real == 1 def test_classical_get(simulation_instance_classical: SimulationInstance) -> None: diff --git a/test/test_custom_code.cpp b/test/test_custom_code.cpp index 1de969c..dbeccad 100644 --- a/test/test_custom_code.cpp +++ b/test/test_custom_code.cpp @@ -121,3 +121,49 @@ TEST_F(CustomCodeTest, DestructiveInterference) { ASSERT_EQ(state->runAll(state, &numErrors), OK); ASSERT_EQ(numErrors, 0); } + +TEST_F(CustomCodeTest, IllegalSubstateSVEqualityAssertion) { + loadCode(3, 0, + "x q[0];" + "h q[0];" + "h q[1];" + "cx q[1], q[2];" + "assert-eq 0.9, q[0], q[1] { 0.5, 0.5, 0.5, 0.5 }"); + size_t numErrors = 0; + ASSERT_EQ(state->runAll(state, &numErrors), ERROR); +} + +TEST_F(CustomCodeTest, LegalSubstateSVEqualityAssertion) { + loadCode(3, 0, + "x q[0];" + "h q[0];" + "h q[1];" + "cx q[1], q[2];" + "assert-eq 0.9, q[1], q[2] { 0.707, 0, 0, 0.707 }"); + size_t numErrors = 0; + ASSERT_EQ(state->runAll(state, &numErrors), OK); + ASSERT_EQ(numErrors, 0); +} + +TEST_F(CustomCodeTest, IllegalSubstateCircuitEqualityAssertion) { + loadCode(3, 0, + "x q[0];" + "h q[0];" + "h q[1];" + "cx q[1], q[2];" + "assert-eq 0.9, q[0], q[1] { qreg q[2]; h q[0]; h q[1]; }"); + size_t numErrors = 0; + ASSERT_EQ(state->runAll(state, &numErrors), ERROR); +} + +TEST_F(CustomCodeTest, LegalSubstateCircuitEqualityAssertion) { + loadCode(3, 0, + "x q[0];" + "h q[0];" + "h q[1];" + "cx q[1], q[2];" + "assert-eq 0.9, q[2], q[1] { qreg q[2]; h q[0]; cx q[0], q[1]; }"); + size_t numErrors = 0; + ASSERT_EQ(state->runAll(state, &numErrors), OK); + ASSERT_EQ(numErrors, 0); +} diff --git a/test/test_data_retrieval.cpp b/test/test_data_retrieval.cpp index 1490cc1..7fe1557 100644 --- a/test/test_data_retrieval.cpp +++ b/test/test_data_retrieval.cpp @@ -169,7 +169,7 @@ TEST_F(DataRetrievalTest, GetStateVectorSub) { qubits[0] = 1; ASSERT_EQ(state->getStateVectorSub(state, 2, qubits.data(), &sv), OK); ASSERT_TRUE(complexEquality(amplitudes[0], 0.0, 0.0)); - ASSERT_TRUE(complexEquality(amplitudes[1], 0.0, + ASSERT_TRUE(complexEquality(amplitudes[1], 1.0, 0.0)); // 0 because of destructive interference ASSERT_TRUE(complexEquality(amplitudes[2], 0.0, 0.0)); ASSERT_TRUE(complexEquality(amplitudes[3], 0.0, 0.0));