Merge pull request #80 from qognitive/feature/SPO_apply

[Feature] Adding `apply` to `SummedPauliOp`

docs/python_api.rst

@@ -27,3 +27,10 @@ PauliOp
 .. autoclass:: fast_pauli.PauliOp
    :members:
    :special-members:
+
+SummedPauliOp
+-------------
+
+.. autoclass:: fast_pauli.SummedPauliOp
+   :members:
+   :special-members:

fast_pauli/cpp/examples/03_summed_pauli_op.cpp

@@ -13,38 +13,35 @@
  */
 
 #include <algorithm>
+#include <chrono>
 #include <experimental/mdspan>
 #include <random>
 
 #include "__pauli.hpp"
 #include "__pauli_string.hpp"
 #include "fast_pauli.hpp"
 
-using namespace fast_pauli;
+namespace fp = fast_pauli;
 
 int main()
 {
     //
     // User settings
     //
-    size_t const n_operators = 10000;
-    //   size_t const n_paulis_per_operator = 631;
-    size_t const n_qubits = 14;
+    size_t const n_operators = 1000;
+    size_t const n_qubits = 12;
     size_t const weight = 2;
     size_t const n_states = 1000;
-    //   std::string fake_pauli_string = "XYZXYZXYZXYZ";
-    using fp_type = double;
+    using fp_type = float;
 
     //
     // Setup the summed pauli operator
     //
-    //   std::vector<PauliString> pauli_strings(n_paulis_per_operator,
-    //                                          PauliString(fake_pauli_string));
-    std::vector<PauliString> pauli_strings = calculate_pauli_strings_max_weight(n_qubits, weight);
+    std::vector<fp::PauliString> pauli_strings = fp::calculate_pauli_strings_max_weight(n_qubits, weight);
 
     size_t const n_paulis_per_operator = pauli_strings.size();
     std::vector<std::complex<fp_type>> coeff_raw(n_paulis_per_operator * n_operators, 1);
-    SummedPauliOp<fp_type> summed_op{pauli_strings, coeff_raw};
+    fp::SummedPauliOp<fp_type> summed_op{pauli_strings, coeff_raw};
 
     //
     // Setup states
@@ -67,8 +64,17 @@ int main()
     //
     // Apply the states
     //
-    //   summed_op.apply(new_states, states, weights);
-    summed_op.apply_weighted(new_states, states, weights);
+    // auto start_seq = std::chrono::high_resolution_clock::now();
+    // summed_op.apply_weighted(std::execution::seq, new_states, states, weights);
+    // auto end_seq = std::chrono::high_resolution_clock::now();
+    // fmt::println("Time taken for sequential execution: {} seconds",
+    //              std::chrono::duration_cast<std::chrono::seconds>(end_seq - start_seq).count());
+
+    auto start_par = std::chrono::high_resolution_clock::now();
+    summed_op.apply_weighted(std::execution::par, new_states, states, weights);
+    auto end_par = std::chrono::high_resolution_clock::now();
+    fmt::println("Time taken for parallel execution:   {} seconds",
+                 std::chrono::duration_cast<std::chrono::seconds>(end_par - start_par).count());
 
     return 0;
-}
+}

fast_pauli/cpp/examples/04_get_sparse_repr.cpp

@@ -0,0 +1,17 @@
+#include <string>
+
+#include "fast_pauli.hpp"
+
+namespace fp = fast_pauli;
+
+int main()
+{
+    size_t const n_qubits = 24;
+    std::string pauli_string(n_qubits, 'X');
+    fp::PauliString ps(pauli_string);
+    std::vector<size_t> k;
+    std::vector<std::complex<double>> m;
+    std::tie(k, m) = get_sparse_repr<double>(ps.paulis);
+
+    return 0;
+}

fast_pauli/cpp/examples/05_summed_pauli_op_sq.cpp

@@ -0,0 +1,56 @@
+/**
+ * This code is part of Fast Pauli.
+ *
+ * (C) Copyright Qognitive Inc 2024.
+ *
+ * This code is licensed under the BSD 2-Clause License. You may
+ * obtain a copy of this license in the LICENSE.txt file in the root directory
+ * of this source tree.
+ *
+ * Any modifications or derivative works of this code must retain this
+ * copyright notice, and modified files need to carry a notice indicating
+ * that they have been altered from the originals.
+ */
+
+#include <algorithm>
+#include <chrono>
+#include <experimental/mdspan>
+#include <random>
+
+#include "__pauli.hpp"
+#include "__pauli_string.hpp"
+#include "fast_pauli.hpp"
+
+namespace fp = fast_pauli;
+
+int main()
+{
+    //
+    // User settings
+    //
+    size_t const n_operators = 1000;
+    size_t const n_qubits = 12;
+    size_t const weight = 2;
+    using fp_type = float;
+
+    //
+    // Setup the summed pauli operator
+    //
+    std::vector<fp::PauliString> pauli_strings = fp::calculate_pauli_strings_max_weight(n_qubits, weight);
+
+    size_t const n_paulis_per_operator = pauli_strings.size();
+    std::vector<std::complex<fp_type>> coeff_raw(n_paulis_per_operator * n_operators, 1);
+    fp::SummedPauliOp<fp_type> summed_op{pauli_strings, coeff_raw};
+
+    //
+    // Apply the states
+    //
+
+    auto start_par = std::chrono::high_resolution_clock::now();
+    auto sq_op = summed_op.square();
+    auto end_par = std::chrono::high_resolution_clock::now();
+    fmt::println("Time taken for parallel execution:   {} seconds",
+                 std::chrono::duration_cast<std::chrono::seconds>(end_par - start_par).count());
+
+    return 0;
+}

fast_pauli/cpp/include/__pauli_helpers.hpp

@@ -112,9 +112,6 @@ std::vector<PauliString> calcutate_pauli_strings(size_t const n_qubits, size_t c
     size_t n_pauli_strings = nontrivial_paulis.size() * idx_combos.size();
     std::vector<PauliString> result(n_pauli_strings);
 
-    fmt::println("n_qubits = {}  weight = {}  n_nontrivial_paulis = {}  n_combos = {}", n_qubits, weight,
-                 nontrivial_paulis.size(), idx_combos.size());
-
     // Iterate through all the nontrivial paulis and all the combinations
     for (size_t i = 0; i < nontrivial_paulis.size(); ++i)
     {
@@ -152,8 +149,6 @@ std::vector<PauliString> calculate_pauli_strings_max_weight(size_t n_qubits, siz
         auto ps = calcutate_pauli_strings(n_qubits, i);
         result.insert(result.end(), ps.begin(), ps.end());
     }
-
-    fmt::println("n_qubits = {}  weight = {}  n_pauli_strings = {}", n_qubits, weight, result.size());
     return result;
 }
 

fast_pauli/cpp/include/__pauli_string.hpp

@@ -64,7 +64,7 @@ std::tuple<std::vector<size_t>, std::vector<std::complex<T>>> get_sparse_repr(st
 
     // Helper function that let's us know if a pauli matrix has diagonal (or
     // conversely off-diagonal) elements
-    auto diag = [](Pauli const &p) {
+    auto diag = [](Pauli const &p) -> size_t {
         if (p.code == 0 || p.code == 3)
         {
             return 0UL;
@@ -98,18 +98,13 @@ std::tuple<std::vector<size_t>, std::vector<std::complex<T>>> get_sparse_repr(st
     }
     m[0] = initial_value();
 
-    // Populate the rest of the values in a recursive-like manner
     for (size_t l = 0; l < n; ++l)
     {
         Pauli const &po = ps[l];
 
-        T eps = 1.0;
-        if (po.code == 2 || po.code == 3)
-        {
-            eps = -1;
-        }
+        T const eps = (po.code == 2 || po.code == 3) ? -1.0 : 1.0;
 
-        T sign = diag(po) ? -1.0 : 1.0;
+        int sign = diag(po) ? -1 : 1;
 
         auto const lower_bound = 1UL << l;
         for (size_t li = lower_bound; li < (lower_bound << 1); li++)