committing work on Observable::to_sparse_matrix() plus PauliOperator implementation using armadillo

Signed-off-by: Alex McCaskey <mccaskeyaj@ornl.gov>

quantum/annealing/ir/AnnealingProgram.hpp

@@ -22,7 +22,7 @@ namespace xacc {
 namespace quantum {
 
 class AnnealingProgram : public CompositeInstruction,
-                   public std::enable_shared_from_this<AnnealingProgram> {
+                         public std::enable_shared_from_this<AnnealingProgram> {
 
 protected:
   std::vector<InstPtr> instructions;
@@ -61,15 +61,17 @@ protected:
     print_backtrace();
     exit(0);
   }
+
 public:
-  AnnealingProgram(std::string kernelName)
-      : _name(kernelName) {}
+  AnnealingProgram(std::string kernelName) : _name(kernelName) {}
 
   AnnealingProgram(std::string kernelName, std::vector<std::string> p)
       : _name(kernelName), variables(p) {}
 
-
   void applyRuntimeArguments() override {
+    for (auto &i : instructions) {
+      i->applyRuntimeArguments();
+    }
   }
 
   std::shared_ptr<CompositeInstruction> enabledView() override {
@@ -84,25 +86,28 @@ public:
   }
 
   // ising or qubo
-  const std::string getTag() override {return tag;}
-  void setTag(const std::string& t) override {tag = t;return;}
+  const std::string getTag() override { return tag; }
+  void setTag(const std::string &t) override {
+    tag = t;
+    return;
+  }
 
   const int nInstructions() override { return instructions.size(); }
-  const int nChildren() override {return instructions.size();}
+  const int nChildren() override { return instructions.size(); }
 
   void mapBits(std::vector<std::size_t> bitMap) override {
     xacc::error("AnnealingProgrma.mapBits not implemented");
   }
 
-  void removeDisabled() override {
-  }
+  void removeDisabled() override {}
 
- void setBitExpression(const std::size_t bit_idx, const std::string expr) override {}
- std::string getBitExpression(const std::size_t bit_idx) override {return "";}
-  const int nRequiredBits() const override {
-    return 0;
+  void setBitExpression(const std::size_t bit_idx,
+                        const std::string expr) override {}
+  std::string getBitExpression(const std::size_t bit_idx) override {
+    return "";
   }
-InstPtr getInstruction(const std::size_t idx) override {
+  const int nRequiredBits() const override { return 0; }
+  InstPtr getInstruction(const std::size_t idx) override {
     validateInstructionIndex(idx);
     return instructions[idx];
   }
@@ -140,9 +145,7 @@ InstPtr getInstruction(const std::size_t idx) override {
     xacc::error("AnnealingProgram graph is undirected, cannot compute depth.");
     return 0;
   }
-  void clear() override {
-      instructions.clear();
-  }
+  void clear() override { instructions.clear(); }
 
   const std::string persistGraph() override {
     std::stringstream s;
@@ -152,26 +155,24 @@ InstPtr getInstruction(const std::size_t idx) override {
 
   std::shared_ptr<Graph> toGraph() override;
 
-
   const std::string name() const override { return _name; }
   const std::string description() const override { return ""; }
-  void setName(const std::string name) override {
+  void setName(const std::string name) override {}
+  const std::vector<std::size_t> bits() override {
+    return std::vector<std::size_t>{};
   }
-  const std::vector<std::size_t> bits() override { return std::vector<std::size_t>{}; }
 
   const std::string toString() override {
-       std::stringstream ss;
+    std::stringstream ss;
     for (auto i : instructions) {
       ss << i->toString() << ";\n";
     }
     return ss.str();
-   }
+  }
 
   const std::size_t nLogicalBits() override { return 0; }
   const std::size_t nPhysicalBits() override { return 0; }
-  const std::set<std::size_t> uniqueBits() override {
-      return {};
-  }
+  const std::set<std::size_t> uniqueBits() override { return {}; }
 
   std::vector<double> getAllBiases() {
     std::vector<double> biases;
@@ -184,9 +185,10 @@ InstPtr getInstruction(const std::size_t idx) override {
     return biases;
   }
 
-  std::string getBufferName(const std::size_t bitIdx) override {return "";}
-  void setBufferNames(const std::vector<std::string> bufferNamesPerIdx) override {}
-  std::vector<std::string> getBufferNames() override {return {};}
+  std::string getBufferName(const std::size_t bitIdx) override { return ""; }
+  void
+  setBufferNames(const std::vector<std::string> bufferNamesPerIdx) override {}
+  std::vector<std::string> getBufferNames() override { return {}; }
 
   std::vector<double> getAllCouplers() {
     std::vector<double> weights;
@@ -202,7 +204,7 @@ InstPtr getInstruction(const std::size_t idx) override {
   void persist(std::ostream &outStream) override;
   void load(std::istream &inStream) override;
 
-   const InstructionParameter
+  const InstructionParameter
   getParameter(const std::size_t idx) const override {
     errorCircuitParameter();
     return InstructionParameter(0);
@@ -224,18 +226,13 @@ InstPtr getInstruction(const std::size_t idx) override {
   std::shared_ptr<CompositeInstruction>
   operator()(const std::vector<double> &params) override;
 
-  void setCoefficient(const std::complex<double> c) override {
-  }
+  void setCoefficient(const std::complex<double> c) override {}
   const std::complex<double> getCoefficient() override { return 1.0; }
-    const std::vector<std::string> requiredKeys() override {
-    return {};
-  }
+  const std::vector<std::string> requiredKeys() override { return {}; }
   void setBits(const std::vector<std::size_t> bits) override {}
 
   bool hasChildren() const override { return !instructions.empty(); }
-  bool expand(const HeterogeneousMap &runtimeOptions) override {
-    return true;
-  }
+  bool expand(const HeterogeneousMap &runtimeOptions) override { return true; }
   void addVariable(const std::string variableName) override {
     variables.push_back(variableName);
   }
@@ -252,12 +249,10 @@ InstPtr getInstruction(const std::size_t idx) override {
 
   const std::size_t nVariables() override { return variables.size(); }
 
-  const std::string accelerator_signature() override {
-    return "dwave";
-  }
+  const std::string accelerator_signature() override { return "dwave"; }
   void set_accelerator_signature(const std::string signature) override {}
 
-std::shared_ptr<Instruction> clone() override {
+  std::shared_ptr<Instruction> clone() override {
     return std::make_shared<AnnealingProgram>(*this);
   }
   EMPTY_DEFINE_VISITABLE()

quantum/gate/GateQuantumActivator.cpp

@@ -50,6 +50,7 @@ public:
     auto z = std::make_shared<xacc::quantum::Z>();
     auto sw = std::make_shared<xacc::quantum::Swap>();
     auto u = std::make_shared<xacc::quantum::U>();
+    auto anneal = std::make_shared<xacc::quantum::AnnealingInstruction>();
 
     auto s = std::make_shared<xacc::quantum::S>();
     auto sdg = std::make_shared<xacc::quantum::Sdg>();
@@ -80,6 +81,7 @@ public:
     context.RegisterService<xacc::Instruction>(sw);
     context.RegisterService<xacc::Instruction>(u);
     context.RegisterService<xacc::Instruction>(ifstmt);
+    context.RegisterService<xacc::Instruction>(anneal);
 
     context.RegisterService<xacc::Instruction>(s);
     context.RegisterService<xacc::Instruction>(sdg);

quantum/gate/ir/CommonGates.hpp

@@ -58,8 +58,7 @@ public:
       for (auto &i : instructions) {
         i->enable();
       }
-    }
-    else {
+    } else {
       // Note: although sub-instructions are initially disabled,
       // we need to disable here as well just in case we run multiple shots
       // and they may be enabled in the previous run.
@@ -116,6 +115,50 @@ public:
   DEFINE_VISITABLE()
 };
 
+class AnnealingInstruction : public Gate {
+public:
+  AnnealingInstruction()
+      : Gate("AnnealingInstruction",
+             std::vector<InstructionParameter>{InstructionParameter(0.0)}) {}
+  AnnealingInstruction(std::size_t qbit, InstructionParameter &&bias)
+      : Gate("AnnealingInstruction", std::vector<std::size_t>{qbit},
+             std::vector<InstructionParameter>{bias}) {}
+  AnnealingInstruction(std::size_t qbit1, std::size_t qbit2,
+                       InstructionParameter &&weight)
+      : Gate("AnnealingInstruction", std::vector<std::size_t>{qbit1, qbit2},
+             std::vector<InstructionParameter>{weight}) {}
+  AnnealingInstruction(std::vector<std::size_t> qbits)
+      : Gate("AnnealingInstruction", qbits,
+             std::vector<InstructionParameter>{InstructionParameter(0.0)}) {
+    // make it easier to specify single qubit bias
+    if (qbits.size() == 1) {
+      qbits.push_back(qbits[0]);
+    }
+  }
+  void setBits(const std::vector<std::size_t> bts) override {
+    if (bts.size() == 1) {
+      qbits.push_back(bts[0]);
+      qbits.push_back(bts[0]);
+    } else {
+        qbits = bts;
+    }
+  }
+  const int nRequiredBits() const override { return 2; }
+  void
+  setBufferNames(const std::vector<std::string> bufferNamesPerIdx) override {
+    std::vector<std::string> tmp;
+    tmp.push_back(bufferNamesPerIdx[0]);
+    if (bufferNamesPerIdx.size() == 1) {
+      tmp.push_back(bufferNamesPerIdx[0]);
+    } else {
+      tmp.push_back(bufferNamesPerIdx[1]);
+    }
+    Gate::setBufferNames(tmp);
+  }
+  DEFINE_CLONE(AnnealingInstruction)
+  DEFINE_VISITABLE()
+};
+
 class Rx : public Gate {
 public:
   Rx()
@@ -193,7 +236,8 @@ public:
       : Gate("U", std::vector<InstructionParameter>{
                       InstructionParameter(0.0), InstructionParameter(0.0),
                       InstructionParameter(0.0)}) {}
-  U(std::size_t qbit, std::vector<xacc::InstructionParameter> params) : Gate("U", params){}
+  U(std::size_t qbit, std::vector<xacc::InstructionParameter> params)
+      : Gate("U", params) {}
   U(std::size_t qbit, double theta, double phi, double lambda)
       : Gate("U", std::vector<std::size_t>{qbit},
              std::vector<InstructionParameter>{InstructionParameter(theta),

quantum/observable/pauli/CMakeLists.txt

@@ -25,7 +25,7 @@ usFunctionGenerateBundleInit(TARGET ${LIBRARY_NAME} OUT SRC)
 add_library(${LIBRARY_NAME} SHARED ${SRC})
 target_include_directories(${LIBRARY_NAME}
                             PUBLIC
-                            .
+                            . ${CMAKE_SOURCE_DIR}/tpls/armadillo
                             ${CMAKE_SOURCE_DIR}/tpls/antlr/runtime/src
                             ${CMAKE_CURRENT_SOURCE_DIR}
                             ${CMAKE_CURRENT_SOURCE_DIR}/generated
@@ -54,7 +54,7 @@ usFunctionEmbedResources(TARGET
                          manifest.json)
 
 if(APPLE)
-	target_link_libraries(${LIBRARY_NAME} PUBLIC xacc ${CMAKE_SOURCE_DIR}/dist/libantlr4-runtime.dylib)
+   target_link_libraries(${LIBRARY_NAME} PUBLIC xacc ${CMAKE_SOURCE_DIR}/dist/libantlr4-runtime.dylib)
    set_target_properties(${LIBRARY_NAME} PROPERTIES INSTALL_RPATH "@loader_path")
    set_target_properties(${LIBRARY_NAME} PROPERTIES LINK_FLAGS "-undefined dynamic_lookup")
 else()

quantum/observable/pauli/PauliOperator.cpp

@@ -15,6 +15,7 @@
 #include <regex>
 #include <set>
 #include <iostream>
+#include "Observable.hpp"
 #include "xacc.hpp"
 #include "xacc_service.hpp"
 
@@ -23,10 +24,92 @@
 #include "PauliOperatorLexer.h"
 #include "PauliListenerImpl.hpp"
 
+#include <armadillo>
+
 namespace xacc {
 namespace quantum {
-// const std::map<std::string, std::pair<c, std::string>> Term:: pauliProducts =
-// Term::create_map();
+
+std::vector<SparseTriplet> PauliOperator::to_sparse_matrix() {
+  auto n_qubits = nQubits();
+  auto n_hilbert = std::pow(2, n_qubits);
+  using SparseMatrix = arma::SpMat<std::complex<double>>;
+
+  SparseMatrix x(2, 2), y(2, 2), z(2, 2);
+  x(0, 1) = 1.0;
+  x(1, 0) = 1.0;
+  y(0, 1) = std::complex<double>(0, -1);
+  y(1, 0) = std::complex<double>(0, 1);
+  z(0, 0) = 1.;
+  z(1, 1) = -1.;
+
+  SparseMatrix i = arma::speye<SparseMatrix>(2, 2);
+
+  std::map<std::string, SparseMatrix> mat_map{
+      {"I", i}, {"X", x}, {"Y", y}, {"Z", z}};
+
+  auto kron_ops = [](std::vector<SparseMatrix> &ops) {
+    auto first = ops[0];
+    for (int i = 1; i < ops.size(); i++) {
+      first = arma::kron(first, ops[i]);
+    }
+    return first;
+  };
+
+  SparseMatrix total(n_hilbert, n_hilbert);
+  for (auto &term : terms) {
+    auto tensor_factor = 0;
+    auto coeff = term.second.coeff();
+
+    std::vector<SparseMatrix> sparse_mats;
+
+    if (term.second.ops().empty()) {
+      // this was I term
+      auto id = arma::speye<SparseMatrix>(n_hilbert, n_hilbert);
+      sparse_mats.push_back(id);
+    } else {
+      for (auto &pauli : term.second.ops()) {
+        if (pauli.first > tensor_factor) {
+
+          auto id_qbits = pauli.first - tensor_factor;
+          auto id = arma::speye<SparseMatrix>((int)std::pow(2, id_qbits),
+                                              (int)std::pow(2, id_qbits));
+          sparse_mats.push_back(id);
+        }
+
+        sparse_mats.push_back(mat_map[pauli.second]);
+        tensor_factor = pauli.first + 1;
+      }
+
+      for (int i = tensor_factor; i < n_qubits; i++) {
+        auto id = arma::speye<SparseMatrix>(2, 2);
+        sparse_mats.push_back(id);
+      }
+    }
+
+    
+    auto sp_matrix = kron_ops(sparse_mats);
+    sp_matrix *= coeff;
+    total += sp_matrix;
+  }
+
+//   arma::vec eigval;
+//   arma::mat eigvec;
+
+//   arma::sp_mat test(total.n_rows, total.n_cols);
+//   for (auto i = total.begin(); i != total.end(); ++i) {
+//     test(i.row(), i.col()) = (*i).real();
+//   }
+
+//   arma::eigs_sym(eigval, eigvec, test, 1);
+
+//   std::cout << "EIGS:\n" << eigval << "\n";
+
+  std::vector<SparseTriplet> trips;
+  for (auto iter = total.begin(); iter != total.end(); ++iter) {
+    trips.emplace_back(iter.row(), iter.col(), *iter);
+  }
+  return trips;
+}
 
 PauliOperator::PauliOperator() {}
 
@@ -119,9 +202,9 @@ PauliOperator::observe(std::shared_ptr<CompositeInstruction> function) {
     }
 
     for (auto arg : function->getArguments()) {
-       gateFunction->addArgument(arg, 0);
+      gateFunction->addArgument(arg, 0);
     }
-    
+
     // Loop over all terms in the Spin Instruction
     // and create instructions to run on the Gate QPU.
     std::vector<std::shared_ptr<xacc::Instruction>> measurements;
@@ -191,29 +274,6 @@ Term::toBinaryVector(const int nQubits) {
   return {v, w};
 }
 
-std::vector<Triplet> PauliOperator::getSparseMatrixElements() {
-
-  // Get number of qubits
-  std::set<int> distinctSites;
-  for (auto &kv : terms) {
-    for (auto &kv2 : kv.second.ops()) {
-      distinctSites.insert(kv2.first);
-    }
-  }
-
-  auto nQubits = distinctSites.size();
-
-  std::vector<Triplet> triplets;
-
-  for (auto &kv : terms) {
-    auto termTrips = kv.second.getSparseMatrixElements(nQubits);
-    triplets.insert(std::end(triplets), std::begin(termTrips),
-                    std::end(termTrips));
-  }
-
-  return triplets;
-}
-
 ActionResult Term::action(const std::string &bitString, ActionType type) {
 
   auto _coeff = coeff();
@@ -470,7 +530,7 @@ PauliOperator::operator*=(const std::complex<double> v) noexcept {
   return *this;
 }
 
-std::vector<Triplet> Term::getSparseMatrixElements(const int nQubits) {
+std::vector<SparseTriplet> Term::getSparseMatrixElements(const int nQubits) {
 
   // X = |1><0| + |0><1|
   // Y = -i|1><0| + i|0><1|
@@ -532,7 +592,7 @@ std::vector<Triplet> Term::getSparseMatrixElements(const int nQubits) {
   auto ket = zeroStr;
   auto bra = zeroStr;
 
-  std::vector<Triplet> triplets;
+  std::vector<SparseTriplet> triplets;
   for (auto &combo : termCombinations) {
 
     std::complex<double> coeff(1, 0), i(0, 1);

quantum/observable/pauli/PauliOperator.hpp

@@ -47,17 +47,6 @@ using TermTuple =
 using c = std::complex<double>;
 using ActionResult = std::pair<std::string, c>;
 enum ActionType { Bra, Ket };
-class Triplet : std::tuple<std::uint64_t, std::uint64_t, std::complex<double>> {
-public:
-  Triplet(std::uint64_t r, std::uint64_t c, std::complex<double> coeff) {
-    std::get<0>(*this) = r;
-    std::get<1>(*this) = c;
-    std::get<2>(*this) = coeff;
-  }
-  const std::uint64_t row() { return std::get<0>(*this); }
-  const std::uint64_t col() { return std::get<1>(*this); }
-  const std::complex<double> coeff() { return std::get<2>(*this); }
-};
 
 class Term : public TermTuple,
              public tao::operators::commutative_multipliable<Term>,
@@ -213,7 +202,7 @@ public:
     return (std::get<1>(*this) == std::get<1>(v) && ops() == std::get<2>(v));
   }
 
-  std::vector<Triplet> getSparseMatrixElements(const int nQubits);
+  std::vector<SparseTriplet> getSparseMatrixElements(const int nQubits);
 
   ActionResult action(const std::string &bitString, ActionType type);
 
@@ -326,9 +315,12 @@ public:
 
   std::unordered_map<std::string, Term> getTerms() const { return terms; }
 
-  std::vector<Triplet> getSparseMatrixElements();
+  std::vector<SparseTriplet> getSparseMatrixElements() {return to_sparse_matrix();}
   std::vector<std::complex<double>> toDenseMatrix(const int nQubits);
 
+  std::vector<SparseTriplet>
+  to_sparse_matrix() override;
+
   std::shared_ptr<IR> toXACCIR();
   void fromXACCIR(std::shared_ptr<IR> ir);
   PauliOperator

quantum/plugins/dwave/CMakeLists.txt

@@ -50,6 +50,7 @@ target_include_directories(${LIBRARY_NAME}
 target_link_libraries(${LIBRARY_NAME}
                       PUBLIC xacc
                              xacc-quantum-annealing
+                             xacc-quantum-gate
                              ${ANTLR_LIB} dwave_sapi
                              )
 

quantum/plugins/dwave/DWaveActivator.cpp

@@ -45,6 +45,8 @@ public:
 
     auto rbm = std::make_shared<xacc::dwave::RBM>();
     context.RegisterService<xacc::Instruction>(rbm);
+    auto rbmascirc = std::make_shared<xacc::dwave::RBMAsCircuitType>();
+    context.RegisterService<xacc::Instruction>(rbmascirc);
     // auto c = std::make_shared<xacc::quantum::DWQMICompiler>();
     // context.RegisterService<xacc::Compiler>(c);
     // context.RegisterService<xacc::OptionsProvider>(c);

quantum/plugins/dwave/generators/rbm.cpp

@@ -11,10 +11,13 @@
  *   Alexander J. McCaskey - initial API and implementation
  *******************************************************************************/
 #include "rbm.hpp"
+#include "CommonGates.hpp"
+#include "Instruction.hpp"
 #include "xacc.hpp"
 #include "xacc_service.hpp"
 
 #include "DWQMI.hpp"
+#include "CommonGates.hpp"
 
 using namespace xacc;
 namespace xacc {
@@ -62,5 +65,62 @@ bool RBM::expand(const xacc::HeterogeneousMap &runtimeOptions) {
   return true;
 } // namespace instructions
 
+void RBMAsCircuitType::applyRuntimeArguments() {
+
+  // We expect args like this
+  // rbm(v, x, nv, nh);
+
+  std::string variable_name = arguments[0]->name;
+
+  std::vector<double> x_vals;
+  if (arguments[0]->type.find("std::vector<double>") != std::string::npos) {
+    x_vals = arguments[0]->runtimeValue.get<std::vector<double>>(
+        INTERNAL_ARGUMENT_VALUE_KEY);
+  } else {
+    xacc::error("Has to be a vector of double");
+  }
+
+  auto n_visible = arguments[1]->runtimeValue.get<int>(INTERNAL_ARGUMENT_VALUE_KEY);
+  auto n_hidden = arguments[2]->runtimeValue.get<int>(INTERNAL_ARGUMENT_VALUE_KEY);
+
+  std::cout << "HELLO WORLD: " << x_vals << ", " << n_visible << ", " << n_hidden << "\n";
+
+  if (nInstructions() > 0) {
+      clear();
+      parameters.clear();
+  }
+
+  for (std::size_t i = 0; i < n_visible; i++) {
+    // std::string paramName = "v" + std::to_string(i);
+    // addVariable(paramName);
+    auto qmi = std::make_shared<xacc::quantum::AnnealingInstruction>(i, x_vals[i]);
+    addInstruction(qmi);
+  }
+
+  for (std::size_t i = n_visible; i < n_visible + n_hidden; i++) {
+    // std::string paramName = "h" + std::to_string(i);
+    // addVariable(paramName);
+    // xacc::InstructionParameter hidParam(paramName);
+    auto qmi = std::make_shared<xacc::quantum::AnnealingInstruction>(i, x_vals[i]);
+    addInstruction(qmi);
+  }
+
+  int counter = n_visible+n_hidden;
+  for (std::size_t i = 0; i < n_visible; i++) {
+    for (std::size_t j = n_visible; j < n_visible + n_hidden; j++) {
+    //   std::string paramName = "w" + std::to_string(i) + std::to_string(j);
+    //   addVariable(paramName);
+    //   xacc::InstructionParameter wParam(paramName);
+      auto qmi = std::make_shared<xacc::quantum::AnnealingInstruction>(i, j, x_vals[counter]);
+      addInstruction(qmi);
+      counter++;
+    }
+  }
+
+
+  parameters.emplace_back(n_visible);
+  parameters.emplace_back(n_hidden);
+}
+
 } // namespace dwave
 } // namespace xacc
\ No newline at end of file

quantum/plugins/dwave/generators/rbm.hpp

@@ -14,6 +14,7 @@
 #define XACC_DWAVE_GENERATORS_RBM_HPP_
 
 #include "AnnealingProgram.hpp"
+#include "Circuit.hpp"
 
 namespace xacc {
 namespace dwave {
@@ -22,7 +23,7 @@ protected:
   std::vector<InstructionParameter> parameters;
 
 public:
-  RBM() : AnnealingProgram("rbm") {}
+  RBM() : AnnealingProgram("rbm-ap") {}
   bool expand(const xacc::HeterogeneousMap &runtimeOptions) override;
   const std::vector<std::string> requiredKeys() override;
  std::shared_ptr<Instruction> clone() override {
@@ -34,6 +35,21 @@ public:
     return parameters[idx];
   }