Feature/pauli string cpp

CMakeLists.txt

@@ -1,13 +1,27 @@
 #
 # Boilerplate CMakeLists.txt for C++ projects
 #
-cmake_minimum_required(VERSION 3.25)
+cmake_minimum_required(VERSION 3.27)
 
 set(CMAKE_EXPORT_COMPILE_COMMANDS
     TRUE
     CACHE BOOL "Export compile commands to build directory" FORCE)
 
 include(cmake/CPM.cmake)
+
+# Set C++ standard
+set(CMAKE_CXX_STANDARD 20)
+set(CMAKE_CXX_STANDARD_REQUIRED ON)
+set(CMAKE_CXX_EXTENSIONS OFF)
+# Need to enforce -fPIC across whole project to build shared libraries
+set(CMAKE_POSITION_INDEPENDENT_CODE ON)
+
+# TODO different compilation flags for Debug/Release modes; also pull in
+# target_compile_options for fast_pauli target defined below
+if(NOT CMAKE_BUILD_TYPE)
+  set(CMAKE_BUILD_TYPE Release)
+endif()
+
 # TODO adding a package lock to help with SBOM
 # cpmusepackagelock(package-lock.cmake)
 
@@ -38,11 +52,6 @@ cpmaddpackage("gh:kokkos/mdspan#b885a2c60ad42f9e1aaa0d317a38105b950cbed0")
 
 project(fast_pauli LANGUAGES CXX)
 
-# Set C++ standard
-set(CMAKE_CXX_STANDARD 20)
-set(CMAKE_CXX_STANDARD_REQUIRED ON)
-set(CMAKE_CXX_EXTENSIONS OFF)
-
 # Set up OpenMP
 find_package(OpenMP REQUIRED)
 
@@ -66,8 +75,10 @@ target_compile_options(
 # Testing
 include(CTest)
 enable_testing()
-# TODO use proper variable for project root
 add_subdirectory(fast_pauli/cpp/tests)
 
 # Examples
 add_subdirectory(fast_pauli/cpp/examples)
+
+# Python bindings
+add_subdirectory(fast_pauli/cpp/src)

Makefile

@@ -26,6 +26,8 @@ lint:
 
 # run with make -i to ignore errors and run all three formatters
 format:
+	pre-commit run --all-files -- cmake-format
+	pre-commit run --all-files -- clang-format
 	pre-commit run --all-files -- ruff-format
 	pre-commit run --all-files -- yamlfmt
 	pre-commit run --all-files -- trailing-whitespace

README.md

@@ -41,9 +41,11 @@ pre-commit install # installs the checks as pre-commit hooks
 
 ```bash
 cmake -B build -DCMAKE_CXX_COMPILER=clang++
-cmake --build build
+cmake --build build --parallel
+cmake --install build
 ctest --test-dir build
 ```
+Compiled `_fast_pauli` python module gets installed into `fast_pauli` directory.
 
 ## Design Choices
 

fast_pauli/cpp/include/__pauli.hpp

@@ -26,7 +26,7 @@ struct Pauli {
    * @brief Default constructor, initializes to I.
    *
    */
-  Pauli() : code(0) {}
+  constexpr Pauli() : code(0) {}
 
   /**
    * @brief Constructor given a numeric code.
@@ -37,11 +37,36 @@ struct Pauli {
    */
   template <class T>
     requires std::convertible_to<T, uint8_t>
-  Pauli(T const code) : code(code) {
+  constexpr Pauli(T const code) : code(code) {
     if (code < 0 || code > 3)
       throw std::invalid_argument("Pauli code must be 0, 1, 2, or 3");
   }
 
+  /**
+   * @brief Constructor given Pauli matrix symbol.
+   *
+   * @param symbol pauli matrix - I, X, Y, Z
+   * @return
+   */
+  constexpr Pauli(char const symbol) {
+    switch (symbol) {
+    case 'I':
+      this->code = 0;
+      break;
+    case 'X':
+      this->code = 1;
+      break;
+    case 'Y':
+      this->code = 2;
+      break;
+    case 'Z':
+      this->code = 3;
+      break;
+    default:
+      throw std::invalid_argument("Invalid Pauli matrix symbol");
+    }
+  }
+
   // Copy ctor
   Pauli(Pauli const &other) = default;
 
@@ -58,7 +83,7 @@ struct Pauli {
    * @param rhs
    * @return  std::pair<std::complex<double>, Pauli>
    */
-  friend std::pair<std::complex<double>, Pauli> operator*(Pauli lhs,
+  friend std::pair<std::complex<double>, Pauli> operator*(Pauli const &lhs,
                                                           Pauli const &rhs) {
     switch (lhs.code) {
     case 0:

fast_pauli/cpp/include/__pauli_string.hpp

@@ -1,12 +1,10 @@
 #ifndef __PAULI_STRING_HPP
 #define __PAULI_STRING_HPP
 
-#include <cstddef>
 #include <fmt/format.h>
 #include <fmt/ranges.h>
 
 #include <algorithm>
-#include <cstring>
 #include <experimental/mdspan>
 #include <ranges>
 #include <string>
@@ -71,7 +69,8 @@ struct PauliString {
         weight += 1;
         break;
       default:
-        throw std::invalid_argument(std::string("Invalid Pauli character") + c);
+        throw std::invalid_argument(std::string("Invalid Pauli character ") +
+                                    c);
       }
     }
   }
@@ -114,21 +113,25 @@ struct PauliString {
    *
    * PauliStrings are always sparse and have only single non-zero element per
    * row. It's N non-zero elements for NxN matrix where N is 2^n_qubits.
-   * Therefore j, k, and m will always have N elements.
+   * Therefore k and m will always have N elements.
    *
-   * TODO remove j because it's unused (and redundant).
    * See Algorithm 1 in https://arxiv.org/pdf/2301.00560.pdf for details about
    * the algorithm.
    *
    * @tparam T The floating point type (i.e. std::complex<T> for the values of
    * the PauliString matrix)
-   * @param j (unused)
    * @param k The column index of the matrix
    * @param m The values of the matrix
    */
   template <std::floating_point T>
-  void get_sparse_repr(std::vector<size_t> &j, std::vector<size_t> &k,
+  void get_sparse_repr(std::vector<size_t> &k,
                        std::vector<std::complex<T>> &m) const {
+    // TODO should we simply return a tuple of k,m instead of modifying args?
+    // for now there are no use-cases where we'd re-use already allocated k, m
+    // TODO also should get_sparse_repr be a non-member helper function
+    // that operates on PauliString?
+    // gonna be like: auto [k, m] = fast_pauli::internal::get_sparse_repr(ps);
+
     // We reverse the view here because the tensor product is taken from right
     // to left
     auto ps = paulis | std::views::reverse;
@@ -139,8 +142,6 @@ struct PauliString {
     size_t const dim = 1 << n;
 
     // Safe, but expensive, we overwrite the vectors
-    j = std::vector<size_t>(dim);
-    // j.clear();
     k = std::vector<size_t>(dim);
     m = std::vector<std::complex<T>>(dim);
 
@@ -153,13 +154,27 @@ struct PauliString {
         return 1UL;
       }
     };
+    // Helper function that resolves first value of pauli string
+    auto inital_value = [&nY]() -> std::complex<T> {
+      switch (nY % 4) {
+      case 0:
+        return 1.0;
+      case 1:
+        return {0.0, -1.0};
+      case 2:
+        return -1.0;
+      case 3:
+        return {0.0, 1.0};
+      }
+      return {};
+    };
 
     // Populate the initial values of our output
     k[0] = 0;
     for (size_t i = 0; i < ps.size(); ++i) {
       k[0] += (1UL << i) * diag(ps[i]);
     }
-    m[0] = std::pow(-1i, nY % 4);
+    m[0] = inital_value();
 
     // Populate the rest of the values in a recursive-like manner
     for (size_t l = 0; l < n; ++l) {
@@ -170,19 +185,18 @@ struct PauliString {
         eps = -1;
       }
 
-      T sign = std::pow(-1., diag(po)); // Keep this outside the loop
+      T sign = diag(po) ? -1.0 : 1.0;
 
-      for (size_t li = (1UL << l); li < (1UL << (l + 1)); li++) {
-        k[li] = k[li - (1UL << l)] + (1UL << l) * sign;
-        m[li] = m[li - (1UL << l)] * eps;
+      auto const lower_bound = 1UL << l;
+      for (size_t li = lower_bound; li < (lower_bound << 1); li++) {
+        k[li] = k[li - lower_bound] + lower_bound * sign;
+        m[li] = m[li - lower_bound] * eps;
       }
     }
   }
 
   /**
-   * @brief Apply the PauliString (using the sparse representation) to a vector.
-   * This performs following matrix-vector multiplication \f$ \mathcal{\hat{P}}
-   * \ket{\psi} \f$
+   * @brief @copybrief PauliString::apply(std::mdspan)
    *
    * @tparam T The floating point base to use for all the complex numbers
    * @param v The input vector to apply the PauliString to. Must be the same
@@ -193,26 +207,14 @@ struct PauliString {
   template <std::floating_point T>
   std::vector<std::complex<T>>
   apply(std::vector<std::complex<T>> const &v) const {
-    // Input check
-    if (v.size() != dims()) {
-      throw std::invalid_argument(
-          "Input vector size must match the number of qubits");
-    }
-
-    std::vector<size_t> j, k; // TODO reminder that j is unused
-    std::vector<std::complex<T>> m;
-    get_sparse_repr(j, k, m);
-
-    std::vector<std::complex<T>> result(v.size(), 0);
-    for (size_t i = 0; i < k.size(); ++i) {
-      result[i] += m[i] * v[k[i]];
-    }
-
-    return result;
+    // route this to implementation we have for mdspan specialization
+    return this->apply(std::mdspan(v.data(), v.size()));
   }
 
   /**
-   * @brief @copybrief PauliString::apply(std::vector<std::complex<T>>)
+   * @brief Apply the PauliString (using the sparse representation) to a vector.
+   * This performs following matrix-vector multiplication \f$ \mathcal{\hat{P}}
+   * \ket{\psi} \f$
    *
    * @tparam T The floating point base to use for all the complex numbers
    * @param v The input vector to apply the PauliString to. Must be the same
@@ -222,16 +224,16 @@ struct PauliString {
    */
   template <std::floating_point T>
   std::vector<std::complex<T>>
-  apply(std::mdspan<std::complex<T>, std::dextents<size_t, 1>> v) const {
+  apply(std::mdspan<const std::complex<T>, std::dextents<size_t, 1>> v) const {
     // Input check
     if (v.size() != dims()) {
       throw std::invalid_argument(
           "Input vector size must match the number of qubits");
     }
 
-    std::vector<size_t> j, k; // TODO reminder that j is unused
+    std::vector<size_t> k;
     std::vector<std::complex<T>> m;
-    get_sparse_repr(j, k, m);
+    get_sparse_repr(k, m);
 
     std::vector<std::complex<T>> result(v.size(), 0);
     for (size_t i = 0; i < k.size(); ++i) {
@@ -253,9 +255,9 @@ struct PauliString {
    *
    * @tparam T The floating point base to use for all the complex numbers
    * @param new_states_T The output states after applying the PauliString
-   * (n_data x n_dim)
-   * @param states_T THe original states to apply the PauliString to (n_data x
-   * n_dim)
+   * (n_dim x n_states)
+   * @param states_T THe original states to apply the PauliString to
+   * (n_dim x n_states)
    * @param c Multiplication factor to apply to the PauliString
    */
   template <std::floating_point T>
@@ -279,20 +281,15 @@ struct PauliString {
           "[PauliString] new_states must have the same dimensions as states");
     }
 
-    std::vector<size_t> j, k; // TODO reminder that j is unused
+    std::vector<size_t> k;
     std::vector<std::complex<T>> m;
-    get_sparse_repr(j, k, m);
+    get_sparse_repr(k, m);
 
-    // TODO we have bad memory access patterns
-    // std::vector<std::complex<T>> col(states_T.extent(1), 0);
     for (size_t i = 0; i < states_T.extent(0); ++i) {
-      std::complex<T> c_m_i = c * m[i];
-      size_t k_i = k[i];
-      std::memcpy(&new_states_T(i, 0), &states_T(k_i, 0),
-                  states_T.extent(1) * sizeof(std::complex<T>));
-
+      std::copy_n(&states_T(k[i], 0), states_T.extent(1), &new_states_T(i, 0));
+      const std::complex<T> c_m_i = c * m[i];
       for (size_t t = 0; t < states_T.extent(1); ++t) {
-        new_states_T(i, t) *= c_m_i; // * states_T(k_i, t);
+        new_states_T(i, t) *= c_m_i;
       }
     }
   }
@@ -309,14 +306,13 @@ struct PauliString {
   template <std::floating_point T>
   std::vector<std::vector<std::complex<T>>> get_dense_repr() const {
     //
-    std::vector<size_t> j, k;
+    std::vector<size_t> k;
     std::vector<std::complex<T>> m;
-    get_sparse_repr(j, k, m);
+    get_sparse_repr(k, m);
 
     // Convert to dense representation
-    size_t const dim = 1UL << paulis.size();
     std::vector<std::vector<std::complex<T>>> result(
-        dim, std::vector<std::complex<T>>(dim, 0));
+        dims(), std::vector<std::complex<T>>(dims(), 0));
 
     for (size_t i = 0; i < k.size(); ++i) {
       result[i][k[i]] = m[i];
@@ -329,6 +325,7 @@ struct PauliString {
 //
 // Helper
 //
+// TODO arrange following functions in a separate header __pauli_utils.hpp
 
 /**
  * @brief Get the nontrivial sets of pauli matrices given a weight.
@@ -352,7 +349,7 @@ std::vector<std::string> get_nontrivial_paulis(size_t const weight) {
         updated_set_of_nontrivial_paulis.push_back(str + pauli);
       }
     }
-    set_of_nontrivial_paulis = updated_set_of_nontrivial_paulis;
+    set_of_nontrivial_paulis = std::move(updated_set_of_nontrivial_paulis);
   }
   return set_of_nontrivial_paulis;
 }

fast_pauli/cpp/src/CMakeLists.txt

@@ -0,0 +1,12 @@
+# Grab all *.cpp files in the directory
+file(GLOB FAST_PAULI_BINDINGS CONFIGURE_DEPENDS "*.cpp")
+
+find_package(pybind11 REQUIRED)
+
+# TODO use different set of parameters depending on CMAKE_BUILD_TYPE
+pybind11_add_module(_fast_pauli SHARED ${FAST_PAULI_BINDINGS})
+target_link_libraries(_fast_pauli PUBLIC fast_pauli)
+# TODO recheck from where it is pulling in numpy headers
+
+install(TARGETS _fast_pauli
+        LIBRARY DESTINATION "${PROJECT_SOURCE_DIR}/fast_pauli")