Merge pull request #25 from kharazit/mixed_op

Mixed operator algebra methods for fermion ops and pauli ops.

examples/fermi/fermi_op.cpp

+#include "qcor.hpp"
+using namespace qcor;
+__qpu__ void ansatz(qreg q, double theta) {
+  X(q[0]);
+  X(q[2]);
+  Rx(q[0], 1.5708);
+  H(q[1]);
+  H(q[2]);
+  H(q[3]);
+  CNOT(q[0], q[1]);
+  CNOT(q[1], q[2]);
+  CNOT(q[2], q[3]);
+  Rz(q[3], theta);
+  CNOT(q[2], q[3]);
+  CNOT(q[1], q[2]);
+  CNOT(q[0], q[1]);
+  Rx(q[0], -1.5708);
+  H(q[1]);
+  H(q[2]);
+  H(q[3]);
+}
+
+int main(int argc, char **argv) {
+  auto q = qalloc(4);
+
+  auto H = -0.165606823582 * adag(1) * adag(2) * a(1) * a(2) +
+           0.120200490713 * adag(1) * adag(0) * a(0) * a(1) +
+           -0.0454063328691 * adag(0) * adag(3) * a(1) * a(2) +
+           0.168335986252 * adag(2) * adag(0) * a(0) * a(2) +
+           0.0454063328691 * adag(1) * adag(2) * a(3) * a(0) +
+           0.168335986252 * adag(0) * adag(2) * a(2) * a(0) +
+           0.165606823582 * adag(0) * adag(3) * a(3) * a(0) -
+           0.0454063328691 * adag(3) * adag(0) * a(2) * a(1) +
+           -0.0454063328691 * adag(1) * adag(3) * a(0) * a(2) -
+           0.0454063328691 * adag(3) * adag(1) * a(2) * a(0) +
+           0.165606823582 * adag(1) * adag(2) * a(2) * a(1) -
+           0.165606823582 * adag(0) * adag(3) * a(0) * a(3) -
+           0.0454063328691 * adag(1) * adag(2) * a(0) * a(3) -
+           0.479677813134 * adag(3) * a(3) -
+           0.174072892497 * adag(1) * adag(3) * a(1) * a(3) -
+           0.0454063328691 * adag(0) * adag(2) * a(1) * a(3) +
+           0.120200490713 * adag(0) * adag(1) * a(1) * a(0) +
+           0.0454063328691 * adag(0) * adag(2) * a(3) * a(1) +
+           0.174072892497 * adag(1) * adag(3) * a(3) * a(1) +
+           0.165606823582 * adag(2) * adag(1) * a(1) * a(2) +
+           -0.0454063328691 * adag(2) * adag(1) * a(3) * a(0) -
+           0.120200490713 * adag(2) * adag(3) * a(2) * a(3) +
+           0.120200490713 * adag(2) * adag(3) * a(3) * a(2) -
+           0.168335986252 * adag(0) * adag(2) * a(0) * a(2) +
+           0.120200490713 * adag(3) * adag(2) * a(2) * a(3) -
+           0.120200490713 * adag(3) * adag(2) * a(3) * a(2) +
+           0.0454063328691 * adag(1) * adag(3) * a(2) * a(0) -
+           1.2488468038 * adag(0) * a(0) +
+           0.0454063328691 * adag(3) * adag(1) * a(0) * a(2) -
+           0.168335986252 * adag(2) * adag(0) * a(2) * a(0) +
+           0.0165606823582 * adag(3) * adag(0) * a(0) * a(3) +
+           2.0 * 0.0454063328691 * adag(2) * adag(0) * a(1) * a(3) -
+           1.2488468038 * adag(2) * a(2) +
+           0.0454063328691 * adag(2) * adag(1) * a(0) * a(3) +
+           0.174072892497 * adag(3) * adag(1) * a(1) * a(3) -
+           0.174072892497 * adag(3) * adag(1) * a(3) * a(1) +
+           0.0454063328691 * adag(3) * adag(0) * a(1) * a(2) -
+           0.479677813134 * adag(1) * a(1) -
+           0.165606823582 * adag(3) * adag(0) * a(3) * a(0) +
+           0.0454063328691 * adag(0) * adag(3) * a(2) * a(1) +
+           -0.165606823582 * adag(2) * adag(1) * a(2) * a(1) -
+           0.120200490713 * adag(0) * adag(1) * a(0) * a(1) +
+           -0.120200490713 * adag(1) * adag(0) * a(1) * a(0) + 0.7080240981;
+
+  auto objective = qcor::createObjectiveFunction("vqe", ansatz, H);
+  auto optimizer = qcor::createOptimizer("nlopt");
+  qcor::OptFunction f(
+      [&](const std::vector<double> &x, std::vector<double> &grad) {
+        return (*objective)(q, x[0]);
+      },
+      10);
+  auto results = optimizer->optimize(f);
+  printf("vqe energy = %f\n", results.first);
+  std::cout << results.second[0] << std::endl;
+  return 0;
+}
+
+/*
+    auto H =  createObservable("fermion", "(-0.165606823582,-0)  1^ 2^ 1 2 +
+   (0.120200490713,0)  1^ 0^ 0 1 + "
+                          "(-0.0454063328691,-0)  0^ 3^ 1 2 + (0.168335986252,0)
+   2^ 0^ 0 2 + "
+                          "(0.0454063328691,0)  1^ 2^ 3 0 + (0.168335986252,0)
+   0^ 2^ 2 0 + "
+                          "(0.165606823582,0)  0^ 3^ 3 0 + (-0.0454063328691,-0)
+   3^ 0^ 2 1 + "
+                          "(-0.0454063328691,-0)  1^ 3^ 0 2 +
+   (-0.0454063328691,-0)  3^ 1^ 2 0 + "
+                          "(0.165606823582,0)  1^ 2^ 2 1 + (-0.165606823582,-0)
+   0^ 3^ 0 3 + "
+                          "(-0.479677813134,-0)  3^ 3 + (-0.0454063328691,-0) 1^
+   2^ 0 3 + "
+                          "(-0.174072892497,-0)  1^ 3^ 1 3 +
+   (-0.0454063328691,-0)  0^ 2^ 1 3 + "
+                          "(0.120200490713,0)  0^ 1^ 1 0 + (0.0454063328691,0)
+   0^ 2^ 3 1 + "
+                          "(0.174072892497,0)  1^ 3^ 3 1 + (0.165606823582,0) 2^
+   1^ 1 2 + "
+                          "(-0.0454063328691,-0)  2^ 1^ 3 0 +
+   (-0.120200490713,-0)  2^ 3^ 2 3 + "
+                          "(0.120200490713,0)  2^ 3^ 3 2 + (-0.168335986252,-0)
+   0^ 2^ 0 2 + "
+                          "(0.120200490713,0)  3^ 2^ 2 3 + (-0.120200490713,-0)
+   3^ 2^ 3 2 + "
+                          "(0.0454063328691,0)  1^ 3^ 2 0 + (-1.2488468038,-0)
+   0^ 0 + "
+                          "(0.0454063328691,0)  3^ 1^ 0 2 + (-0.168335986252,-0)
+   2^ 0^ 2 0 + "
+                          "(0.165606823582,0)  3^ 0^ 0 3 + (-0.0454063328691,-0)
+   2^ 0^ 3 1 + "
+                          "(0.0454063328691,0)  2^ 0^ 1 3 + (-1.2488468038,-0)
+   2^ 2 + "
+                          "(0.0454063328691,0)  2^ 1^ 0 3 + (0.174072892497,0)
+   3^ 1^ 1 3 + "
+                          "(-0.479677813134,-0)  1^ 1 + (-0.174072892497,-0)  3^
+   1^ 3 1 + "
+                          "(0.0454063328691,0)  3^ 0^ 1 2 + (-0.165606823582,-0)
+   3^ 0^ 3 0 + "
+                          "(0.0454063328691,0)  0^ 3^ 2 1 + (-0.165606823582,-0)
+   2^ 1^ 2 1 + "
+                          "(-0.120200490713,-0)  0^ 1^ 0 1 +
+   (-0.120200490713,-0)  1^ 0^ 1 0 + (0.7080240981,0)");
+
+*/
\ No newline at end of file

examples/fermi/mixed_op.cpp

+#include "qcor.hpp"
+
+using namespace qcor;
+
+__qpu__ void ansatz(qreg q, double theta, std::shared_ptr<xacc::Observable> H) {
+  X(q[0]);
+  X(q[1]);
+}
+
+int main(int argc, char **argv) {
+  using qcor::X; // Due to ambiguous call with xacc::quantum::X
+
+  auto q = qalloc(2);
+
+  std::cout << "Fermi op " << std::endl;
+
+  auto fermi_H = adag(0) * a(1) + adag(1) * a(0);
+
+  std::cout << fermi_H.toString() << std::endl;
+
+  std::cout << "fermi -> pauli transf " << std::endl;
+
+  auto pauli_H = transform(fermi_H);
+
+  std::cout << pauli_H.toString() << std::endl;
+
+  std::cout << "mixed sum: " << std::endl;
+
+  auto mixed_sum = fermi_H + pauli_H;
+
+  std::cout << mixed_sum.toString() << std::endl;
+
+  std::cout << "mixed product " << std::endl;
+
+  auto mixed_prod = adag(0) * X(1) + X(0) * a(1);
+
+  std::cout << mixed_prod.toString() << std::endl;
+
+  auto H = X(0) + X(0) * a(1) + adag(1) * X(0) + X(0) * adag(0) * adag(1);
+
+  std::cout << "Mixed sum and products: " << std::endl;
+
+  std::cout << H.toString() << std::endl;
+}
\ No newline at end of file

runtime/CMakeLists.txt

@@ -6,7 +6,7 @@ add_library(${LIBRARY_NAME} SHARED ${SRC})
 
 target_include_directories(${LIBRARY_NAME} PUBLIC . qrt)
 
-target_link_libraries(${LIBRARY_NAME} PUBLIC xacc::xacc xacc::quantum_gate qrt xacc::pauli)
+target_link_libraries(${LIBRARY_NAME} PUBLIC xacc::xacc xacc::quantum_gate qrt utils)
 
 xacc_configure_library_rpath(${LIBRARY_NAME})
 
@@ -19,4 +19,5 @@ if (QCOR_BUILD_TESTS)
 endif()
 
 add_subdirectory(objectives)
-add_subdirectory(qrt)
\ No newline at end of file
+add_subdirectory(qrt)
+add_subdirectory(utils)
\ No newline at end of file

runtime/qcor.cpp

@@ -2,12 +2,12 @@
 
 #include "Optimizer.hpp"
 
+#include "qalloc.hpp"
+#include "qrt.hpp"
 #include "xacc.hpp"
 #include "xacc_quantum_gate_api.hpp"
 #include "xacc_service.hpp"
-
-#include "qalloc.hpp"
-#include "qrt.hpp"
+#include <memory>
 
 namespace qcor {
 void set_verbose(bool verbose) { xacc::set_verbose(verbose); }
@@ -63,6 +63,7 @@ double observe(std::shared_ptr<CompositeInstruction> program, Observable &obs,
 }
 } // namespace __internal__
 
+
 std::shared_ptr<xacc::Optimizer> createOptimizer(const std::string &type,
                                                  HeterogeneousMap &&options) {
   if (!xacc::isInitialized())
@@ -76,6 +77,13 @@ std::shared_ptr<xacc::Observable> createObservable(const std::string &repr) {
   return xacc::quantum::getObservable("pauli", std::string(repr));
 }
 
+std::shared_ptr<Observable> createObservable(const std::string &type,
+                                             const std::string &repr) {
+  if (!xacc::isInitialized())
+    xacc::internal_compiler::compiler_InitializeXACC();
+  return xacc::quantum::getObservable(type, std::string(repr));
+}
+
 std::shared_ptr<xacc::CompositeInstruction> compile(const std::string &src) {
   return xacc::getCompiler("xasm")->compile(src)->getComposites()[0];
 }

runtime/qcor.hpp

@@ -10,14 +10,12 @@
 #include <vector>
 
 #include "CompositeInstruction.hpp"
-#include "Observable.hpp"
 #include "Optimizer.hpp"
-
-#include "PauliOperator.hpp"
+#include "Observable.hpp"
 #include "qalloc"
-#include "xacc_internal_compiler.hpp"
-
 #include "qrt.hpp"
+#include "operator_algebra.hpp"
+#include "xacc_internal_compiler.hpp"
 
 namespace qcor {
 
@@ -26,46 +24,6 @@ using HeterogeneousMap = xacc::HeterogeneousMap;
 using Observable = xacc::Observable;
 using Optimizer = xacc::Optimizer;
 using CompositeInstruction = xacc::CompositeInstruction;
-using PauliOperator = xacc::quantum::PauliOperator;
-
-PauliOperator X(int idx) { return PauliOperator({{idx, "X"}}); }
-
-PauliOperator Y(int idx) { return PauliOperator({{idx, "Y"}}); }
-
-PauliOperator Z(int idx) { return PauliOperator({{idx, "Z"}}); }
-
-PauliOperator allZs(const int nQubits) {
-  auto ret = Z(0);
-  for (int i = 1; i < nQubits; i++) {
-    ret *= Z(i);
-  }
-  return ret;
-}
-
-template <typename T> PauliOperator operator+(T coeff, PauliOperator &op) {
-  return PauliOperator(coeff) + op;
-}
-template <typename T> PauliOperator operator+(PauliOperator &op, T coeff) {
-  return PauliOperator(coeff) + op;
-}
-
-template <typename T> PauliOperator operator-(T coeff, PauliOperator &op) {
-  return -1.0 * coeff + op;
-}
-
-template <typename T> PauliOperator operator-(PauliOperator &op, T coeff) {
-  return -1.0 * coeff + op;
-}
-
-PauliOperator SP(int idx) {
-  std::complex<double> imag(0.0, 1.0);
-  return X(idx) + imag * Y(idx);
-}
-
-PauliOperator SM(int idx) {
-  std::complex<double> imag(0.0, 1.0);
-  return X(idx) - imag * Y(idx);
-}
 
 class ResultsBuffer {
 public:
@@ -406,6 +364,8 @@ createOptimizer(const std::string &type, HeterogeneousMap &&options = {});
 // Create an observable from a string representation
 std::shared_ptr<Observable> createObservable(const std::string &repr);
 
+std::shared_ptr<Observable> createObservable(const std::string &type,
+                                             const std::string &repr);
 std::shared_ptr<ObjectiveFunction> createObjectiveFunction(
     const std::string &obj_name, std::shared_ptr<CompositeInstruction> kernel,
     std::shared_ptr<Observable> observable, HeterogeneousMap &&options = {}) {

runtime/utils/CMakeLists.txt

+set(LIBRARY_NAME utils)
+
+file(GLOB SRC *.cpp)
+
+add_library(${LIBRARY_NAME} SHARED ${SRC})
+
+target_include_directories(${LIBRARY_NAME} PUBLIC .)
+
+target_link_libraries(${LIBRARY_NAME} PUBLIC xacc::xacc xacc::quantum_gate xacc::pauli xacc::fermion)
+
+xacc_configure_library_rpath(${LIBRARY_NAME})
+
+file(GLOB HEADERS operator_algebra.hpp)
+install(FILES ${HEADERS} DESTINATION include/qcor)
+install(TARGETS ${LIBRARY_NAME} DESTINATION lib)

runtime/utils/operator_algebra.cpp

+#include "operator_algebra.hpp"
+#include "xacc.hpp"
+#include "xacc_service.hpp"
+#include "xacc_internal_compiler.hpp"
+
+
+namespace qcor{
+    PauliOperator transform(FermionOperator &obs, std::string transf) {
+  if (!xacc::isInitialized())
+    xacc::internal_compiler::compiler_InitializeXACC();
+  auto obsv = std::make_shared<FermionOperator>(obs);
+  auto terms = std::dynamic_pointer_cast<PauliOperator>(
+      xacc::getService<xacc::ObservableTransform>(transf)->transform(obsv));
+  return *terms;
+}
+}
+
+

runtime/utils/operator_algebra.hpp

+#ifndef RUNTIME_QCOR_OP_ALG_HPP_
+#define RUNTIME_QCOR_OP_ALG_HPP_
+
+#include "qalloc.hpp"
+#include <CompositeInstruction.hpp>
+#include <memory>
+
+#include "FermionOperator.hpp"
+#include "ObservableTransform.hpp"
+#include "PauliOperator.hpp"
+#include "Observable.hpp"
+
+namespace qcor{
+using PauliOperator = xacc::quantum::PauliOperator;
+using FermionOperator = xacc::quantum::FermionOperator;
+
+PauliOperator X(int idx) { return PauliOperator({{idx, "X"}}); }
+PauliOperator Y(int idx) { return PauliOperator({{idx, "Y"}}); }
+PauliOperator Z(int idx) { return PauliOperator({{idx, "Z"}}); }
+PauliOperator SP(int idx) {
+  std::complex<double> imag(0.0, 1.0);
+  return X(idx) + imag * Y(idx);
+}
+PauliOperator SM(int idx) {
+  std::complex<double> imag(0.0, 1.0);
+  return X(idx) - imag * Y(idx);
+}
+FermionOperator a(int idx) {
+  std::string s("(1.0, 0) " + std::to_string(idx));
+  return FermionOperator(s);
+}
+FermionOperator adag(int idx) {
+  std::string s("(1.0, 0) " + std::to_string(idx) + "^");
+  return FermionOperator(s);
+}
+
+PauliOperator allZs(const int nQubits) {
+  auto ret = Z(0);
+  for (int i = 1; i < nQubits; i++) {
+    ret *= Z(i);
+  }
+  return ret;
+}
+
+// transform FermionOperator to PauliOperator
+PauliOperator transform(FermionOperator &obs, std::string transf = "jw");
+
+template <typename T, typename = typename std::enable_if<
+                          std::is_arithmetic<T>::value, T>::type>
+PauliOperator operator+(T coeff, PauliOperator &op) {
+  return PauliOperator(coeff) + op;
+}
+template <typename T, typename = typename std::enable_if<
+                          std::is_arithmetic<T>::value, T>::type>
+PauliOperator operator+(PauliOperator &op, T coeff) {
+  return PauliOperator(coeff) + op;
+}
+template <typename T, typename = typename std::enable_if<
+                          std::is_arithmetic<T>::value, T>::type>
+PauliOperator operator-(T coeff, PauliOperator &op) {
+  return -1.0 * coeff + op;
+}
+template <typename T, typename = typename std::enable_if<
+                          std::is_arithmetic<T>::value, T>::type>
+PauliOperator operator-(PauliOperator &op, T coeff) {
+  return -1.0 * coeff + op;
+}
+template <typename T, typename = typename std::enable_if<
+                          std::is_arithmetic<T>::value, T>::type>
+FermionOperator operator+(T coeff, FermionOperator &op) {
+  return FermionOperator(coeff) + op;
+}
+template <typename T, typename = typename std::enable_if<
+                          std::is_arithmetic<T>::value, T>::type>
+FermionOperator operator+(FermionOperator &op, T coeff) {
+  return FermionOperator(coeff) + op;
+}
+template <typename T, typename = typename std::enable_if<
+                          std::is_arithmetic<T>::value, T>::type>
+FermionOperator operator-(T coeff, FermionOperator &op) {
+  return -1.0 * coeff + op;
+}
+template <typename T, typename = typename std::enable_if<
+                          std::is_arithmetic<T>::value, T>::type>
+FermionOperator operator-(FermionOperator &op, T coeff) {
+  return -1.0 * coeff + op;
+}
+
+template <typename T, typename = typename std::enable_if<
+                          std::is_arithmetic<T>::value, T>::type>
+PauliOperator operator+(T coeff, PauliOperator &&op) {
+  return PauliOperator(coeff) + op;
+}
+template <typename T, typename = typename std::enable_if<
+                          std::is_arithmetic<T>::value, T>::type>
+PauliOperator operator+(PauliOperator &&op, T coeff) {
+  return PauliOperator(coeff) + op;
+}
+template <typename T, typename = typename std::enable_if<
+                          std::is_arithmetic<T>::value, T>::type>
+PauliOperator operator-(T coeff, PauliOperator &&op) {
+  return -1.0 * coeff + op;
+}
+template <typename T, typename = typename std::enable_if<
+                          std::is_arithmetic<T>::value, T>::type>
+PauliOperator operator-(PauliOperator &&op, T coeff) {
+  return -1.0 * coeff + op;
+}
+template <typename T, typename = typename std::enable_if<
+                          std::is_arithmetic<T>::value, T>::type>
+FermionOperator operator+(T coeff, FermionOperator &&op) {
+  return FermionOperator(coeff) + op;
+}
+template <typename T, typename = typename std::enable_if<
+                          std::is_arithmetic<T>::value, T>::type>
+FermionOperator operator+(FermionOperator &&op, T coeff) {
+  return FermionOperator(coeff) + op;
+}
+template <typename T, typename = typename std::enable_if<
+                          std::is_arithmetic<T>::value, T>::type>
+FermionOperator operator-(T coeff, FermionOperator &&op) {
+  return -1.0 * coeff + op;
+}
+template <typename T, typename = typename std::enable_if<
+                          std::is_arithmetic<T>::value, T>::type>
+FermionOperator operator-(FermionOperator &&op, T coeff) {
+  return -1.0 * coeff + op;
+}
+
+PauliOperator operator+(FermionOperator &&fop, PauliOperator &&pop) {
+  auto pfop = transform(fop);
+  return pfop + pop;
+}
+
+PauliOperator operator+(PauliOperator &&pop, FermionOperator &&fop) {
+  auto pfop = transform(fop);
+  return pop + pfop;
+}
+
+PauliOperator operator*(PauliOperator &&pop, FermionOperator &&fop) {
+  auto pfop = transform(fop);
+  return pfop * pop;
+}
+
+PauliOperator operator*(FermionOperator &&fop, PauliOperator &&pop) {
+  auto pfop = transform(fop);
+  return pop * pfop;
+}
+
+PauliOperator operator+(FermionOperator &fop, PauliOperator &pop) {
+  auto pfop = transform(fop);
+  return pfop + pop;
+}
+
+PauliOperator operator+(PauliOperator &pop, FermionOperator &fop) {
+  auto pfop = transform(fop);
+  return pop + pfop;
+}
+
+PauliOperator operator*(PauliOperator &pop, FermionOperator &fop) {
+  auto pfop = transform(fop);
+  return pfop * pop;
+}
+
+PauliOperator operator*(FermionOperator &fop, PauliOperator &pop) {
+  auto pfop = transform(fop);
+  return pop * pfop;
+}
+
+}
+#endif
\ No newline at end of file

tools/driver/qcor.in

@@ -8,7 +8,7 @@ def main(argv=None):
                             '-L', '@CMAKE_INSTALL_PREFIX@/lib', '-lxacc',
                             '-lqcor', '-lqrt', '-lqcor-hybrid', '-lCppMicroServices',
                             '-lxacc-quantum-gate',
-                            '-lxacc-pauli', '-lpthread']
+                            '-lxacc-pauli', '-lpthread', '-lxacc-fermion']
     baseIncludes = ['-I', '@CMAKE_INSTALL_PREFIX@/include/xacc', '-I', '@CMAKE_INSTALL_PREFIX@/include/qcor', '-I', '@CMAKE_INSTALL_PREFIX@/include/quantum/gate']
     defaultFlags = ['-std=c++17', '-fplugin=@CMAKE_INSTALL_PREFIX@/clang-plugins/libqcor-syntax-handler.so']