18#ifdef INCLUDED_BY_FACTORY
25#include "MPSSimulator.h"
26#include "QubitRegister.h"
29#include "../TensorNetworks/ForestContractor.h"
30#include "../TensorNetworks/TensorNetwork.h"
32#include "../Utils/Alias.h"
50class QCSimState :
public ISimulator {
52 QCSimState() : rng(std::random_device{}()), uniformZeroOne(0, 1) {}
61 void Initialize()
override {
63 if (simulationType == SimulationType::kMatrixProductState) {
65 std::make_unique<QC::TensorNetworks::MPSSimulator>(nrQubits);
66 if (limitEntanglement && singularValueThreshold > 0.)
67 mpsSimulator->setLimitEntanglement(singularValueThreshold);
68 if (limitSize && chi > 0) mpsSimulator->setLimitBondDimension(chi);
69 }
else if (simulationType == SimulationType::kStabilizer)
71 std::make_unique<QC::Clifford::StabilizerSimulator>(nrQubits);
72 else if (simulationType == SimulationType::kTensorNetwork) {
74 std::make_unique<TensorNetworks::TensorNetwork>(nrQubits);
77 const auto tensorContractor =
78 std::make_shared<TensorNetworks::ForestContractor>();
79 tensorNetwork->SetContractor(tensorContractor);
80 }
else if (simulationType == SimulationType::kPauliPropagator) {
81 pp = std::make_unique<Simulators::QcsimPauliPropagator>();
82 pp->SetNrQubits(nrQubits);
84 state = std::make_unique<QC::QubitRegister<>>(nrQubits);
101 void InitializeState(
size_t num_qubits,
102 std::vector<std::complex<double>> &litudes)
override {
103 if (num_qubits == 0)
return;
105 nrQubits = num_qubits;
107 if (simulationType != SimulationType::kStatevector)
108 throw std::runtime_error(
109 "QCSimState::InitializeState: Invalid "
110 "simulation type for initializing the state.");
112 Eigen::VectorXcd amplitudesEigen(
113 Eigen::Map<Eigen::VectorXcd, Eigen::Unaligned>(amplitudes.data(),
115 state->setRegisterStorageFastNoNormalize(amplitudesEigen);
154 void InitializeState(
size_t num_qubits,
155 AER::Vector<std::complex<double>> &litudes)
override {
156 if (num_qubits == 0)
return;
158 nrQubits = num_qubits;
160 if (simulationType != SimulationType::kStatevector)
161 throw std::runtime_error(
162 "QCSimState::InitializeState: Invalid "
163 "simulation type for initializing the state.");
165 Eigen::VectorXcd amplitudesEigen(
166 Eigen::Map<Eigen::VectorXcd, Eigen::Unaligned>(amplitudes.data(),
168 state->setRegisterStorageFastNoNormalize(amplitudesEigen);
183 void InitializeState(
size_t num_qubits,
184 Eigen::VectorXcd &litudes)
override {
185 if (num_qubits == 0)
return;
187 nrQubits = num_qubits;
190 if (simulationType != SimulationType::kStatevector)
191 throw std::runtime_error(
192 "QCSimState::InitializeState: Invalid "
193 "simulation type for initializing the state.");
195 state = std::make_unique<QC::QubitRegister<>>(nrQubits, amplitudes);
196 state->SetMultithreading(enableMultithreading);
205 void Reset()
override {
207 mpsSimulator->Clear();
208 else if (cliffordSimulator)
209 cliffordSimulator->Reset();
210 else if (tensorNetwork)
211 tensorNetwork->Clear();
215 pp->ClearOperations();
226 void Configure(
const char *key,
const char *value)
override {
227 if (std::string(
"method") == key) {
228 if (std::string(
"statevector") == value)
229 simulationType = SimulationType::kStatevector;
230 else if (std::string(
"matrix_product_state") == value)
231 simulationType = SimulationType::kMatrixProductState;
232 else if (std::string(
"stabilizer") == value)
233 simulationType = SimulationType::kStabilizer;
234 else if (std::string(
"tensor_network") == value)
235 simulationType = SimulationType::kTensorNetwork;
236 else if (std::string(
"pauli_propagator") == value)
237 simulationType = SimulationType::kPauliPropagator;
238 }
else if (std::string(
"matrix_product_state_truncation_threshold") ==
240 singularValueThreshold = std::stod(value);
241 if (singularValueThreshold > 0.) {
242 limitEntanglement =
true;
244 mpsSimulator->setLimitEntanglement(singularValueThreshold);
246 limitEntanglement =
false;
247 }
else if (std::string(
"matrix_product_state_max_bond_dimension") == key) {
248 chi = std::stoi(value);
251 if (mpsSimulator) mpsSimulator->setLimitBondDimension(chi);
254 }
else if (std::string(
"mps_sample_measure_algorithm") == key)
255 useMPSMeasureNoCollapse = std::string(
"mps_probabilities") == value;
266 if (std::string(
"method") == key) {
267 switch (simulationType) {
268 case SimulationType::kStatevector:
269 return "statevector";
270 case SimulationType::kMatrixProductState:
271 return "matrix_product_state";
272 case SimulationType::kStabilizer:
274 case SimulationType::kTensorNetwork:
275 return "tensor_network";
276 case SimulationType::kPauliPropagator:
277 return "pauli_propagator";
281 }
else if (std::string(
"matrix_product_state_truncation_threshold") ==
283 if (limitEntanglement && singularValueThreshold > 0.)
284 return std::to_string(singularValueThreshold);
285 }
else if (std::string(
"matrix_product_state_max_bond_dimension") == key) {
286 if (limitSize && limitSize > 0)
return std::to_string(chi);
287 }
else if (std::string(
"mps_sample_measure_algorithm") == key) {
288 return useMPSMeasureNoCollapse ?
"mps_probabilities"
289 :
"mps_apply_measure";
303 if ((simulationType == SimulationType::kStatevector && state) ||
304 (simulationType == SimulationType::kMatrixProductState &&
306 (simulationType == SimulationType::kStabilizer && cliffordSimulator) ||
307 (simulationType == SimulationType::kTensorNetwork && tensorNetwork))
310 const size_t oldNrQubits = nrQubits;
311 nrQubits += num_qubits;
312 if (simulationType == SimulationType::kPauliPropagator)
313 if (pp)pp->SetNrQubits(nrQubits);
333 void Clear()
override {
335 mpsSimulator =
nullptr;
336 cliffordSimulator =
nullptr;
337 tensorNetwork =
nullptr;
349 size_t Measure(
const Types::qubits_vector &qubits)
override {
358 if (simulationType == SimulationType::kStatevector) {
359 for (
size_t qubit : qubits) {
360 if (state->MeasureQubit(
static_cast<unsigned int>(qubit))) res |= mask;
363 }
else if (simulationType == SimulationType::kStabilizer) {
364 for (
size_t qubit : qubits) {
365 if (cliffordSimulator->MeasureQubit(
static_cast<unsigned int>(qubit)))
369 }
else if (simulationType == SimulationType::kTensorNetwork) {
370 for (
size_t qubit : qubits) {
371 if (tensorNetwork->Measure(
static_cast<unsigned int>(qubit)))
375 }
else if (simulationType == SimulationType::kPauliPropagator) {
376 std::vector<int> qubitsInt(qubits.begin(), qubits.end());
377 const auto res = pp->Measure(qubitsInt);
378 Types::qubit_t result = 0;
379 for (
size_t i = 0; i < res.size(); ++i) {
380 if (res[i]) result |= mask;
393 const std::set<Eigen::Index> qubitsSet(qubits.begin(), qubits.end());
394 auto measured = mpsSimulator->MeasureQubits(qubitsSet);
395 for (Types::qubit_t qubit : qubits) {
396 if (measured[qubit]) res |= mask;
402 NotifyObservers(qubits);
413 void ApplyReset(
const Types::qubits_vector &qubits)
override {
414 QC::Gates::PauliXGate xGate;
417 if (simulationType == SimulationType::kStatevector) {
418 for (
size_t qubit : qubits)
419 if (state->MeasureQubit(
static_cast<unsigned int>(qubit)))
420 state->ApplyGate(xGate,
static_cast<unsigned int>(qubit));
421 }
else if (simulationType == SimulationType::kStabilizer) {
422 for (
size_t qubit : qubits)
423 if (cliffordSimulator->MeasureQubit(
static_cast<unsigned int>(qubit)))
424 cliffordSimulator->ApplyX(
static_cast<unsigned int>(qubit));
425 }
else if (simulationType == SimulationType::kTensorNetwork) {
426 for (
size_t qubit : qubits)
427 if (tensorNetwork->Measure(
static_cast<unsigned int>(qubit)))
428 tensorNetwork->AddGate(xGate,
static_cast<unsigned int>(qubit));
429 }
else if (simulationType == SimulationType::kPauliPropagator) {
430 std::vector<int> qubitsInt(qubits.begin(), qubits.end());
431 const auto res = pp->Measure(qubitsInt);
432 for (
size_t i = 0; i < res.size(); ++i) {
433 if (res[i]) pp->ApplyX(qubitsInt[i]);
436 for (
size_t qubit : qubits)
437 if (mpsSimulator->MeasureQubit(
static_cast<unsigned int>(qubit)))
438 mpsSimulator->ApplyGate(xGate,
static_cast<unsigned int>(qubit));
442 NotifyObservers(qubits);
456 double Probability(Types::qubit_t outcome)
override {
457 if (simulationType == SimulationType::kMatrixProductState)
458 return mpsSimulator->getBasisStateProbability(
459 static_cast<unsigned int>(outcome));
460 else if (simulationType == SimulationType::kStabilizer)
461 return cliffordSimulator->getBasisStateProbability(
462 static_cast<unsigned int>(outcome));
463 else if (simulationType == SimulationType::kTensorNetwork)
464 return tensorNetwork->getBasisStateProbability(outcome);
465 else if (simulationType == SimulationType::kPauliPropagator)
466 return pp->Probability(outcome);
468 return state->getBasisStateProbability(
static_cast<unsigned int>(outcome));
481 std::complex<double>
Amplitude(Types::qubit_t outcome)
override {
482 if (simulationType == SimulationType::kMatrixProductState)
483 return mpsSimulator->getBasisStateAmplitude(
484 static_cast<unsigned int>(outcome));
485 else if (simulationType == SimulationType::kStabilizer)
486 throw std::runtime_error(
487 "QCSimState::Amplitude: Invalid simulation type for obtaining the "
488 "amplitude of the specified outcome.");
489 else if (simulationType == SimulationType::kTensorNetwork)
490 throw std::runtime_error(
491 "QCSimState::Amplitude: Not supported for the "
492 "tensor network simulator.");
493 else if (simulationType == SimulationType::kPauliPropagator)
494 throw std::runtime_error(
495 "QCSimState::Amplitude: Invalid simulation type for obtaining the "
496 "amplitude of the specified outcome.");
498 return state->getBasisStateAmplitude(
static_cast<unsigned int>(outcome));
513 if (simulationType == SimulationType::kTensorNetwork)
514 throw std::runtime_error(
515 "QCSimState::AllProbabilities: Invalid "
516 "simulation type for obtaining probabilities.");
517 else if (simulationType == SimulationType::kStabilizer)
518 return cliffordSimulator->AllProbabilities();
519 else if (simulationType == SimulationType::kPauliPropagator) {
521 std::vector<double> result(nrBasisStates);
522 for (
size_t i = 0; i < nrBasisStates; ++i)
523 result[i] = pp->Probability(i);
527 const Eigen::VectorXcd probs =
528 simulationType == SimulationType::kMatrixProductState
529 ? mpsSimulator->getRegisterStorage().cwiseAbs2()
530 : state->getRegisterStorage().cwiseAbs2();
532 std::vector<double> result(probs.size());
534 for (
int i = 0; i < probs.size(); ++i) result[i] = probs[i].real();
551 const Types::qubits_vector &qubits)
override {
552 if (simulationType == SimulationType::kStabilizer)
553 throw std::runtime_error(
554 "QCSimState::Probabilities: Invalid simulation "
555 "type for obtaining probabilities.");
556 else if (simulationType == SimulationType::kTensorNetwork) {
558 throw std::runtime_error(
559 "QCSimState::Probabilities: Not implemented yet "
560 "for the tensor network simulator.");
563 std::vector<double> result(qubits.size());
565 if (simulationType == SimulationType::kMatrixProductState) {
566 for (
int i = 0; i < static_cast<int>(qubits.size()); ++i)
567 result[i] = mpsSimulator->getBasisStateProbability(qubits[i]);
568 }
else if (simulationType == SimulationType::kPauliPropagator) {
569 for (
int i = 0; i < static_cast<int>(qubits.size()); ++i)
570 result[i] = pp->Probability(qubits[i]);
572 const Eigen::VectorXcd ® = state->getRegisterStorage();
574 for (
int i = 0; i < static_cast<int>(qubits.size()); ++i)
575 result[i] = std::norm(reg[qubits[i]]);
594 std::unordered_map<Types::qubit_t, Types::qubit_t>
SampleCounts(
595 const Types::qubits_vector &qubits,
size_t shots = 1000)
override {
596 if (qubits.empty() || shots == 0)
return {};
601 std::unordered_map<Types::qubit_t, Types::qubit_t> result;
605 if (simulationType == SimulationType::kMatrixProductState) {
607 if (useMPSMeasureNoCollapse) {
609 std::unordered_set qset(qubits.begin(), qubits.end());
613 for (
size_t shot = 0; shot < shots; ++shot) {
619 for (
auto q : qubits) {
620 const size_t qubitMask = 1ULL << q;
621 if (measRaw & qubitMask) meas |= mask;
631 auto savedState = mpsSimulator->getState();
632 for (
size_t shot = 0; shot < shots; ++shot) {
633 const size_t meas =
Measure(qubits);
635 mpsSimulator->setState(savedState);
638 }
else if (simulationType == SimulationType::kStabilizer) {
639 cliffordSimulator->SaveState();
640 for (
size_t shot = 0; shot < shots; ++shot) {
641 const size_t meas =
Measure(qubits);
643 cliffordSimulator->RestoreState();
645 cliffordSimulator->ClearSavedState();
646 }
else if (simulationType == SimulationType::kTensorNetwork) {
647 tensorNetwork->SaveState();
648 for (
size_t shot = 0; shot < shots; ++shot) {
649 const size_t meas =
Measure(qubits);
651 tensorNetwork->RestoreState();
653 tensorNetwork->ClearSavedState();
654 }
else if (simulationType == SimulationType::kPauliPropagator) {
655 std::vector<int> qubitsInt(qubits.begin(), qubits.end());
656 for (
size_t shot = 0; shot < shots; ++shot) {
657 const auto res = pp->Sample(qubitsInt);
660 for (
size_t i = 0; i < qubits.size(); ++i) {
661 if (res[i]) meas |= (1ULL << i);
668 const auto &statev = state->getRegisterStorage();
672 for (
size_t shot = 0; shot < shots; ++shot) {
673 const double prob = 1. - uniformZeroOne(rng);
674 const size_t measRaw = alias.
Sample(prob);
678 for (
auto q : qubits) {
679 const size_t qubitMask = 1ULL << q;
680 if ((measRaw & qubitMask) != 0) meas |= mask;
687 for (
size_t shot = 0; shot < shots; ++shot) {
692 for (
auto q : qubits) {
693 const size_t qubitMask = 1ULL << q;
694 if ((measRaw & qubitMask) != 0) meas |= mask;
704 NotifyObservers(qubits);
720 double ExpectationValue(
const std::string &pauliStringOrig)
override {
721 if (pauliStringOrig.empty())
return 1.0;
723 std::string pauliString = pauliStringOrig;
726 const auto pauliOp = toupper(pauliString[i]);
727 if (pauliOp !=
'I' && pauliOp !=
'Z')
return 0.0;
733 if (simulationType == SimulationType::kStabilizer)
734 return cliffordSimulator->ExpectationValue(pauliString);
735 else if (simulationType == SimulationType::kTensorNetwork)
736 return tensorNetwork->ExpectationValue(pauliString);
737 else if (simulationType == SimulationType::kPauliPropagator)
738 return pp->ExpectationValue(pauliString);
741 static const QC::Gates::PauliXGate<> xgate;
742 static const QC::Gates::PauliYGate<> ygate;
743 static const QC::Gates::PauliZGate<> zgate;
745 std::vector<QC::Gates::AppliedGate<Eigen::MatrixXcd>> pauliStringVec;
746 pauliStringVec.reserve(pauliString.size());
748 for (
size_t q = 0; q < pauliString.size(); ++q) {
749 switch (toupper(pauliString[q])) {
751 QC::Gates::AppliedGate<Eigen::MatrixXcd> ag(
752 xgate.getRawOperatorMatrix(),
static_cast<Types::qubit_t
>(q));
753 pauliStringVec.emplace_back(std::move(ag));
756 QC::Gates::AppliedGate<Eigen::MatrixXcd> ag(
757 ygate.getRawOperatorMatrix(),
static_cast<Types::qubit_t
>(q));
758 pauliStringVec.emplace_back(std::move(ag));
761 QC::Gates::AppliedGate<Eigen::MatrixXcd> ag(
762 zgate.getRawOperatorMatrix(),
static_cast<Types::qubit_t
>(q));
763 pauliStringVec.emplace_back(std::move(ag));
772 if (pauliStringVec.empty())
return 1.0;
774 if (simulationType == SimulationType::kMatrixProductState)
775 return mpsSimulator->ExpectationValue(pauliStringVec).real();
777 return state->ExpectationValue(pauliStringVec).real();
787 SimulatorType GetType()
const override {
return SimulatorType::kQCSim; }
807 void Flush()
override {}
838 if (simulationType == SimulationType::kMatrixProductState)
839 mpsSimulator->SaveState();
840 else if (simulationType == SimulationType::kStabilizer)
841 cliffordSimulator->SaveState();
842 else if (simulationType == SimulationType::kTensorNetwork)
843 tensorNetwork->SaveState();
844 else if (simulationType == SimulationType::kPauliPropagator)
859 if (simulationType == SimulationType::kMatrixProductState)
860 mpsSimulator->RestoreState();
861 else if (simulationType == SimulationType::kStabilizer)
862 cliffordSimulator->RestoreState();
863 else if (simulationType == SimulationType::kTensorNetwork)
864 tensorNetwork->RestoreState();
865 else if (simulationType == SimulationType::kPauliPropagator)
868 state->RestoreState();
878 std::complex<double> AmplitudeRaw(Types::qubit_t outcome)
override {
891 enableMultithreading = multithreading;
892 if (state) state->SetMultithreading(multithreading);
893 if (cliffordSimulator) cliffordSimulator->SetMultithreading(multithreading);
894 if (tensorNetwork) tensorNetwork->SetMultithreading(multithreading);
897 if (multithreading) pp->EnableParallel();
898 else pp->DisableParallel();
921 bool IsQcsim()
const override {
return true; }
937 if (simulationType == SimulationType::kStatevector)
938 return state->MeasureNoCollapse();
939 else if (simulationType == SimulationType::kMatrixProductState) {
940 const auto measured = mpsSimulator->MeasureNoCollapse();
941 Types::qubit_t result = 0;
942 Types::qubit_t mask = 1;
943 for (
const auto &meas : measured) {
944 if (meas) result |= mask;
948 }
else if (simulationType == SimulationType::kPauliPropagator) {
950 std::iota(qubitsInt.begin(), qubitsInt.end(), 0);
951 const auto res = pp->Sample(qubitsInt);
952 Types::qubit_t result = 0;
953 for (
size_t i = 0; i < res.size(); ++i) {
954 if (res[i]) result |= (1ULL << i);
959 throw std::runtime_error(
960 "QCSimState::MeasureNoCollapse: Invalid simulation type for measuring "
961 "all the qubits without collapsing the state.");
967 SimulationType simulationType =
968 SimulationType::kStatevector;
970 std::unique_ptr<QC::QubitRegister<>> state;
971 std::unique_ptr<QC::TensorNetworks::MPSSimulator>
973 std::unique_ptr<QC::Clifford::StabilizerSimulator>
975 std::unique_ptr<TensorNetworks::TensorNetwork>
977 std::unique_ptr<QcsimPauliPropagator>
981 bool limitSize =
false;
982 bool limitEntanglement =
false;
983 Eigen::Index chi = 10;
984 double singularValueThreshold = 0.;
985 bool enableMultithreading =
true;
986 bool useMPSMeasureNoCollapse =
990 std::uniform_real_distribution<double> uniformZeroOne;
double Probability(void *sim, unsigned long long int outcome)
char * GetConfiguration(void *sim, const char *key)
int RestoreState(void *sim)
int ApplyReset(void *sim, const unsigned long int *qubits, unsigned long int nrQubits)
unsigned long int AllocateQubits(void *sim, unsigned long int nrQubits)
unsigned long int GetNumberOfQubits(void *sim)
double * AllProbabilities(void *sim)
unsigned long long int MeasureNoCollapse(void *sim)
int GetMultithreading(void *sim)
unsigned long long int Measure(void *sim, const unsigned long int *qubits, unsigned long int nrQubits)
double * Amplitude(void *sim, unsigned long long int outcome)
double * Probabilities(void *sim, const unsigned long long int *qubits, unsigned long int nrQubits)
int SetMultithreading(void *sim, int multithreading)
int SaveStateToInternalDestructive(void *sim)
int GetSimulationType(void *sim)
unsigned long long int * SampleCounts(void *sim, const unsigned long long int *qubits, unsigned long int nrQubits, unsigned long int shots)
int RestoreInternalDestructiveSavedState(void *sim)
size_t Sample(double v) const