Got the one-qubit base case worked

Need to enhance the Pauli algebra handler to support multiple qubits
Signed-off-by: Thien Nguyen <nguyentm@ornl.gov>

python/examples/qite_example.py

+# Simple 1-qubit demonstration of the Quatum Imaginary Time Evolution algorithm
+# Reference: https://www.nature.com/articles/s41567-019-0704-4
+# Target H = 1/sqrt(2)(X + Z)
+# Expected minimum value: -1.0
+import xacc,sys, numpy as np
+import matplotlib.pyplot as plt
+
+# Get access to the desired QPU and
+# allocate some qubits to run on
+qpu = xacc.getAccelerator('qpp')
+
+# Construct the Hamiltonian as an XACC PauliOperator
+ham = xacc.getObservable('pauli', '0.70710678118 X0 + 0.70710678118 Z0')
+
+# We just need 1 qubit
+buffer = xacc.qalloc(1)
+
+# See Fig. 2 (e) of https://www.nature.com/articles/s41567-019-0704-4
+# Horizontal axis: 0 -> 2.5
+# The number of Trotter steps 
+nbSteps = 25
+
+# The Trotter step size
+stepSize = 0.1
+
+# Create the QITE algorithm
+qite = xacc.getAlgorithm('qite', {
+                        'accelerator': qpu,
+                        'observable': ham,
+                        'step-size': stepSize,
+                        'steps': nbSteps
+                        })
+
+result = qite.execute(buffer)
+
+# Expected result: ~ -1
+print('Min energy value = ', buffer.getInformation('opt-val'))
+E = buffer.getInformation('exp-vals')
+# Reproduce Fig. 2(e) of https://www.nature.com/articles/s41567-019-0704-4
+plt.plot(np.arange(0, nbSteps + 1) * stepSize, E, 'ro', label = 'XACC QITE')
+plt.grid()
+plt.show()
\ No newline at end of file

quantum/plugins/algorithms/qite/qite.cpp

@@ -64,6 +64,8 @@ bool QITE::initialize(const HeterogeneousMap &parameters)
     m_observable = xacc::as_shared_ptr(parameters.getPointerLike<Observable>("observable"));
   }
 
+  m_approxOps.clear();
+  m_energyAtStep.clear();
   return initializeOk;
 }
 
@@ -118,10 +120,58 @@ std::shared_ptr<CompositeInstruction> QITE::constructPropagateCircuit() const
     }
   }
 
-  std::cout << "Progagated kernel:\n" << propagateKernel->toString() << "\n";
+  // std::cout << "Progagated kernel:\n" << propagateKernel->toString() << "\n";
   return propagateKernel;
 }
 
+double QITE::calcCurrentEnergy(int in_nbQubits) const
+{
+  // Trotter kernel up to this point
+  auto propagateKernel = constructPropagateCircuit();
+  auto kernels = m_observable->observe(propagateKernel);
+  std::vector<double> coefficients;
+  std::vector<std::string> kernelNames;
+  std::vector<std::shared_ptr<CompositeInstruction>> fsToExec;
+
+  double identityCoeff = 0.0;
+  for (auto &f : kernels) 
+  {
+    kernelNames.push_back(f->name());
+    std::complex<double> coeff = f->getCoefficient();
+    int nFunctionInstructions = 0;
+    if (f->getInstruction(0)->isComposite()) 
+    {
+      nFunctionInstructions = propagateKernel->nInstructions() + f->nInstructions() - 1;
+    } 
+    else 
+    {
+      nFunctionInstructions = f->nInstructions();
+    }
+
+    if (nFunctionInstructions > propagateKernel->nInstructions()) 
+    {
+      fsToExec.push_back(f);
+      coefficients.push_back(std::real(coeff));
+    } 
+    else 
+    {
+      identityCoeff += std::real(coeff);
+    }
+  }
+
+  auto tmpBuffer = xacc::qalloc(in_nbQubits);
+  m_accelerator->execute(tmpBuffer, fsToExec);
+  auto buffers = tmpBuffer->getChildren();
+
+  double energy = identityCoeff;
+  for (int i = 0; i < buffers.size(); ++i) 
+  {
+    auto expval = buffers[i]->getExpectationValueZ();
+    energy += expval * coefficients[i];
+  }
+  std::cout << "Energy = " << energy << "\n";
+  return energy;
+}
 
 void QITE::execute(const std::shared_ptr<AcceleratorBuffer> buffer) const 
 {
@@ -347,6 +397,8 @@ void QITE::execute(const std::shared_ptr<AcceleratorBuffer> buffer) const
     return updatedAham; 
   };
   
+  m_energyAtStep.emplace_back(calcCurrentEnergy(buffer->size()));
+  
   // Time stepping:
   for (int i = 0; i < m_nbSteps; ++i)
   {
@@ -358,12 +410,22 @@ void QITE::execute(const std::shared_ptr<AcceleratorBuffer> buffer) const
       auto nextAOps = calcAOps(kernel, hamTerm);
       m_approxOps.emplace_back(nextAOps);
     }
+    m_energyAtStep.emplace_back(calcCurrentEnergy(buffer->size()));
   }
+
+  assert(m_energyAtStep.size() == m_nbSteps + 1);
+  // Last energy value
+  buffer->addExtraInfo("opt-val", ExtraInfo(m_energyAtStep.back()));
+  // Also returns the full list of energy values 
+  // at each Trotter step.
+  buffer->addExtraInfo("exp-vals", ExtraInfo(m_energyAtStep));
 }
 
 std::vector<double> QITE::execute(const std::shared_ptr<AcceleratorBuffer> buffer, const std::vector<double>& x) 
 {
-  // TODO
+  // We don't do any parameter optimization here,
+  // hence don't support this!
+  xacc::error("This method is unsupported!");
   return {};
 }
 } // namespace algorithm

quantum/plugins/algorithms/qite/qite.hpp

@@ -30,7 +30,9 @@ public:
 private:
   // Construct the Trotter propagate circuit up to current time step.
   std::shared_ptr<CompositeInstruction> constructPropagateCircuit() const;
-
+  // Calculate the current energy, i.e.
+  // the value of the observable at the current Trotter step.
+  double calcCurrentEnergy(int in_nbQubits) const;
 private:
   // Number of Trotter steps
   int m_nbSteps;
@@ -47,6 +49,8 @@ private:
   // of the Hamiltonian observable
   // i.e. exp(-iAt) -> exp(-Ht)
   mutable std::vector<std::shared_ptr<Observable>> m_approxOps;
+  // Energy value achieved at a Trotter step
+  mutable std::vector<double> m_energyAtStep;
 };
 } // namespace algorithm
 } // namespace xacc

quantum/plugins/algorithms/qite/tests/QITETester.cpp

@@ -39,6 +39,13 @@ TEST(QITETester, checkSimple)
   
   auto buffer = xacc::qalloc(1);
   qite->execute(buffer);
+  const double finalEnergy = (*buffer)["opt-val"].as<double>();
+  std::cout << "Final Energy: " << finalEnergy << "\n";
+  // Fig (2.e) of https://www.nature.com/articles/s41567-019-0704-4
+  // Minimal Energy = -1
+  EXPECT_NEAR(finalEnergy, -1.0, 1e-3);
+  const std::vector<double> energyValues = (*buffer)["exp-vals"].as<std::vector<double>>();
+  EXPECT_EQ(energyValues.size(), nbSteps + 1);
 }
 
 int main(int argc, char **argv) {