Add an option to set the tolerance for jordan_wigner (NVIDIA#23)

add option to set jordan-wigner's tolerance

Signed-off-by: Melody Ren <[email protected]>

libs/core/include/cuda-qx/core/heterogeneous_map.h

@@ -46,10 +46,6 @@ class heterogeneous_map {
   /// @param _other The map to copy from
   heterogeneous_map(const heterogeneous_map &_other) { *this = _other; }
 
-  /// @brief Move constructor
-  /// @param _other The map to move from
-  heterogeneous_map(heterogeneous_map &_other) { *this = _other; }
-
   /// @brief Constructor from initializer list
   /// @param list The initializer list of key-value pairs
   heterogeneous_map(
@@ -65,8 +61,10 @@ class heterogeneous_map {
   /// @param _other The map to assign from
   /// @return Reference to this map
   heterogeneous_map &operator=(const heterogeneous_map &_other) {
-    clear();
-    items = _other.items;
+    if (this != &_other) {
+      clear();
+      items = _other.items;
+    }
     return *this;
   }
 
@@ -76,20 +74,14 @@ class heterogeneous_map {
   /// @param value The value
   template <typename T>
   void insert(const std::string &key, const T &value) {
-    auto iter = items.find(key);
-    if (iter == items.end()) {
+
+    if constexpr (is_bounded_char_array<T>{}) {
       // Never insert a raw char array or char ptr,
       // auto conver to a string
-      if constexpr (is_bounded_char_array<T>{}) {
-        items.insert({key, std::string(value)});
-        return;
-      }
-
-      items.insert({key, value});
-      return;
+      items.insert_or_assign(key, std::string(value));
+    } else {
+      items.insert_or_assign(key, value);
     }
-
-    items.at(key) = value;
   }
 
   /// @brief Get a value from the map

libs/solvers/include/cudaq/solvers/operators/molecule/fermion_compiler.h

@@ -9,6 +9,7 @@
 #pragma once
 
 #include "cuda-qx/core/extension_point.h"
+#include "cuda-qx/core/heterogeneous_map.h"
 #include "cuda-qx/core/tensor.h"
 
 #include "cudaq/spin_op.h"
@@ -23,9 +24,10 @@ class fermion_compiler : public cudaqx::extension_point<fermion_compiler> {
 public:
   /// @brief Given a fermionic representation of an operator
   /// generate an equivalent operator on spins.
-  virtual cudaq::spin_op generate(const double constant,
-                                  const cudaqx::tensor<> &hpq,
-                                  const cudaqx::tensor<> &hpqrs) = 0;
+  virtual cudaq::spin_op
+  generate(const double constant, const cudaqx::tensor<> &hpq,
+           const cudaqx::tensor<> &hpqrs,
+           const cudaqx::heterogeneous_map &options = {}) = 0;
   virtual ~fermion_compiler() {}
 };
 

libs/solvers/include/cudaq/solvers/operators/molecule/fermion_compilers/jordan_wigner.h

@@ -15,7 +15,8 @@ namespace cudaq::solvers {
 class jordan_wigner : public fermion_compiler {
 public:
   cudaq::spin_op generate(const double constant, const cudaqx::tensor<> &hpq,
-                          const cudaqx::tensor<> &hpqrs) override;
+                          const cudaqx::tensor<> &hpqrs,
+                          const cudaqx::heterogeneous_map &options) override;
 
   CUDAQ_EXTENSION_CREATOR_FUNCTION(fermion_compiler, jordan_wigner)
 };

libs/solvers/lib/operators/molecule/fermion_compilers/jordan_wigner.cpp

@@ -382,12 +382,16 @@ cudaq::spin_op two_body(std::size_t p, std::size_t q, std::size_t r,
 
 cudaq::spin_op jordan_wigner::generate(const double constant,
                                        const tensor<> &hpq,
-                                       const tensor<> &hpqrs) {
+                                       const tensor<> &hpqrs,
+                                       const heterogeneous_map &options) {
   assert(hpq.rank() == 2 && "hpq must be a rank-2 tensor");
   assert(hpqrs.rank() == 4 && "hpqrs must be a rank-4 tensor");
   auto spin_hamiltonian = constant * cudaq::spin_op();
   std::size_t nqubit = hpq.shape()[0];
-  double tolerance = 1e-15;
+
+  double tolerance =
+      options.get<double>(std::vector<std::string>{"tolerance", "tol"}, 1e-15);
+
   for (auto p : cudaq::range(nqubit)) {
     auto coef = hpq.at({p, p});
     if (std::fabs(coef) > tolerance)

libs/solvers/lib/operators/operator_pools/uccsd_operator_pool.cpp

@@ -18,7 +18,8 @@ using excitation_list = std::vector<std::vector<std::size_t>>;
 std::vector<cudaq::spin_op>
 uccsd::generate(const heterogeneous_map &config) const {
 
-  auto numQubits = config.get<int>({"num-qubits", "num_qubits"});
+  auto numQubits =
+      config.get<int>({"num-qubits", "num_qubits", "n-qubits", "n_qubits"});
   auto numElectrons = config.get<int>({"num-electrons", "num_electrons"});
   std::size_t spin = 0;
   if (config.contains("spin"))

libs/solvers/python/bindings/solvers/py_solvers.cpp

@@ -203,7 +203,8 @@ void bindOperators(py::module &mod) {
 
   mod.def(
       "jordan_wigner",
-      [](py::buffer hpq, py::buffer hpqrs, double core_energy = 0.0) {
+      [](py::buffer hpq, py::buffer hpqrs, double core_energy = 0.0,
+         py::kwargs options) {
         auto hpqInfo = hpq.request();
         auto hpqrsInfo = hpqrs.request();
         auto *hpqData = reinterpret_cast<std::complex<double> *>(hpqInfo.ptr);
@@ -216,7 +217,7 @@ void bindOperators(py::module &mod) {
                       {hpqrsInfo.shape.begin(), hpqrsInfo.shape.end()});
 
         return fermion_compiler::get("jordan_wigner")
-            ->generate(core_energy, hpqT, hpqrsT);
+            ->generate(core_energy, hpqT, hpqrsT, hetMapFromKwargs(options));
       },
       py::arg("hpq"), py::arg("hpqrs"), py::arg("core_energy") = 0.0,
       R"#(
@@ -235,6 +236,11 @@ hpqrs : numpy.ndarray
     Shape should be (N, N, N, N) where N is the number of spin molecular orbitals.
 core_energy : float, optional
     The core energy of the system when using active space Hamiltonian, nuclear energy otherwise. Default is 0.0.
+tolerance : float, optional
+    The threshold value for ignoring small coefficients.
+    Can also be specified using 'tol'.
+    Coefficients with absolute values smaller than this tolerance are considered as zero.
+    Default is 1e-15.
 
 Returns:
 --------
@@ -254,7 +260,7 @@ RuntimeError
 >>> h1 = np.array([[0, 1], [1, 0]], dtype=np.complex128)
 >>> h2 = np.zeros((2, 2, 2, 2), dtype=np.complex128)
 >>> h2[0, 1, 1, 0] = h2[1, 0, 0, 1] = 0.5
->>> qubit_op = jordan_wigner(h1, h2, core_energy=0.1)
+>>> qubit_op = jordan_wigner(h1, h2, core_energy=0.1, tolerance=1e-14)
 
 Notes:
 ------
@@ -267,7 +273,7 @@ RuntimeError
 
   mod.def(
       "jordan_wigner",
-      [](py::buffer buffer, double core_energy = 0.0) {
+      [](py::buffer buffer, double core_energy = 0.0, py::kwargs options) {
         auto info = buffer.request();
         auto *data = reinterpret_cast<std::complex<double> *>(info.ptr);
         std::size_t size = 1;
@@ -279,14 +285,14 @@ RuntimeError
           cudaqx::tensor hpq, hpqrs({dim, dim, dim, dim});
           hpq.borrow(data, {info.shape.begin(), info.shape.end()});
           return fermion_compiler::get("jordan_wigner")
-              ->generate(core_energy, hpq, hpqrs);
+              ->generate(core_energy, hpq, hpqrs, hetMapFromKwargs(options));
         }
 
         std::size_t dim = info.shape[0];
         cudaqx::tensor hpq({dim, dim}), hpqrs;
         hpqrs.borrow(data, {info.shape.begin(), info.shape.end()});
         return fermion_compiler::get("jordan_wigner")
-            ->generate(core_energy, hpq, hpqrs);
+            ->generate(core_energy, hpq, hpqrs, hetMapFromKwargs(options));
       },
       py::arg("hpq"), py::arg("core_energy") = 0.0,
       R"#(
@@ -304,6 +310,11 @@ hpq : numpy.ndarray
     where N is the number of orbitals.
 core_energy : float, optional
     The core energy of the system. Default is 0.0.
+tolerance : float, optional
+    The threshold value for ignoring small coefficients.
+    Can also be specified using 'tol'.
+    Coefficients with absolute values smaller than this tolerance are considered as zero.
+    Default is 1e-15.
 
 Returns:
 --------
@@ -322,7 +333,7 @@ RuntimeError
 >>> import numpy as np
 >>> # One-body integrals
 >>> h1 = np.array([[0, 1], [1, 0]], dtype=np.complex128)
->>> qubit_op1 = jordan_wigner(h1, core_energy=0.1)
+>>> qubit_op1 = jordan_wigner(h1, core_energy=0.1, tolerance=1e-14)
 
 >>> # Two-body integrals
 >>> h2 = np.zeros((2, 2, 2, 2), dtype=np.complex128)

libs/solvers/python/tests/test_molecule.py

@@ -53,8 +53,11 @@ def test_jordan_wigner():
                                        0,
                                        verbose=True,
                                        casci=True)
-    op = solvers.jordan_wigner(molecule.hpq, molecule.hpqrs,
-                               molecule.energies['nuclear_energy'])
+
+    op = solvers.jordan_wigner(molecule.hpq,
+                               molecule.hpqrs,
+                               core_energy=molecule.energies['nuclear_energy'],
+                               tol=1e-15)
     assert molecule.hamiltonian == op
     hpq = np.array(molecule.hpq)
     hpqrs = np.array(molecule.hpqrs)
@@ -115,7 +118,9 @@ def test_jordan_wigner_as():
 
     hpq = np.array(molecule.hpq)
     hpqrs = np.array(molecule.hpqrs)
-    hpqJw = solvers.jordan_wigner(hpq, molecule.energies['core_energy'])
+    hpqJw = solvers.jordan_wigner(hpq,
+                                  core_energy=molecule.energies['core_energy'],
+                                  tolerance=1e-15)
     hpqrsJw = solvers.jordan_wigner(hpqrs)
     op2 = hpqJw + hpqrsJw