[feature/SPO_apply] Rounding out python tests

fast_pauli/cpp/examples/03_summed_pauli_op.cpp

@@ -67,4 +67,4 @@ int main()
     summed_op.apply_weighted(std::execution::par, new_states, states, weights);
 
     return 0;
-}
+}

fast_pauli/cpp/examples/04_get_sparse_repr.cpp

@@ -14,4 +14,4 @@ int main()
     std::tie(k, m) = get_sparse_repr<double>(ps.paulis);
 
     return 0;
-}
+}

fast_pauli/cpp/include/__summed_pauli_op.hpp

@@ -542,6 +542,8 @@ template <std::floating_point T> struct SummedPauliOp
 #pragma omp parallel for schedule(static)
         for (size_t k = 0; k < n_operators(); ++k)
         {
+            // This is a lot of wasted work, but in some quick tests with n_ops = 1000, n_qubits = 10, weight = 2 this
+            // strategy was slightly faster than have the loop over i being the outer-most loop.
             for (size_t i = 0; i < pauli_strings.size(); ++i)
             {
                 PauliString const &ps = pauli_strings[i];

tests/fast_pauli/test_summed_pauli_op.py

@@ -30,36 +30,32 @@
     "n_states,n_operators,n_qubits",
     [(s, o, q) for s in [1, 10, 1000] for o in [1, 10, 100] for q in [1, 2, 6]],
 )
-def test_expectation_values(
+def test_apply(
     summed_pauli_op: type[fp.SummedPauliOp],
     n_states: int,
     n_operators: int,
     n_qubits: int,
 ) -> None:
-    """Test expectation value calculation."""
-    # TODO This will break on python version
+    """Test applying the summed pauli operator method."""
     pauli_strings = fp.helpers.calculate_pauli_strings_max_weight(n_qubits, 2)
-    pauli_strings_str = [str(s) for s in pauli_strings]
     n_strings = len(pauli_strings)
 
     coeffs_2d = np.random.rand(n_strings, n_operators).astype(np.complex128)
     psi = np.random.rand(2**n_qubits, n_states).astype(np.complex128)
 
     op = summed_pauli_op(pauli_strings, coeffs_2d)
 
-    # The expectation values we want to test
-    expectation_vals = op.expectation_value(psi)
-
-    # The "trusted" expectation_values
-    expectation_vals_naive = np.zeros((n_operators, n_states), dtype=np.complex128)
+    # The new_states we want to check
+    new_states = op.apply(psi)
 
-    # Calculate using brute force
+    # Trusted new_states
+    new_states_naive = np.zeros((2**n_qubits, n_states), dtype=np.complex128)
     for k in range(n_operators):
-        A_k = pp.helpers.naive_pauli_operator(coeffs_2d[:, k], pauli_strings_str)
-        expectation_vals_naive[k] = np.einsum("it,ij,jt->t", psi.conj(), A_k, psi)
+        A_k = fp.PauliOp(coeffs_2d[:, k].copy(), pauli_strings)
+        new_states_naive += A_k.apply(psi)
 
     # Check
-    np.testing.assert_allclose(expectation_vals, expectation_vals_naive, atol=1e-13)
+    np.testing.assert_allclose(new_states, new_states_naive, atol=1e-13)
 
 
 @pytest.mark.parametrize(
@@ -69,7 +65,7 @@ def test_expectation_values(
     "n_states,n_operators,n_qubits",
     [(s, o, q) for s in [1, 10, 1000] for o in [1, 10, 100] for q in [1, 2, 6]],
 )
-def test_apply(
+def test_apply_weighted(
     summed_pauli_op: type[fp.SummedPauliOp],
     n_states: int,
     n_operators: int,
@@ -81,22 +77,62 @@ def test_apply(
 
     coeffs_2d = np.random.rand(n_strings, n_operators).astype(np.complex128)
     psi = np.random.rand(2**n_qubits, n_states).astype(np.complex128)
+    data_weights = np.random.rand(n_operators, n_states).astype(np.float64)
 
     op = summed_pauli_op(pauli_strings, coeffs_2d)
 
     # The new_states we want to check
-    new_states = op.apply(psi)
+    new_states = op.apply_weighted(psi, data_weights)
 
     # Trusted new_states
     new_states_naive = np.zeros((2**n_qubits, n_states), dtype=np.complex128)
     for k in range(n_operators):
         A_k = fp.PauliOp(coeffs_2d[:, k].copy(), pauli_strings)
-        new_states_naive += A_k.apply(psi)
+        new_states_naive += A_k.apply(psi) * data_weights[k]
 
     # Check
     np.testing.assert_allclose(new_states, new_states_naive, atol=1e-13)
 
 
+@pytest.mark.parametrize(
+    "summed_pauli_op", [fp.SummedPauliOp], ids=resolve_parameter_repr
+)
+@pytest.mark.parametrize(
+    "n_states,n_operators,n_qubits",
+    [(s, o, q) for s in [1, 10, 1000] for o in [1, 10, 100] for q in [1, 2, 6]],
+)
+def test_expectation_values(
+    summed_pauli_op: type[fp.SummedPauliOp],
+    n_states: int,
+    n_operators: int,
+    n_qubits: int,
+) -> None:
+    """Test expectation value calculation."""
+    # TODO This will break on python version
+    pauli_strings = fp.helpers.calculate_pauli_strings_max_weight(n_qubits, 2)
+    pauli_strings_str = [str(s) for s in pauli_strings]
+    n_strings = len(pauli_strings)
+
+    coeffs_2d = np.random.rand(n_strings, n_operators).astype(np.complex128)
+    psi = np.random.rand(2**n_qubits, n_states).astype(np.complex128)
+
+    op = summed_pauli_op(pauli_strings, coeffs_2d)
+
+    # The expectation values we want to test
+    expectation_vals = op.expectation_value(psi)
+
+    # The "trusted" expectation_values
+    expectation_vals_naive = np.zeros((n_operators, n_states), dtype=np.complex128)
+
+    # Calculate using brute force
+    for k in range(n_operators):
+        A_k = pp.helpers.naive_pauli_operator(coeffs_2d[:, k], pauli_strings_str)
+        expectation_vals_naive[k] = np.einsum("it,ij,jt->t", psi.conj(), A_k, psi)
+
+    # Check
+    np.testing.assert_allclose(expectation_vals, expectation_vals_naive, atol=1e-13)
+
+
 @pytest.mark.parametrize(
     "summed_pauli_op", [fp.SummedPauliOp], ids=resolve_parameter_repr
 )