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 <[email protected]>

quantum/plugins/ibm/aer/accelerator/aer_accelerator.cpp

@@ -141,6 +141,32 @@ double AerAccelerator::calcExpectationValueZ(
   return result;
 }
 
+double AerAccelerator::calcExpectationValueZFromDensityMatrix(
+    const std::vector<std::vector<std::pair<double, double>>> &in_densityMat,
+    const std::vector<std::size_t> &in_bits) {
+  const auto hasEvenParity =
+      [](size_t x, const std::vector<size_t> &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<AcceleratorBuffer> buffer,
     const std::shared_ptr<CompositeInstruction> 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<AER::Simulator::QasmController>(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

quantum/plugins/ibm/aer/accelerator/aer_accelerator.hpp

@@ -66,6 +66,10 @@ class AerAccelerator : public Accelerator {
       const std::vector<std::pair<double, double>> &in_stateVec,
       const std::vector<std::size_t> &in_bits);
 
+  static double calcExpectationValueZFromDensityMatrix(
+      const std::vector<std::vector<std::pair<double, double>>> &in_densityMat,
+      const std::vector<std::size_t> &in_bits);
+
   std::shared_ptr<Compiler> xacc_to_qobj;
   int m_shots = 1024;
   std::string m_simtype = "qasm";

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<xacc::NoiseModel>("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<double> 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<xacc::NoiseModel>("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<double>(), -1.74886, 0.1);
+}
+
 int main(int argc, char **argv) {
   xacc::Initialize(argc, argv);
   ::testing::InitGoogleTest(&argc, argv);