From a1557ba2aa35e5fad885147a0f0423cef1447c73 Mon Sep 17 00:00:00 2001 From: Thien Nguyen Date: Fri, 22 Jan 2021 12:11:41 -0500 Subject: [PATCH] Aer able to compute exp-val-z when simulating density matrix The procedure is similar to state-vector but using the diagonal elements of the density matrix rather than computing the square norm of the state vector. Signed-off-by: Thien Nguyen --- .../ibm/aer/accelerator/aer_accelerator.cpp | 32 ++++++- .../ibm/aer/accelerator/aer_accelerator.hpp | 4 + .../tests/JsonNoiseModelTester.cpp | 88 +++++++++++++++++++ 3 files changed, 122 insertions(+), 2 deletions(-) diff --git a/quantum/plugins/ibm/aer/accelerator/aer_accelerator.cpp b/quantum/plugins/ibm/aer/accelerator/aer_accelerator.cpp index 331674fc6..fe9fb8d55 100644 --- a/quantum/plugins/ibm/aer/accelerator/aer_accelerator.cpp +++ b/quantum/plugins/ibm/aer/accelerator/aer_accelerator.cpp @@ -141,6 +141,32 @@ double AerAccelerator::calcExpectationValueZ( return result; } +double AerAccelerator::calcExpectationValueZFromDensityMatrix( + const std::vector>> &in_densityMat, + const std::vector &in_bits) { + const auto hasEvenParity = + [](size_t x, const std::vector &in_qubitIndices) -> bool { + size_t count = 0; + for (const auto &bitIdx : in_qubitIndices) { + if (x & (1ULL << bitIdx)) { + count++; + } + } + return (count % 2) == 0; + }; + + double result = 0.0; + for (uint64_t i = 0; i < in_densityMat.size(); ++i) { + const auto &diag_elem = in_densityMat[i][i]; + // The diag. elements of the DM should be real. + assert(std::abs(diag_elem.second) < 1e-3); + // When using DM elements, don't need the square-norm. + result += ((hasEvenParity(i, in_bits) ? 1.0 : -1.0) * diag_elem.first); + } + + return result; +} + void AerAccelerator::execute( std::shared_ptr buffer, const std::shared_ptr program) { @@ -271,10 +297,10 @@ void AerAccelerator::execute( snapshotInst["snapshot_type"] = "density_matrix"; auto& exprJson = *(j["experiments"].begin()); exprJson["instructions"].push_back(snapshotInst); - std::cout << "Qobj:\n" << j.dump(2); + // std::cout << "Qobj:\n" << j.dump(2); auto results_json = nlohmann::json::parse( AER::controller_execute_json(j.dump())); - std::cout << "Result:\n" << results_json.dump() << "\n"; + // std::cout << "Result:\n" << results_json.dump() << "\n"; auto results = *results_json["results"].begin(); auto dm_mat = (*(results["data"]["snapshots"]["density_matrix"]["dm_snapshot"] @@ -285,6 +311,8 @@ void AerAccelerator::execute( flattenDm.insert(flattenDm.end(), row.begin(), row.end()); } buffer->addExtraInfo("density_matrix", flattenDm); + auto exp_val = calcExpectationValueZFromDensityMatrix(dm_mat, measured_bits); + buffer->addExtraInfo("exp-val-z", exp_val); } else { // statevector diff --git a/quantum/plugins/ibm/aer/accelerator/aer_accelerator.hpp b/quantum/plugins/ibm/aer/accelerator/aer_accelerator.hpp index a51d9593b..b8cdd6751 100644 --- a/quantum/plugins/ibm/aer/accelerator/aer_accelerator.hpp +++ b/quantum/plugins/ibm/aer/accelerator/aer_accelerator.hpp @@ -66,6 +66,10 @@ class AerAccelerator : public Accelerator { const std::vector> &in_stateVec, const std::vector &in_bits); + static double calcExpectationValueZFromDensityMatrix( + const std::vector>> &in_densityMat, + const std::vector &in_bits); + std::shared_ptr xacc_to_qobj; int m_shots = 1024; std::string m_simtype = "qasm"; diff --git a/quantum/plugins/noise_model/tests/JsonNoiseModelTester.cpp b/quantum/plugins/noise_model/tests/JsonNoiseModelTester.cpp index 732f06332..0aa179bce 100644 --- a/quantum/plugins/noise_model/tests/JsonNoiseModelTester.cpp +++ b/quantum/plugins/noise_model/tests/JsonNoiseModelTester.cpp @@ -2,6 +2,8 @@ #include "xacc.hpp" #include "xacc_service.hpp" #include "NoiseModel.hpp" +#include "PauliOperator.hpp" +#include "xacc_observable.hpp" namespace { // A sample Json for testing @@ -260,6 +262,92 @@ TEST(JsonNoiseModelTester, checkIbmNoiseJsonIndexing) { EXPECT_EQ(densityMatrix.size(), 16); } +// Test that when doing noisy simulation w/ Aer using density matrix, +// the result is consistent (no sampling involved) +TEST(JsonNoiseModelTester, testVqeModeWithDensityMatrixSim) { + auto noiseModel = xacc::getService("json"); + noiseModel->initialize({{"noise-model", depol_json_2q}}); + const std::string ibmNoiseJson = noiseModel->toJson(); + std::cout << "IBM Equiv: \n" << ibmNoiseJson << "\n"; + auto accelerator = xacc::getAccelerator( + "aer", {{"noise-model", ibmNoiseJson}, {"sim-type", "density_matrix"}}); + auto xasmCompiler = xacc::getCompiler("xasm"); + + auto ir = + xasmCompiler->compile(R"(__qpu__ void ansatz_temp(qbit q, double t) { + X(q[0]); + Ry(q[1], t); + CX(q[1], q[0]); + H(q[0]); + H(q[1]); + Measure(q[0]); + Measure(q[1]); + })", + accelerator); + + auto program = ir->getComposite("ansatz_temp"); + // Data for first run: + std::vector firstRuns; + const auto angles = + xacc::linspace(-xacc::constants::pi, xacc::constants::pi, 20); + for (size_t i = 0; i < angles.size(); ++i) { + auto buffer = xacc::qalloc(2); + auto evaled = program->operator()({angles[i]}); + accelerator->execute(buffer, evaled); + firstRuns.emplace_back(buffer->getExpectationValueZ()); + } + + // Run a second-time, the exp-val must be consistent... + for (size_t i = 0; i < angles.size(); ++i) { + auto buffer = xacc::qalloc(2); + auto evaled = program->operator()({angles[i]}); + accelerator->execute(buffer, evaled); + EXPECT_NEAR(buffer->getExpectationValueZ(), firstRuns[i], 1e-9); + } +} + +// Make sure the exp-val calc. is correct by running VQE opt. +TEST(JsonNoiseModelTester, testVqe) { + auto noiseModel = xacc::getService("json"); + noiseModel->initialize({{"noise-model", depol_json_2q}}); + const std::string ibmNoiseJson = noiseModel->toJson(); + std::cout << "IBM Equiv: \n" << ibmNoiseJson << "\n"; + auto accelerator = xacc::getAccelerator( + "aer", {{"noise-model", ibmNoiseJson}, {"sim-type", "density_matrix"}}); + + // Create the N=2 deuteron Hamiltonian + auto H_N_2 = xacc::quantum::getObservable( + "pauli", std::string("5.907 - 2.1433 X0X1 " + "- 2.1433 Y0Y1" + "+ .21829 Z0 - 6.125 Z1")); + + auto optimizer = xacc::getOptimizer("nlopt"); + xacc::qasm(R"( + .compiler xasm + .circuit deuteron_ansatz + .parameters theta + .qbit q + X(q[0]); + Ry(q[1], theta); + CNOT(q[1],q[0]); + )"); + auto ansatz = xacc::getCompiled("deuteron_ansatz"); + + // Get the VQE Algorithm and initialize it + auto vqe = xacc::getAlgorithm("vqe"); + vqe->initialize({std::make_pair("ansatz", ansatz), + std::make_pair("observable", H_N_2), + std::make_pair("accelerator", accelerator), + std::make_pair("optimizer", optimizer)}); + + // Allocate some qubits and execute + auto buffer = xacc::qalloc(2); + vqe->execute(buffer); + + // Expected result: -1.74886 + EXPECT_NEAR((*buffer)["opt-val"].as(), -1.74886, 0.1); +} + int main(int argc, char **argv) { xacc::Initialize(argc, argv); ::testing::InitGoogleTest(&argc, argv);