Updated Accelerator and VQE to work with measurement bases

quantum/observable/pauli/PauliOperator.cpp

@@ -212,7 +212,6 @@ PauliOperator::observe(std::shared_ptr<CompositeInstruction> function) {
   auto it1 = basisRotations.begin();
   auto it2 = terms.begin();
   for (; it1 != basisRotations.end() && it2 != terms.end(); ++it1, ++it2) {
-
     Term spinInst = (*it2).second;
 
     auto gateFunction =
@@ -1138,7 +1137,9 @@ std::vector<std::shared_ptr<CompositeInstruction>> PauliOperator::getMeasurement
         terms.push_back({kv.first, kv.second});
       }
     }
+
     auto rotationGates = gateRegistry->createComposite(inst.first);
+    rotationGates->setCoefficient(spinInst.coeff());
     std::vector<std::size_t> measureIdxs;
     for (int i = terms.size() - 1; i >= 0; i--) {
 
@@ -1167,9 +1168,8 @@ std::vector<std::shared_ptr<CompositeInstruction>> PauliOperator::getMeasurement
           meas->setBufferNames({buf_name});
         rotationGates->addInstruction(meas);
       }
-    basisRotations.push_back(rotationGates);
+      basisRotations.push_back(rotationGates);
   }
-
   return basisRotations;
 }
 

quantum/plugins/algorithms/vqe/vqe.cpp

@@ -18,9 +18,7 @@
 #include "xacc_service.hpp"
 #include "AcceleratorDecorator.hpp"
 
-#include <memory>
 #include <iomanip>
-#include <string>
 
 using namespace xacc;
 
@@ -55,8 +53,6 @@ bool VQE::initialize(const HeterogeneousMap &parameters) {
     gradientStrategy =
         parameters.get<std::shared_ptr<AlgorithmGradientStrategy>>(
             "gradient_strategy");
-    // gradientStrategy->initialize({std::make_pair("observable",
-    // xacc::as_shared_ptr(observable))});
   }
 
   if (parameters.stringExists("gradient_strategy") &&
@@ -82,6 +78,11 @@ bool VQE::initialize(const HeterogeneousMap &parameters) {
     gradientStrategy = xacc::getService<AlgorithmGradientStrategy>("autodiff");
     gradientStrategy->initialize(parameters);
   }
+
+  cacheMeasurements = false;
+  if (parameters.keyExists<bool>("cache-measurement-basis")) {
+    cacheMeasurements = parameters.get<bool>("cache-measurement-basis");
+  }
   return true;
 }
 
@@ -98,47 +99,67 @@ void VQE::execute(const std::shared_ptr<AcceleratorBuffer> buffer) const {
 
   std::vector<std::shared_ptr<AcceleratorBuffer>> min_child_buffers;
 
-  // auto kernels = observable->observe(xacc::as_shared_ptr(kernel));
   // Cache of energy values during iterations.
   std::vector<double> energies;
   double last_energy = std::numeric_limits<double>::max();
 
+  if (cacheMeasurements && basisRotations.empty())
+    basisRotations = observable->getMeasurementBasisRotations();
+
   // Here we just need to make a lambda kernel
   // to optimize that makes calls to the targeted QPU.
   OptFunction f(
       [&, this](const std::vector<double> &x, std::vector<double> &dx) {
         std::vector<double> coefficients;
         std::vector<std::string> kernelNames;
         std::vector<std::shared_ptr<CompositeInstruction>> fsToExec;
+        double identityCoeff = 0.0;
+        int nInstructionsEnergy, nInstructionsGradient = 0;
 
         // call CompositeInstruction::operator()()
         auto evaled = kernel->operator()(x);
-        // observe
-        auto kernels = observable->observe(evaled);
 
-        double identityCoeff = 0.0;
-        int nInstructionsEnergy = kernels.size(), nInstructionsGradient = 0;
-        for (auto &f : kernels) {
-          kernelNames.push_back(f->name());
-          std::complex<double> coeff = f->getCoefficient();
-
-          int nFunctionInstructions;
-          if (f->getInstruction(0)->isComposite()) {
-            nFunctionInstructions =
-                kernel->nInstructions() + f->nInstructions() - 1;
-          } else {
-            nFunctionInstructions = f->nInstructions();
-          }
+        // only deal with the measurement basis instead of entire circuits
+        if (cacheMeasurements) {
 
-          if (nFunctionInstructions > kernel->nInstructions()) {
-            fsToExec.push_back(f);
-            coefficients.push_back(std::real(coeff));
-          } else {
-            identityCoeff += std::real(coeff);
+          nInstructionsEnergy = basisRotations.size() - 1;
+          for (auto it = basisRotations.begin(); it != basisRotations.end();) {
+
+            kernelNames.push_back((*it)->name());
+            std::complex<double> coeff = (*it)->getCoefficient();
+            if ((*it)->name() == "I") {
+              identityCoeff += std::real(coeff);
+              it = basisRotations.erase(it);
+            }
             coefficients.push_back(std::real(coeff));
+            ++it;
           }
-        }
 
+        } else {
+
+          // observe
+          auto kernels = observable->observe(evaled);
+          for (auto &f : kernels) {
+            kernelNames.push_back(f->name());
+            std::complex<double> coeff = f->getCoefficient();
+
+            int nFunctionInstructions;
+            if (f->getInstruction(0)->isComposite()) {
+              nFunctionInstructions =
+                  kernel->nInstructions() + f->nInstructions() - 1;
+            } else {
+              nFunctionInstructions = f->nInstructions();
+            }
+
+            if (nFunctionInstructions > kernel->nInstructions()) {
+              fsToExec.push_back(f);
+              coefficients.push_back(std::real(coeff));
+            } else {
+              identityCoeff += std::real(coeff);
+              coefficients.push_back(std::real(coeff));
+            }
+          }
+        }
         // Retrieve instructions for gradient, if a pointer of type
         // AlgorithmGradientStrategy is given
         if (gradientStrategy) {
@@ -158,7 +179,11 @@ void VQE::execute(const std::shared_ptr<AcceleratorBuffer> buffer) const {
         }
 
         auto tmpBuffer = xacc::qalloc(buffer->size());
-        accelerator->execute(tmpBuffer, fsToExec);
+        if (cacheMeasurements) {
+          accelerator->execute(tmpBuffer, evaled, basisRotations);
+        } else {
+          accelerator->execute(tmpBuffer, fsToExec);
+        }
         auto buffers = tmpBuffer->getChildren();
 
         // Tag any gradient buffers;
@@ -170,7 +195,6 @@ void VQE::execute(const std::shared_ptr<AcceleratorBuffer> buffer) const {
         for (auto &[k, v] : tmp_buffer_extra_info) {
           buffer->addExtraInfo(k, v);
         }
-
         // Create buffer child for the Identity term
         auto idBuffer = xacc::qalloc(buffer->size());
         idBuffer->addExtraInfo("coefficient", identityCoeff);
@@ -181,15 +205,13 @@ void VQE::execute(const std::shared_ptr<AcceleratorBuffer> buffer) const {
         if (accelerator->name() == "ro-error")
           idBuffer->addExtraInfo("ro-fixed-exp-val-z", 1.0);
         buffer->appendChild("I", idBuffer);
-
         // Add information about the variational parameters to the child
         // buffers.
         // Other energy (observable-related) information will be populated by
         // the Observable::postProcess().
         for (auto &childBuffer : buffers) {
           childBuffer->addExtraInfo("parameters", x);
         }
-
         // Special key to indicate that the buffer was processed by a
         // HPC virtualization decorator.
         const std::string aggregate_key =
@@ -216,7 +238,6 @@ void VQE::execute(const std::shared_ptr<AcceleratorBuffer> buffer) const {
             return observable->postProcess(tmpBuffer);
           }
         }();
-
         // Compute the variance as well as populate any variance-related
         // information to the child buffers
         const double variance = [&]() {
@@ -233,7 +254,6 @@ void VQE::execute(const std::shared_ptr<AcceleratorBuffer> buffer) const {
                 tmpBuffer, Observable::PostProcessingTask::VARIANCE_CALC);
           }
         }();
-
         if (gradientStrategy) {
           // gradient-based optimization
           // If gradientStrategy is numerical, pass the energy
@@ -257,7 +277,6 @@ void VQE::execute(const std::shared_ptr<AcceleratorBuffer> buffer) const {
         for (auto &b : buffers) {
           buffer->appendChild(b->name(), b);
         }
-
         std::stringstream ss;
         ss << "E(" << (!x.empty() ? std::to_string(x[0]) : "");
         for (int i = 1; i < x.size(); i++)
@@ -275,7 +294,6 @@ void VQE::execute(const std::shared_ptr<AcceleratorBuffer> buffer) const {
           }
           last_energy = energy;
         }
-
         return energy;
       },
       kernel->nVariables());
@@ -314,31 +332,61 @@ VQE::execute(const std::shared_ptr<AcceleratorBuffer> buffer,
   std::vector<double> coefficients;
   std::vector<std::string> kernelNames;
   std::vector<std::shared_ptr<CompositeInstruction>> fsToExec;
-
   double identityCoeff = 0.0;
-  auto evaled = xacc::as_shared_ptr(kernel)->operator()(x);
-  auto kernels = observable->observe(evaled);
-  for (auto &f : kernels) {
-    kernelNames.push_back(f->name());
-    std::complex<double> coeff = f->getCoefficient();
-
-    int nFunctionInstructions = 0;
-    if (f->getInstruction(0)->isComposite()) {
-      nFunctionInstructions = kernel->nInstructions() + f->nInstructions() - 1;
-    } else {
-      nFunctionInstructions = f->nInstructions();
-    }
+  int nInstructionsEnergy, nInstructionsGradient = 0;
 
-    if (nFunctionInstructions > kernel->nInstructions()) {
-      fsToExec.push_back(f);
+  // call CompositeInstruction::operator()()
+  auto evaled = kernel->operator()(x);
+
+  // if we want to only deal with the measurement basis instead of entire
+  // circuits
+  if (cacheMeasurements) {
+
+    nInstructionsEnergy = basisRotations.size() - 1;
+    for (auto it = basisRotations.begin(); it != basisRotations.end();) {
+
+      kernelNames.push_back((*it)->name());
+      std::complex<double> coeff = (*it)->getCoefficient();
+      if ((*it)->name() == "I") {
+        identityCoeff += std::real(coeff);
+        it = basisRotations.erase(it);
+      }
       coefficients.push_back(std::real(coeff));
-    } else {
-      identityCoeff += std::real(coeff);
+      ++it;
+    }
+
+  } else {
+
+    // observe
+    auto kernels = observable->observe(evaled);
+    for (auto &f : kernels) {
+      kernelNames.push_back(f->name());
+      std::complex<double> coeff = f->getCoefficient();
+
+      int nFunctionInstructions;
+      if (f->getInstruction(0)->isComposite()) {
+        nFunctionInstructions =
+            kernel->nInstructions() + f->nInstructions() - 1;
+      } else {
+        nFunctionInstructions = f->nInstructions();
+      }
+
+      if (nFunctionInstructions > kernel->nInstructions()) {
+        fsToExec.push_back(f);
+        coefficients.push_back(std::real(coeff));
+      } else {
+        identityCoeff += std::real(coeff);
+        coefficients.push_back(std::real(coeff));
+      }
     }
   }
 
   auto tmpBuffer = xacc::qalloc(buffer->size());
-  accelerator->execute(tmpBuffer, fsToExec);
+  if (cacheMeasurements) {
+    accelerator->execute(tmpBuffer, evaled, basisRotations);
+  } else {
+    accelerator->execute(tmpBuffer, fsToExec);
+  }
   auto buffers = tmpBuffer->getChildren();
   for (auto &b : buffers) {
     b->addExtraInfo("parameters", x);

quantum/plugins/algorithms/vqe/vqe.hpp

@@ -26,7 +26,8 @@ class VQE : public Algorithm {
   Accelerator * accelerator;
   std::vector<double> initial_params;
   std::shared_ptr<AlgorithmGradientStrategy> gradientStrategy;
-
+  mutable std::vector<std::shared_ptr<CompositeInstruction>> basisRotations;
+  bool cacheMeasurements;
   HeterogeneousMap parameters;
 
 public:

quantum/plugins/qsim/accelerator/QsimAccelerator.cpp

@@ -11,6 +11,7 @@
  *   Thien Nguyen - initial API and implementation
  *******************************************************************************/
 #include "QsimAccelerator.hpp"
+#include "xacc.hpp"
 #include "xacc_plugin.hpp"
 #include "IRUtils.hpp"
 #include <cassert>
@@ -274,6 +275,56 @@ void QsimAccelerator::execute(
   }
 }
 
+void QsimAccelerator::execute(std::shared_ptr<AcceleratorBuffer> buffer,
+          const std::shared_ptr<CompositeInstruction> baseCircuit,
+          const std::vector<std::shared_ptr<CompositeInstruction>> basisRotations) {
+
+  constexpr int MAX_NUMBER_CIRCUITS_TO_ANALYZE = 100;
+  if (!m_vqeMode || basisRotations.size() <= 1 ||
+      basisRotations.size() > MAX_NUMBER_CIRCUITS_TO_ANALYZE) {
+    auto provider = getIRProvider("quantum");
+    // Cannot run VQE mode, just run each composite independently.
+    for (auto &b : basisRotations) {
+      auto f = provider->createComposite(b->name(), baseCircuit->getVariables());
+      auto tmpBuffer =
+          std::make_shared<xacc::AcceleratorBuffer>(b->name(), buffer->size());
+      execute(tmpBuffer, f);
+      buffer->appendChild(f->name(), tmpBuffer);
+    }
+  } else {
+    xacc::info("Running VQE mode");
+    QsimCircuitVisitor visitor(buffer->size());
+    // Walk the base IR tree, and visit each node
+    InstructionIterator it(baseCircuit);
+    while (it.hasNext()) {
+      auto nextInst = it.next();
+      if (nextInst->isEnabled() && !nextInst->isComposite()) {
+        nextInst->accept(&visitor);
+      }
+    }
+
+    // Run the base circuit:
+    auto circuit = visitor.getQsimCircuit();
+    StateSpace stateSpace(m_numThreads);
+    State state = stateSpace.Create(circuit.num_qubits);
+    stateSpace.SetStateZero(state);
+
+    const bool runOk = Runner::Run(m_qsimParam, Factory(m_numThreads), circuit, state);
+    assert(runOk);
+
+    // Now we have a wavefunction that represents execution of the ansatz.
+    // Run the observable sub-circuits (change of basis + measurements)
+    for (int i = 0; i < basisRotations.size(); ++i) {
+      auto tmpBuffer = std::make_shared<xacc::AcceleratorBuffer>(
+          basisRotations[i]->name(), buffer->size());
+      const double e = getExpectationValueZ(basisRotations[i], stateSpace, state);
+      tmpBuffer->addExtraInfo("exp-val-z", e);
+      buffer->appendChild(basisRotations[i]->name(), tmpBuffer);
+    }
+  }
+
+}
+
 double QsimAccelerator::getExpectationValueZ(
     std::shared_ptr<CompositeInstruction> compositeInstruction,
     const StateSpace &stateSpace, const State &state) const {

quantum/plugins/qsim/accelerator/QsimAccelerator.hpp

@@ -319,6 +319,9 @@ class QsimAccelerator : public Accelerator {
   virtual void execute(std::shared_ptr<AcceleratorBuffer> buffer,
                        const std::vector<std::shared_ptr<CompositeInstruction>>
                            compositeInstructions) override;
+  void execute(std::shared_ptr<AcceleratorBuffer> buffer,
+          const std::shared_ptr<CompositeInstruction> baseCircuit,
+          const std::vector<std::shared_ptr<CompositeInstruction>> basisRotations) override;
   virtual void apply(std::shared_ptr<AcceleratorBuffer> buffer,
                      std::shared_ptr<Instruction> inst) override;
 

xacc/accelerator/Accelerator.hpp

@@ -100,6 +100,13 @@ class Accelerator : public Identifiable {
                        const std::vector<std::shared_ptr<CompositeInstruction>>
                            CompositeInstructions) = 0;
 
+  virtual void
+  execute(std::shared_ptr<AcceleratorBuffer> buffer,
+          const std::shared_ptr<CompositeInstruction> baseCircuit,
+          const std::vector<std::shared_ptr<CompositeInstruction>> basisRotations) {
+    return;
+  }
+
   virtual void cancel(){};
 
   virtual std::vector<std::pair<int, int>> getConnectivity() {