Skip to content

Commit

Permalink
Started to add tests for the state of the MPS simulator by executing …
Browse files Browse the repository at this point in the history 
…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.
  • Loading branch information
@aromanro
aromanro committed Aug 24, 2024
1 parent d313885 commit cfcd7f1
Show file tree
Hide file tree
Showing 2 changed files with 155 additions and 45 deletions.
15 changes: 11 additions & 4 deletions QCSim/MPSSimulator.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<int>(qubit) - static_cast<int>(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");


Expand Down Expand Up @@ -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)
Expand All @@ -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 {};
}
Expand Down Expand Up @@ -359,14 +364,15 @@ namespace QC {

SVD.compute(thetaMatrix, Eigen::DecompositionOptions::ComputeThinU | Eigen::DecompositionOptions::ComputeThinV);


const MatrixClass& UmatrixFull = SVD.matrixU();
const MatrixClass& VmatrixFull = SVD.matrixV();

const Eigen::VectorXd& SvaluesFull = SVD.singularValues();

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<long long>(szm, UmatrixFull.cols());
const long long sz = limitSize ? std::min<long long>(chi, szm) : szm;
Expand Down Expand Up @@ -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;
Expand Down
185 changes: 144 additions & 41 deletions QCSim/MPSSimulatorTests.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,71 +4,174 @@

#include "Tests.h"

#include "QubitRegister.h"
#include "MPSSimulator.h"

bool StateSimulationTest()
#define _USE_MATH_DEFINES
#include <math.h>


void FillOneQubitGates(std::vector<std::shared_ptr<QC::Gates::QuantumGateWithOp<>>>& gates)
{
gates.emplace_back(std::make_shared<QC::Gates::HadamardGate<>>());
gates.emplace_back(std::make_shared<QC::Gates::HyGate<>>());
gates.emplace_back(std::make_shared<QC::Gates::PhaseGate<>>());
gates.emplace_back(std::make_shared<QC::Gates::PhaseShiftGate<>>(0.5 * M_PI));
gates.emplace_back(std::make_shared<QC::Gates::PauliXGate<>>());
gates.emplace_back(std::make_shared<QC::Gates::PauliYGate<>>());
gates.emplace_back(std::make_shared<QC::Gates::PauliZGate<>>());
gates.emplace_back(std::make_shared<QC::Gates::SigmaPlusGate<>>());
gates.emplace_back(std::make_shared<QC::Gates::SigmaMinusGate<>>());
gates.emplace_back(std::make_shared<QC::Gates::SquareRootNOTGate<>>());
gates.emplace_back(std::make_shared<QC::Gates::SquareRootNOTDagGate<>>());
gates.emplace_back(std::make_shared<QC::Gates::SplitterGate<>>());
gates.emplace_back(std::make_shared<QC::Gates::RxGate<>>(M_PI / 4));
gates.emplace_back(std::make_shared<QC::Gates::RyGate<>>(M_PI / 4));
gates.emplace_back(std::make_shared<QC::Gates::RzGate<>>(M_PI / 4));
gates.emplace_back(std::make_shared<QC::Gates::UGate<>>(M_PI / 4, M_PI / 8));
}

void FillTwoQubitGates(std::vector<std::shared_ptr<QC::Gates::QuantumGateWithOp<>>>& 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<QC::Gates::SwapGate<>>());
gates.emplace_back(std::make_shared<QC::Gates::iSwapGate<>>());
gates.emplace_back(std::make_shared<QC::Gates::DecrementGate<>>());
gates.emplace_back(std::make_shared<QC::Gates::CNOTGate<>>());
gates.emplace_back(std::make_shared<QC::Gates::ControlledYGate<>>());
gates.emplace_back(std::make_shared<QC::Gates::ControlledZGate<>>());
gates.emplace_back(std::make_shared<QC::Gates::ControlledHadamardGate<>>());
gates.emplace_back(std::make_shared<QC::Gates::ControlledSquareRootNOTGate<>>());
gates.emplace_back(std::make_shared<QC::Gates::ControlledSquareRootNOTDagGate<>>());
gates.emplace_back(std::make_shared<QC::Gates::ControlledPhaseGate<>>());
gates.emplace_back(std::make_shared<QC::Gates::ControlledPhaseShiftGate<>>(M_PI / 3));
gates.emplace_back(std::make_shared<QC::Gates::ControlledUGate<>>(M_PI / 4, M_PI / 3));
gates.emplace_back(std::make_shared<QC::Gates::ControlledRxGate<>>(M_PI / 5));
gates.emplace_back(std::make_shared<QC::Gates::ControlledRyGate<>>(M_PI / 3));
gates.emplace_back(std::make_shared<QC::Gates::ControlledRzGate<>>(M_PI / 4));
}

bool OneQubitGatesTest()
{
std::cout << "\nMPS simulator state test for one qubit gates" << std::endl;

std::vector<std::shared_ptr<QC::Gates::QuantumGateWithOp<>>> gates;
FillOneQubitGates(gates);

std::uniform_int_distribution nrGatesDistr(25, 50);
std::uniform_int_distribution gateDistr(0, static_cast<int>(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<std::shared_ptr<QC::Gates::QuantumGateWithOp<>>> 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<int>(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();
}

0 comments on commit cfcd7f1

Please sign in to comment.