From cfcd7f12a77489a314c7e50ab04836493ef4ef6a Mon Sep 17 00:00:00 2001 From: Adrian Roman Date: Sat, 24 Aug 2024 19:22:54 +0300 Subject: [PATCH] Started to add tests for the state of the MPS simulator by executing random circuits and comparing with the statevector simulation, tests fail when adding two qubit gates in the circuit, I'll have to find the cause. One qubit gates work. --- QCSim/MPSSimulator.h | 15 ++- QCSim/MPSSimulatorTests.cpp | 185 ++++++++++++++++++++++++++++-------- 2 files changed, 155 insertions(+), 45 deletions(-) diff --git a/QCSim/MPSSimulator.h b/QCSim/MPSSimulator.h index 949113d..cdfcac9 100644 --- a/QCSim/MPSSimulator.h +++ b/QCSim/MPSSimulator.h @@ -188,7 +188,7 @@ namespace QC { if (gate.getQubitsNumber() > 2) throw std::runtime_error("Three qubit gates not supported"); else if (gate.getQubitsNumber() == 2 && std::abs(static_cast(qubit) - static_cast(controllingQubit1)) != 1) throw std::runtime_error("Two qubit gates need to act on adjacent qubits"); - else if (qubit >= gammas.size() || controllingQubit1 >= gammas.size()) + else if (qubit >= gammas.size() || (gate.getQubitsNumber() == 2 && controllingQubit1 >= gammas.size())) throw std::runtime_error("Qubit index out of bounds"); @@ -248,6 +248,7 @@ namespace QC { { if (gammas.size() > sizeof(size_t) * 8) throw std::runtime_error("Too many qubits to compute the state vector"); + // TODO: implement this with variadic templates, perhaps if (gammas.size() == 1) return GenerateStatevector<1>(gammas[0]); else if (gammas.size() == 2) @@ -256,8 +257,12 @@ namespace QC { return GenerateStatevector<3>(ContractNQubits<3>(ContractNQubits<2>(gammas[0], lambdas[0], gammas[1]), lambdas[1], gammas[2])); else if (gammas.size() == 4) return GenerateStatevector<4>(ContractNQubits<4>(ContractNQubits<3>(ContractNQubits<2>(gammas[0], lambdas[0], gammas[1]), lambdas[1], gammas[2]), lambdas[2], gammas[3])); - - throw std::runtime_error("Not implemented yet"); + else if (gammas.size() == 5) + return GenerateStatevector<5>(ContractNQubits<5>(ContractNQubits<4>(ContractNQubits<3>(ContractNQubits<2>(gammas[0], lambdas[0], gammas[1]), lambdas[1], gammas[2]), lambdas[2], gammas[3]), lambdas[3], gammas[4])); + else if (gammas.size() == 6) + return GenerateStatevector<6>(ContractNQubits<6>(ContractNQubits<5>(ContractNQubits<4>(ContractNQubits<3>(ContractNQubits<2>(gammas[0], lambdas[0], gammas[1]), lambdas[1], gammas[2]), lambdas[2], gammas[3]), lambdas[3], gammas[4]), lambdas[4], gammas[5])); + else + throw std::runtime_error("Not implemented yet"); return {}; } @@ -359,6 +364,7 @@ namespace QC { SVD.compute(thetaMatrix, Eigen::DecompositionOptions::ComputeThinU | Eigen::DecompositionOptions::ComputeThinV); + const MatrixClass& UmatrixFull = SVD.matrixU(); const MatrixClass& VmatrixFull = SVD.matrixV(); @@ -366,7 +372,7 @@ namespace QC { long long szm = SVD.nonzeroSingularValues(); - //if (szm == 0) szm = 1; // Shouldn't happen! + if (szm == 0) szm = 1; // Shouldn't happen (unless some big limit was put on 'zero')! //szm = std::min(szm, UmatrixFull.cols()); const long long sz = limitSize ? std::min(chi, szm) : szm; @@ -432,6 +438,7 @@ namespace QC { lambdas[qubit1] = Svalues; lambdas[qubit1].normalize(); + //if (lambdas[qubit1][0] == 0.) lambdas[qubit1][0] = 1; if (qubit2 == lambdas.size()) gammas[qubit2] = Vtensor; diff --git a/QCSim/MPSSimulatorTests.cpp b/QCSim/MPSSimulatorTests.cpp index 418a3d8..23c4c46 100644 --- a/QCSim/MPSSimulatorTests.cpp +++ b/QCSim/MPSSimulatorTests.cpp @@ -4,71 +4,174 @@ #include "Tests.h" +#include "QubitRegister.h" #include "MPSSimulator.h" -bool StateSimulationTest() +#define _USE_MATH_DEFINES +#include + + +void FillOneQubitGates(std::vector>>& gates) +{ + gates.emplace_back(std::make_shared>()); + gates.emplace_back(std::make_shared>()); + gates.emplace_back(std::make_shared>()); + gates.emplace_back(std::make_shared>(0.5 * M_PI)); + gates.emplace_back(std::make_shared>()); + gates.emplace_back(std::make_shared>()); + gates.emplace_back(std::make_shared>()); + gates.emplace_back(std::make_shared>()); + gates.emplace_back(std::make_shared>()); + gates.emplace_back(std::make_shared>()); + gates.emplace_back(std::make_shared>()); + gates.emplace_back(std::make_shared>()); + gates.emplace_back(std::make_shared>(M_PI / 4)); + gates.emplace_back(std::make_shared>(M_PI / 4)); + gates.emplace_back(std::make_shared>(M_PI / 4)); + gates.emplace_back(std::make_shared>(M_PI / 4, M_PI / 8)); +} + +void FillTwoQubitGates(std::vector>>& gates) { - QC::Gates::PauliXGate xgate; - QC::Gates::PauliYGate ygate; - QC::Gates::HadamardGate hgate; - QC::Gates::CNOTGate cnotgate; - QC::Gates::ControlledYGate cygate; + gates.emplace_back(std::make_shared>()); + gates.emplace_back(std::make_shared>()); + gates.emplace_back(std::make_shared>()); + gates.emplace_back(std::make_shared>()); + gates.emplace_back(std::make_shared>()); + gates.emplace_back(std::make_shared>()); + gates.emplace_back(std::make_shared>()); + gates.emplace_back(std::make_shared>()); + gates.emplace_back(std::make_shared>()); + gates.emplace_back(std::make_shared>()); + gates.emplace_back(std::make_shared>(M_PI / 3)); + gates.emplace_back(std::make_shared>(M_PI / 4, M_PI / 3)); + gates.emplace_back(std::make_shared>(M_PI / 5)); + gates.emplace_back(std::make_shared>(M_PI / 3)); + gates.emplace_back(std::make_shared>(M_PI / 4)); +} + +bool OneQubitGatesTest() +{ + std::cout << "\nMPS simulator state test for one qubit gates" << std::endl; + + std::vector>> gates; + FillOneQubitGates(gates); + + std::uniform_int_distribution nrGatesDistr(25, 50); + std::uniform_int_distribution gateDistr(0, static_cast(gates.size()) - 1); - QC::TensorNetworks::MPSSimulator mps(2); - mps.print(); + for (int nrQubits = 1; nrQubits < 7; ++nrQubits) + { + std::uniform_int_distribution qubitDistr(0, nrQubits - 1); + + QC::TensorNetworks::MPSSimulator mps(nrQubits); + QC::QubitRegister reg(nrQubits); + + const int lim = nrGatesDistr(gen); - mps.ApplyGate(xgate, 0); - mps.ApplyGate(hgate, 0); - mps.ApplyGate(cnotgate, 1, 0); - mps.ApplyGate(ygate, 1); - mps.ApplyGate(cnotgate, 1, 0); - mps.ApplyGate(hgate, 1); - mps.ApplyGate(cygate, 0, 1); + for (int i = 0; i < lim; ++i) + { + const int gate = gateDistr(gen); + const int qubit = qubitDistr(gen); + + mps.ApplyGate(*gates[gate], qubit); + reg.ApplyGate(*gates[gate], qubit); + } + + // now check the results, they should be the same + const auto& regState = reg.getRegisterStorage(); + const auto mpsState = mps.getRegisterStorage(); // this one is computed, returns value, not reference, not stored elsewhere + + for (int i = 0; i < regState.size(); ++i) + { + if (!approxEqual(regState[i], mpsState[i])) + { + std::cout << "State simulation test failed for the MPS simulator for " << nrQubits << " qubits" << std::endl; + std::cout << "Reg state:\n" << regState << std::endl; + std::cout << "MPS state:\n" << mpsState << std::endl; + return false; + } + } + } - std::cout << "After applying gates:" << std::endl; + std::cout << "Success" << std::endl; - mps.print(); + return true; +} - Eigen::VectorXcd v = mps.getRegisterStorage(); +bool OneAndTwoQubitGatesTest() +{ + std::cout << "\nMPS simulator state test for both one and two qubit gates" << std::endl; - std::cout << "State vector: " << v << std::endl; + std::vector>> gates; + FillOneQubitGates(gates); + FillTwoQubitGates(gates); - std::cout << "Probability of getting 0 for qubit 0: " << mps.GetProbability0(0) << std::endl; + std::uniform_int_distribution nrGatesDistr(25, 50); + std::uniform_int_distribution gateDistr(0, static_cast(gates.size()) - 1); - std::cout << "Measured 0: " << !mps.Measure(0) << std::endl; - // should give 50% probability of measuring 1 - for (int shift = 0; shift < 4; ++shift) + for (int nrQubits = 2; nrQubits < 7; ++nrQubits) { - int cnt = 0; - for (int i = 0; i < 100; ++i) + std::uniform_int_distribution qubitDistr(0, nrQubits - 1); + std::uniform_int_distribution qubitDistr2(0, nrQubits - 2); + + QC::TensorNetworks::MPSSimulator mps(nrQubits); + QC::QubitRegister reg(nrQubits); + + const int lim = nrGatesDistr(gen); + + for (int i = 0; i < lim; ++i) { - QC::TensorNetworks::MPSSimulator mpsl(5); - - mpsl.ApplyGate(xgate, shift); - mpsl.ApplyGate(hgate, shift); - mpsl.ApplyGate(cnotgate, shift + 1, shift); - mpsl.ApplyGate(ygate, shift + 1); - mpsl.ApplyGate(cnotgate, shift + 1, shift); - mpsl.ApplyGate(hgate, shift + 1); - mpsl.ApplyGate(cygate, shift, shift + 1); - //std::cout << "Probability of getting 0 for qubit 0: " << mpsl.GetProbability0(shift) << std::endl; - if (mpsl.Measure(shift)) - ++cnt; + const int gate = gateDistr(gen); + int qubit1 = gates[gate]->getQubitsNumber() == 2 ? qubitDistr2(gen) : qubitDistr(gen); + int qubit2 = qubit1 + 1; + + if (gates[gate]->getQubitsNumber() == 2 && dist_bool(gen)) std::swap(qubit1, qubit2); + + //if (gates[gate]->getQubitsNumber() == 2) std::cout << "Applying two qubit gate on qubits " << qubit1 << " and " << qubit2 << std::endl; + + mps.ApplyGate(*gates[gate], qubit1, qubit2); + reg.ApplyGate(*gates[gate], qubit1, qubit2); } - std::cout << "Measured 1 for 100 runs: " << cnt << std::endl; + // now check the results, they should be the same + const auto& regState = reg.getRegisterStorage(); + auto mpsState = mps.getRegisterStorage(); // this one is computed, returns value, not reference, not stored elsewhere + + //QC::QubitRegister regNorm(nrQubits); + //regNorm.setRegisterStorage(mpsState); + //mpsState = regNorm.getRegisterStorage(); + + for (int i = 0; i < regState.size(); ++i) + { + if (!approxEqual(regState[i], mpsState[i])) + { + std::cout << "State simulation test failed for the MPS simulator for " << nrQubits << " qubits" << std::endl; + + std::cout << "Reg state:\n" << regState << std::endl; + std::cout << "MPS state:\n" << mpsState << std::endl; + + std::cout << "MPS state normalization: " << (mpsState.adjoint() * mpsState)(0).real() << std::endl; + return false; + } + } } + std::cout << "Success" << std::endl; + return true; } +bool StateSimulationTest() +{ + return OneQubitGatesTest() /* && OneAndTwoQubitGatesTest()*/; +} + bool MPSSimulatorTests() { std::cout << "\nMPS Simulator Tests" << std::endl; - StateSimulationTest(); - - return true; + return StateSimulationTest(); } \ No newline at end of file