Fixed https://github.com/eclipse/xacc/issues/509

- Gradient strategies to filter out non-observing circuits associated with identity term.

- For two-sided gradients (central and parameter shift), the identity term has no effect. For single-sided ones (forward and backward), add the identity coefficient in the shifted enery value.

- Fixed an oversight in setting optional 'step' param. It was only set in the ctor, but not reset during initialization. Hence, it may carry the value set by a previous run.

- Added a set of unit test cases covering 'shots' mode
Signed-off-by: Thien Nguyen <thien.md.nguyen@gmail.com>

quantum/plugins/algorithms/gradient_strategies/BackwardDifferenceGradient.hpp

@@ -18,6 +18,7 @@
 #include "xacc_service.hpp"
 #include "AlgorithmGradientStrategy.hpp"
 #include <iomanip>
+#include "InstructionIterator.hpp"
 
 using namespace xacc;
 
@@ -52,6 +53,8 @@ public:
 
     obs = parameters.get<std::shared_ptr<Observable>>("observable");
 
+    // Default step size
+    step = 1.0e-7; 
     // Change step size if need be
     if (parameters.keyExists<double>("step")) {
       step = parameters.get<double>("step");
@@ -70,7 +73,18 @@ public:
   std::vector<std::shared_ptr<CompositeInstruction>>
   getGradientExecutions(std::shared_ptr<CompositeInstruction> circuit,
                         const std::vector<double> &x) override {
-
+    // Check if the composite contains measure gates
+    const auto containMeasureGates =
+        [](std::shared_ptr<CompositeInstruction> f) -> bool {
+      InstructionIterator it(f);
+      while (it.hasNext()) {
+        auto nextInst = it.next();
+        if (nextInst->name() == "Measure") {
+          return true;
+        }
+      }
+      return false;
+    };
     // std::stringstream ss;
     // ss << std::setprecision(5) << "Input parameters: ";
     // for (auto param : x) {
@@ -90,8 +104,10 @@ public:
         auto evaled_base = kernel_evaluator(tmpX);
         kernels = obs->observe(evaled_base);
         for (auto &f : kernels) {
-          coefficients.push_back(std::real(f->getCoefficient()));
-          gradientInstructions.push_back(f);
+          if (containMeasureGates(f)) {
+            coefficients.push_back(std::real(f->getCoefficient()));
+            gradientInstructions.push_back(f);
+          }
         }
       } else {
         kernels = obs->observe(circuit);
@@ -99,13 +115,18 @@ public:
         // loop over circuit instructions
         // and gather coefficients/instructions
         for (auto &f : kernels) {
-          auto evaled = f->operator()(tmpX);
-          coefficients.push_back(std::real(f->getCoefficient()));
-          gradientInstructions.push_back(evaled);
+          if (containMeasureGates(f)) {
+            auto evaled = f->operator()(tmpX);
+            coefficients.push_back(std::real(f->getCoefficient()));
+            gradientInstructions.push_back(evaled);
+          }
         }
       }
 
-      nInstructionsElement.push_back(kernels.size());
+      const auto numberGradKernels = std::count_if(
+          kernels.begin(), kernels.end(),
+          [&containMeasureGates](auto &f) { return containMeasureGates(f); });
+      nInstructionsElement.push_back(numberGradKernels);
     }
 
     return gradientInstructions;
@@ -119,8 +140,9 @@ public:
     int shift = 0;
     // loop over the remaining number of entries in the gradient vector
     for (int gradTerm = 0; gradTerm < dx.size(); gradTerm++) {
-
-      double gradElement = 0.0;
+      auto identityTerm = obs->getIdentitySubTerm();
+      double gradElement =
+          identityTerm ? identityTerm->coefficient().real() : 0.0;
 
       // loop over instructions for a given term, compute <+> and <->
       for (int instElement = 0; instElement < nInstructionsElement[gradTerm];

quantum/plugins/algorithms/gradient_strategies/CentralDifferenceGradient.hpp

@@ -18,6 +18,7 @@
 #include "xacc_service.hpp"
 #include "AlgorithmGradientStrategy.hpp"
 #include <iomanip>
+#include "InstructionIterator.hpp"
 
 using namespace xacc;
 
@@ -52,6 +53,8 @@ public:
 
     obs = parameters.get<std::shared_ptr<Observable>>("observable");
 
+    // Default value
+    step = 1.0e-7; 
     // Change step size if need be
     if (parameters.keyExists<double>("step")) {
       step = parameters.get<double>("step");
@@ -72,7 +75,18 @@ public:
   std::vector<std::shared_ptr<CompositeInstruction>>
   getGradientExecutions(std::shared_ptr<CompositeInstruction> circuit,
                         const std::vector<double> &x) override {
-
+    // Check if the composite contains measure gates
+    const auto containMeasureGates =
+        [](std::shared_ptr<CompositeInstruction> f) -> bool {
+      InstructionIterator it(f);
+      while (it.hasNext()) {
+        auto nextInst = it.next();
+        if (nextInst->name() == "Measure") {
+          return true;
+        }
+      }
+      return false;
+    };
     std::vector<std::shared_ptr<CompositeInstruction>> gradientInstructions;
     for (int op = 0; op < x.size(); op++) { // loop over operators
       for (double sign : {1.0, -1.0}) {     // change sign
@@ -86,8 +100,10 @@ public:
           auto evaled_base = kernel_evaluator(tmpX);
           kernels = obs->observe(evaled_base);
           for (auto &f : kernels) {
-            coefficients.push_back(std::real(f->getCoefficient()));
-            gradientInstructions.push_back(f);
+            if (containMeasureGates(f)) {
+              coefficients.push_back(std::real(f->getCoefficient()));
+              gradientInstructions.push_back(f);
+            }
           }
         } else {
           kernels = obs->observe(circuit);
@@ -95,15 +111,21 @@ public:
           // loop over circuit instructions
           // and gather coefficients/instructions
           for (auto &f : kernels) {
-            auto evaled = f->operator()(tmpX);
-            coefficients.push_back(std::real(f->getCoefficient()));
-            gradientInstructions.push_back(evaled);
+            if (containMeasureGates(f)) {
+              auto evaled = f->operator()(tmpX);
+              coefficients.push_back(std::real(f->getCoefficient()));
+              gradientInstructions.push_back(evaled);
+            }
           }
         }
         // the number of instructions for a given element of x is the same
         // regardless of the parameter sign, so we need only one of this
         if (sign == 1.0) {
-          nInstructionsElement.push_back(kernels.size());
+          const auto numberGradKernels = std::count_if(
+              kernels.begin(), kernels.end(), [&containMeasureGates](auto &f) {
+                return containMeasureGates(f);
+              });
+          nInstructionsElement.push_back(numberGradKernels);
         }
       }
     }

quantum/plugins/algorithms/gradient_strategies/ForwardDifferenceGradient.hpp

@@ -18,6 +18,7 @@
 #include "xacc_service.hpp"
 #include "AlgorithmGradientStrategy.hpp"
 #include <iomanip>
+#include "InstructionIterator.hpp"
 
 using namespace xacc;
 
@@ -52,6 +53,8 @@ public:
 
     obs = parameters.get<std::shared_ptr<Observable>>("observable");
 
+    // Default step size
+    step = 1.0e-7; 
     // Change step size if need be
     if (parameters.keyExists<double>("step")) {
       step = parameters.get<double>("step");
@@ -71,6 +74,18 @@ public:
   getGradientExecutions(std::shared_ptr<CompositeInstruction> circuit,
                         const std::vector<double> &x) override {
 
+    // Check if the composite contains measure gates
+    const auto containMeasureGates =
+        [](std::shared_ptr<CompositeInstruction> f) -> bool {
+      InstructionIterator it(f);
+      while (it.hasNext()) {
+        auto nextInst = it.next();
+        if (nextInst->name() == "Measure") {
+          return true;
+        }
+      }
+      return false;
+    };
     // std::stringstream ss;
     // ss << std::setprecision(5) << "Input parameters: ";
     // for (auto param : x) {
@@ -92,8 +107,10 @@ public:
         auto evaled_base = kernel_evaluator(tmpX);
         kernels = obs->observe(evaled_base);
         for (auto &f : kernels) {
-          coefficients.push_back(std::real(f->getCoefficient()));
-          gradientInstructions.push_back(f);
+          if (containMeasureGates(f)) {
+            coefficients.push_back(std::real(f->getCoefficient()));
+            gradientInstructions.push_back(f);
+          }
         }
       } else {
         kernels = obs->observe(circuit);
@@ -101,13 +118,18 @@ public:
         // loop over circuit instructions
         // and gather coefficients/instructions
         for (auto &f : kernels) {
-          auto evaled = f->operator()(tmpX);
-          coefficients.push_back(std::real(f->getCoefficient()));
-          gradientInstructions.push_back(evaled);
+          if (containMeasureGates(f)) {
+            auto evaled = f->operator()(tmpX);
+            coefficients.push_back(std::real(f->getCoefficient()));
+            gradientInstructions.push_back(evaled);
+          }
         }
       }
 
-      nInstructionsElement.push_back(kernels.size());
+      const auto numberGradKernels = std::count_if(
+          kernels.begin(), kernels.end(),
+          [&containMeasureGates](auto &f) { return containMeasureGates(f); });
+      nInstructionsElement.push_back(numberGradKernels);
     }
 
     return gradientInstructions;
@@ -117,12 +139,13 @@ public:
   void
   compute(std::vector<double> &dx,
           std::vector<std::shared_ptr<AcceleratorBuffer>> results) override {
-
     int shift = 0;
     // loop over the remaining number of entries in the gradient vector
     for (int gradTerm = 0; gradTerm < dx.size(); gradTerm++) {
 
-      double gradElement = 0.0;
+      auto identityTerm = obs->getIdentitySubTerm();
+      double gradElement =
+          identityTerm ? identityTerm->coefficient().real() : 0.0;
 
       // loop over instructions for a given term, compute <+> and <->
       for (int instElement = 0; instElement < nInstructionsElement[gradTerm];
@@ -151,10 +174,7 @@ public:
     return;
   }
 
-
-  const std::string name() const override {
-    return "forward";
-  }
+  const std::string name() const override { return "forward"; }
   const std::string description() const override { return ""; }
 };
 

quantum/plugins/algorithms/gradient_strategies/ParameterShiftGradient.hpp

@@ -14,6 +14,7 @@
 #define XACC_PARAMETER_SHIFT_GRADIENT_HPP_
 
 #include "CompositeInstruction.hpp"
+#include "InstructionIterator.hpp"
 #include "xacc.hpp"
 #include "xacc_service.hpp"
 #include "AlgorithmGradientStrategy.hpp"
@@ -48,7 +49,9 @@ public:
           parameters.get<std::function<std::shared_ptr<CompositeInstruction>(
               std::vector<double>)>>("kernel-evaluator");
     }
-
+    // Default shiftScalar (this is not clonable, hence need to be
+    // reinitialized)
+    shiftScalar = 0.25;
     if (parameters.keyExists<double>("shift-scalar")) {
       shiftScalar = parameters.get<double>("shift-scalar");
     }
@@ -63,9 +66,23 @@ public:
   std::vector<std::shared_ptr<CompositeInstruction>>
   getGradientExecutions(std::shared_ptr<CompositeInstruction> circuit,
                         const std::vector<double> &x) override {
+    // Check if the composite contains measure gates
+    const auto containMeasureGates =
+        [](std::shared_ptr<CompositeInstruction> f) -> bool {
+      InstructionIterator it(f);
+      while (it.hasNext()) {
+        auto nextInst = it.next();
+        if (nextInst->name() == "Measure") {
+          return true;
+        }
+      }
+      return false;
+    };
 
     std::vector<std::shared_ptr<CompositeInstruction>> gradientInstructions;
-
+    // Note: for the purpose of parameter-shift gradient calculation, the
+    // identity term has no effect, i.e., its contributions to the plus and
+    // minus sides will cancel out.
     // loop over parameters
     for (int param = 0; param < x.size(); param++) {
       // parameter shift sign
@@ -81,8 +98,10 @@ public:
           auto evaled_base = kernel_evaluator(tmpX);
           kernels = obs->observe(evaled_base);
           for (auto &f : kernels) {
-            coefficients.push_back(std::real(f->getCoefficient()));
-            gradientInstructions.push_back(f);
+            if (containMeasureGates(f)) {
+              coefficients.push_back(std::real(f->getCoefficient()));
+              gradientInstructions.push_back(f);
+            }
           }
         } else {
           kernels = obs->observe(circuit);
@@ -90,16 +109,22 @@ public:
           // loop over circuit instructions
           // and gather coefficients/instructions
           for (auto &f : kernels) {
-            auto evaled = f->operator()(tmpX);
-            coefficients.push_back(std::real(f->getCoefficient()));
-            gradientInstructions.push_back(evaled);
+            if (containMeasureGates(f)) {
+              auto evaled = f->operator()(tmpX);
+              coefficients.push_back(std::real(f->getCoefficient()));
+              gradientInstructions.push_back(evaled);
+            }
           }
         }
 
         // the number of instructions for a given element of x is the same
-        // regardless of the parameter sign, so we need only one of this
+        // regardless of the parameter sign, so we need only one of this.
         if (sign == 1.0) {
-          nInstructionsElement.push_back(kernels.size());
+          const auto numberGradKernels = std::count_if(
+              kernels.begin(), kernels.end(), [&containMeasureGates](auto &f) {
+                return containMeasureGates(f);
+              });
+          nInstructionsElement.push_back(numberGradKernels);
         }
       }
     }
@@ -121,7 +146,6 @@ public:
       // loop over instructions for a given term, compute <+> and <->
       for (int instElement = 0; instElement < nInstructionsElement[gradTerm];
            instElement++) {
-
         auto plus_expval =
             std::real(results[instElement + shift]->getExpectationValueZ());
         auto minus_expval = std::real(

quantum/plugins/algorithms/gradient_strategies/tests/GradientStrategiesTester.cpp

@@ -200,6 +200,151 @@ TEST(GradientStrategiesTester, checkYanPSproblem) {
   EXPECT_NEAR((*q)["opt-val"].as<double>(), -1.0, 1e-4);
 }
 
+// Check running shots:
+TEST(GradientStrategiesTester, checkParameterShiftShots) {
+  auto accelerator = xacc::getAccelerator("aer", {{"shots", 65536}});
+  auto buffer = xacc::qalloc(1);
+
+  std::shared_ptr<Observable> observable =
+      std::make_shared<xacc::quantum::PauliOperator>();
+  observable->fromString("X0");
+
+  auto provider = xacc::getIRProvider("quantum");
+  auto ansatz = provider->createComposite("testCircuit");
+  ansatz->addVariable("x0");
+  ansatz->addInstruction(provider->createInstruction(
+      "Ry", std::vector<std::size_t>{0}, {InstructionParameter("x0")}));
+
+  auto parameterShift =
+      xacc::getService<AlgorithmGradientStrategy>("parameter-shift");
+  parameterShift->initialize({std::make_pair("observable", observable)});
+  auto gradientInstructions =
+      parameterShift->getGradientExecutions(ansatz, {0.0});
+  accelerator->execute(buffer, gradientInstructions);
+  std::vector<double> dx(1);
+  parameterShift->compute(dx, buffer->getChildren());
+  EXPECT_NEAR(dx[0], std::sqrt(2), 0.1);
+}
+
+TEST(GradientStrategiesTester, checkCentralDifferenceShots) {
+  auto accelerator2 = xacc::getAccelerator("aer", {{"shots", 65536}});
+  // auto accelerator2 = xacc::getAccelerator("qpp");
+  auto buffer2 = xacc::qalloc(1);
+
+  std::shared_ptr<Observable> observable2 =
+      std::make_shared<xacc::quantum::PauliOperator>();
+  observable2->fromString("X0");
+
+  auto provider2 = xacc::getIRProvider("quantum");
+  auto ansatz2 = provider2->createComposite("testCircuit");
+  ansatz2->addVariable("x0");
+  ansatz2->addInstruction(provider2->createInstruction(
+      "Ry", std::vector<std::size_t>{0}, {InstructionParameter("x0")}));
+
+  auto centralDifference =
+      xacc::getService<AlgorithmGradientStrategy>("central");
+  centralDifference->initialize({{"observable", observable2}, {"step", 0.1}});
+  auto gradientInstructions =
+      centralDifference->getGradientExecutions(ansatz2, {0.0});
+  accelerator2->execute(buffer2, gradientInstructions);
+
+  std::vector<double> dx(1);
+  centralDifference->compute(dx, buffer2->getChildren());
+  EXPECT_NEAR(dx[0], 1.0, 0.1);
+}
+
+TEST(GradientStrategiesTester, checkForwardDifferenceShots) {
+  auto accelerator3 = xacc::getAccelerator("aer", {{"shots", 65536}});
+  auto buffer3 = xacc::qalloc(1);
+
+  std::shared_ptr<Observable> observable3 =
+      std::make_shared<xacc::quantum::PauliOperator>();
+  observable3->fromString("X0");
+
+  auto provider3 = xacc::getIRProvider("quantum");
+  auto ansatz3 = provider3->createComposite("testCircuit");
+  ansatz3->addVariable("x0");
+  ansatz3->addInstruction(provider3->createInstruction(
+      "Ry", std::vector<std::size_t>{0}, {InstructionParameter("x0")}));
+
+  auto forwardDifference =
+      xacc::getService<AlgorithmGradientStrategy>("forward");
+  forwardDifference->initialize({{"observable", observable3}, {"step", 0.1}});
+  auto gradientInstructions =
+      forwardDifference->getGradientExecutions(ansatz3, {0.0});
+  accelerator3->execute(buffer3, gradientInstructions);
+
+  std::vector<double> dx(1);
+  forwardDifference->compute(dx, buffer3->getChildren());
+  EXPECT_NEAR(dx[0], 1.0, 0.1);
+}
+
+TEST(GradientStrategiesTester, checkBackwardDifferenceShots) {
+  auto accelerator4 = xacc::getAccelerator("aer", {{"shots", 65536}});
+  auto buffer4 = xacc::qalloc(1);
+
+  std::shared_ptr<Observable> observable4 =
+      std::make_shared<xacc::quantum::PauliOperator>();
+  observable4->fromString("X0");
+
+  auto provider4 = xacc::getIRProvider("quantum");
+  auto ansatz4 = provider4->createComposite("testCircuit");
+  ansatz4->addVariable("x0");
+  ansatz4->addInstruction(provider4->createInstruction(
+      "Ry", std::vector<std::size_t>{0}, {InstructionParameter("x0")}));
+
+  auto backwardDifference =
+      xacc::getService<AlgorithmGradientStrategy>("backward");
+  backwardDifference->initialize({{"observable", observable4}, {"step", 0.1}});
+  auto gradientInstructions =
+      backwardDifference->getGradientExecutions(ansatz4, {0.0});
+  accelerator4->execute(buffer4, gradientInstructions);
+
+  std::vector<double> dx(1);
+  backwardDifference->compute(dx, buffer4->getChildren());
+  EXPECT_NEAR(dx[0], 1.0, 0.1);
+}
+
+TEST(GradientStrategiesTester, checkDeuteronVQEShots) {
+  // Test for: https://github.com/eclipse/xacc/issues/509
+  auto accelerator = xacc::getAccelerator("aer", {{"shots", 65536}});
+
+  // 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", {{"nlopt-optimizer", "l-bfgs"}, {"nlopt-ftol", 0.1}});
+  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");
+  auto centralDifference =
+      xacc::getService<AlgorithmGradientStrategy>("central");
+  centralDifference->initialize({{"observable", H_N_2}, {"step", 0.1}});
+  // Get the VQE Algorithm and initialize it
+  auto vqe = xacc::getAlgorithm("vqe");
+  vqe->initialize({{"ansatz", ansatz},
+                   {"observable", H_N_2},
+                   {"accelerator", accelerator},
+                   {"optimizer", optimizer},
+                   {"gradient_strategy", centralDifference}});
+
+  // 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.25);
+}
+
 int main(int argc, char **argv) {
   xacc::Initialize(argc, argv);
   ::testing::InitGoogleTest(&argc, argv);