Add generate_samples() to lightning

pennylane_lightning/_version.py

@@ -15,5 +15,4 @@
 """Version information.
    Version number (major.minor.patch[-label])
 """
-
-__version__ = "0.22.0-dev15"
+__version__ = "0.22.0-dev16"

pennylane_lightning/lightning_qubit.py

@@ -564,6 +564,33 @@ def probability(self, wires=None, shot_range=None, bin_size=None):
 
         return M.probs(device_wires)
 
+    def generate_samples(self):
+        """Generate samples
+
+        Returns:
+            Array of samples in binary
+        """
+
+        # To support np.complex64 based on the type of self._state
+        dtype = self._state.dtype
+        if dtype == np.complex64:
+            use_csingle = True
+        elif dtype == np.complex128:
+            use_csingle = False
+        else:
+            raise TypeError(f"Unsupported complex Type: {dtype}")
+
+        # Initialization of state
+        ket = np.ravel(self._state)
+
+        if use_csingle:
+            ket = ket.astype(np.complex64)
+
+        state_vector = StateVectorC64(ket) if use_csingle else StateVectorC128(ket)
+        M = MeasuresC64(state_vector) if use_csingle else MeasuresC128(state_vector)
+
+        return M.generate_samples(self.shots)
+
     def expval(self, observable, shot_range=None, bin_size=None):
         """Expectation value of the supplied observable.
 

pennylane_lightning/src/bindings/Bindings.cpp

@@ -324,6 +324,29 @@ void lightning_class_bindings(py::module &m) {
                 const std::vector<size_t> &wires) {
                  return M.expval(operation, wires);
              })
+        .def("generate_samples",
+             [](Measures<PrecisionT> &M, size_t num_shots) {
+                 auto &&result = M.generate_samples(num_shots);
+                 // int num_samples = num_shots;
+                 // std::vector<int> result(num_samples*wires.size(),0);
+
+                 int num_wires = result.size() / num_shots;
+                 ssize_t ndim = 2;
+                 std::vector<ssize_t> shape = {static_cast<ssize_t>(num_shots),
+					       static_cast<ssize_t>(num_wires)};
+                 std::vector<ssize_t> strides = {static_cast<ssize_t>(sizeof(int) * num_wires),
+						 static_cast<ssize_t>(sizeof(int))};
+
+                 // return 2-D NumPy array
+                 return py::array(py::buffer_info(
+                     result.data(), /* data as contiguous array  */
+                     sizeof(int),   /* size of one scalar        */
+                     py::format_descriptor<int>::format(), /* data type */
+                     ndim,   /* number of dimensions      */
+                     shape,  /* shape of the matrix       */
+                     strides /* strides for each axis     */
+                     ));
+             })
         .def("var", [](Measures<PrecisionT> &M, const std::string &operation,
                        const std::vector<size_t> &wires) {
             return M.var(operation, wires);

pennylane_lightning/src/gates/GateUtil.cpp

@@ -59,8 +59,8 @@ template <class OperatorImplementation> struct ImplementedGenerators {
 };
 
 template <class TypeList, class ValueType, template <class> class ValueClass>
-auto ValueForKernelHelper(
-    [[maybe_unused]] Pennylane::Gates::KernelType kernel) {
+auto ValueForKernelHelper([
+    [maybe_unused]] Pennylane::Gates::KernelType kernel) {
     if constexpr (std::is_same_v<TypeList, void>) {
         return std::vector<ValueType>{};
     } else {

pennylane_lightning/src/simulator/Measures.hpp

@@ -23,13 +23,17 @@
 #include <algorithm>
 #include <complex>
 #include <cstdio>
+#include <random>
+#include <stack>
+#include <unordered_map>
 #include <vector>
 
 #include "LinearAlgebra.hpp"
 #include "StateVectorManaged.hpp"
 #include "StateVectorRaw.hpp"
 
 namespace Pennylane {
+
 /**
  * @brief Observable's Measurement Class.
  *
@@ -256,5 +260,92 @@ class Measures {
 
         return expected_value_list;
     };
+
+    /**
+     * @brief Generate samples using the alias method.
+     * Reference: https://en.wikipedia.org/wiki/Alias_method
+     *
+     * @param num_samples The number of samples to generate.
+     * @return 1-D vector of samples in binary, each sample is
+     * separated by a stride equal to the number of qubits.
+     */
+    std::vector<size_t> generate_samples(size_t num_samples) {
+        size_t num_qubits = original_statevector.getNumQubits();
+        auto &&probabilities = probs();
+
+        std::vector<size_t> samples(num_samples * num_qubits, 0);
+        std::mt19937 generator(std::random_device{}());
+        std::uniform_real_distribution<fp_t> distribution(0.0, 1.0);
+        std::unordered_map<size_t, size_t> cache;
+
+        size_t N = probabilities.size();
+        std::vector<double> bucket(N);
+        std::vector<size_t> bucket_partner(N);
+        std::stack<size_t> overfull_bucket_ids;
+        std::stack<size_t> underfull_bucket_ids;
+
+        for (size_t i = 0; i < N; i++) {
+            bucket[i] = N * probabilities[i];
+            bucket_partner[i] = i;
+            if (bucket[i] > 1.0) {
+                overfull_bucket_ids.push(i);
+            }
+            if (bucket[i] < 1.0) {
+                underfull_bucket_ids.push(i);
+            }
+        }
+
+        // Run alias algorithm
+        while (underfull_bucket_ids.size() != 0 and
+               overfull_bucket_ids.size() != 0) {
+            // get an overfull bucket
+            auto i = overfull_bucket_ids.top();
+
+            // get an underfull bucket
+            auto j = underfull_bucket_ids.top();
+            underfull_bucket_ids.pop();
+
+            // underfull bucket is partned with an overfull bucket
+            bucket_partner[j] = i;
+            bucket[i] = bucket[i] + bucket[j] - 1;
+
+            // if overfull bucket is now underfull
+            // put in underfull stack
+            if (bucket[i] < 1) {
+                overfull_bucket_ids.pop();
+                underfull_bucket_ids.push(i);
+            }
+
+            // if overfull bucket is full -> remove
+            else if (bucket[i] == 1.0) {
+                overfull_bucket_ids.pop();
+            }
+        }
+
+        // Pick samples
+        for (size_t i = 0; i < num_samples; i++) {
+            fp_t pct = distribution(generator) * N;
+            size_t idx = static_cast<size_t>(pct);
+            if (pct - idx > bucket[idx])
+                idx = bucket_partner[idx];
+
+            // If cached, retrieve sample 1from cache
+            if (cache.count(idx) != 0) {
+                size_t cache_id = cache[idx];
+                auto it_temp = samples.begin() + cache_id * num_qubits;
+                std::copy(it_temp, it_temp + num_qubits,
+                          samples.begin() + i * num_qubits);
+            }
+            // If not cached, compute
+            else {
+                for (size_t j = 0; j < num_qubits; j++)
+                    samples[i * num_qubits + (num_qubits - 1 - j)] =
+                        (idx >> j) & 1;
+                cache[idx] = i;
+            }
+        }
+
+        return samples;
+    }
 }; // class Measures
 } // namespace Pennylane

pennylane_lightning/src/simulator/StateVectorBase.hpp

@@ -46,7 +46,7 @@
 #define PENNYLANE_STATEVECTOR_DEFINE_GATE(GATE_NAME)                           \
     template <Gates::KernelType kernel, typename... Ts>                        \
     inline void apply##GATE_NAME##_(const std::vector<size_t> &wires,          \
-                                    bool inverse, Ts &&...args) {              \
+                                    bool inverse, Ts &&... args) {             \
         auto *arr = getData();                                                 \
         static_assert(Gates::static_lookup<Gates::GateOperation::GATE_NAME>(   \
                           Gates::Constant::gate_num_params) == sizeof...(Ts),  \
@@ -63,7 +63,7 @@
 #define PENNYLANE_STATEVECTOR_DEFINE_DEFAULT_GATE(GATE_NAME)                   \
     template <typename... Ts>                                                  \
     inline void apply##GATE_NAME(const std::vector<size_t> &wires,             \
-                                 bool inverse, Ts &&...args) {                 \
+                                 bool inverse, Ts &&... args) {                \
         constexpr auto kernel =                                                \
             Gates::static_lookup<Gates::GateOperation::GATE_NAME>(             \
                 Gates::Constant::default_kernel_for_gates);                    \

pennylane_lightning/src/tests/Test_GateImplementations_Nonparam.cpp

@@ -37,9 +37,8 @@ using std::vector;
     template <typename PrecisionT, class GateImplementation>                   \
     struct Apply##GATE_NAME##IsDefined<                                        \
         PrecisionT, GateImplementation,                                        \
-        std::enable_if_t<std::is_pointer_v<                                    \
-            decltype(&GateImplementation::template apply##GATE_NAME<           \
-                     PrecisionT>)>>> {                                         \
+        std::enable_if_t<std::is_pointer_v<decltype(                           \
+            &GateImplementation::template apply##GATE_NAME<PrecisionT>)>>> {   \
         constexpr static bool value = true;                                    \
     };                                                                         \
     template <typename PrecisionT, typename TypeList>                          \

pennylane_lightning/src/tests/Test_GateImplementations_Param.cpp

@@ -30,9 +30,9 @@ using namespace Pennylane;
     template <typename PrecisionT, typename ParamT, class GateImplementation>  \
     struct Apply##GATE_NAME##IsDefined<                                        \
         PrecisionT, ParamT, GateImplementation,                                \
-        std::enable_if_t<std::is_pointer_v<                                    \
-            decltype(&GateImplementation::template apply##GATE_NAME<           \
-                     PrecisionT, ParamT>)>>> {                                 \
+        std::enable_if_t<std::is_pointer_v<decltype(                           \
+            &GateImplementation::template apply##GATE_NAME<PrecisionT,         \
+                                                           ParamT>)>>> {       \
         constexpr static bool value = true;                                    \
     };                                                                         \
     template <typename PrecisionT, typename ParamT, typename TypeList>         \

pennylane_lightning/src/tests/Test_Measures.cpp

@@ -163,6 +163,55 @@ TEST_CASE("Expected Values", "[Measures]") {
     }
 }
 
+TEST_CASE("Sample", "[Measures]") {
+    constexpr int twos[] = {
+        1 << 0,  1 << 1,  1 << 2,  1 << 3,  1 << 4,  1 << 5,  1 << 6,  1 << 7,
+        1 << 8,  1 << 9,  1 << 10, 1 << 11, 1 << 12, 1 << 13, 1 << 14, 1 << 15,
+        1 << 16, 1 << 17, 1 << 18, 1 << 19, 1 << 20, 1 << 21, 1 << 22, 1 << 23,
+        1 << 24, 1 << 25, 1 << 26, 1 << 27, 1 << 28, 1 << 29, 1 << 30, 1 << 31};
+
+    // Defining the State Vector that will be measured.
+    StateVectorManaged<double> Measured_StateVector =
+        Initializing_StateVector();
+
+    // Initializing the measures class.
+    // It will attach to the StateVector, allowing measures to keep been taken.
+    Measures<double, StateVectorManaged<double>> Measurer(Measured_StateVector);
+    vector<double> expected_probabilities = {0.687573, 0.013842, 0.089279,
+                                             0.001797, 0.180036, 0.003624,
+                                             0.023377, 0.000471};
+
+    size_t num_qubits = 3;
+    size_t N = std::pow(2, num_qubits);
+    size_t num_samples = 10000;
+    auto &&samples = Measurer.generate_samples(num_samples);
+
+    std::vector<int> counts(N, 0);
+    std::vector<int> samples_decimal(num_samples, 0);
+
+    // convert samples to decimal and then bin them in counts
+    for (size_t i = 0; i < num_samples; i++) {
+        for (size_t j = 0; j < num_qubits; j++) {
+            if (samples[i * num_qubits + j] != 0) {
+                samples_decimal[i] += twos[(num_qubits - 1 - j)];
+            }
+        }
+        counts[samples_decimal[i]] += 1;
+    }
+
+    // compute estimated probabilities from histogram
+    std::vector<double> probabilities(counts.size());
+    for (size_t i = 0; i < counts.size(); i++) {
+        probabilities[i] = counts[i] / (double)num_samples;
+    }
+
+    // compare estimated probabilities to real probabilities
+    SECTION("No wires provided:") {
+        REQUIRE_THAT(probabilities,
+                     Catch::Approx(expected_probabilities).margin(1e-3));
+    }
+}
+
 TEST_CASE("Variances", "[Measures]") {
     // Defining the State Vector that will be measured.
     StateVectorManaged<double> Measured_StateVector =

pennylane_lightning/src/tests/Test_OpToMemberFuncPtr.cpp

@@ -205,8 +205,8 @@ template <typename PrecisionT, typename ParamT, size_t num_params,
           size_t tuple_idx>
 constexpr auto gateOpFuncPtrPairsWithNumParamsIter() {
     if constexpr (tuple_idx <
-                  std::tuple_size_v<
-                      decltype(gate_op_func_ptr_pairs<PrecisionT, ParamT>)>) {
+                  std::tuple_size_v<decltype(
+                      gate_op_func_ptr_pairs<PrecisionT, ParamT>)>) {
         constexpr auto elt =
             std::get<tuple_idx>(gate_op_func_ptr_pairs<PrecisionT, ParamT>);
         if constexpr (static_lookup<elt.first>(Constant::gate_num_params) ==

tests/test_adjoint_jacobian.py

@@ -207,7 +207,7 @@ def test_pauli_rotation_gradient(self, G, theta, tol, dev, C):
     def test_Rot_gradient(self, theta, tol, dev, C):
         """Tests that the device gradient of an arbitrary Euler-angle-parameterized gate is
         correct."""
-        params = np.array([theta, theta**3, np.sqrt(2) * theta])
+        params = np.array([theta, theta ** 3, np.sqrt(2) * theta])
 
         dev._state = dev._state.astype(C)
 
@@ -462,8 +462,8 @@ def circuit(x, y, z):
 
     thetas = np.linspace(-2 * np.pi, 2 * np.pi, 8)
 
-    @pytest.mark.parametrize("reused_p", thetas**3 / 19)
-    @pytest.mark.parametrize("other_p", thetas**2 / 1)
+    @pytest.mark.parametrize("reused_p", thetas ** 3 / 19)
+    @pytest.mark.parametrize("other_p", thetas ** 2 / 1)
     def test_fanout_multiple_params(self, reused_p, other_p, tol, mocker, dev):
         """Tests that the correct gradient is computed for qnodes which
         use the same parameter in multiple gates."""
@@ -624,7 +624,7 @@ def f(params1, params2):
 
 def circuit_ansatz(params, wires):
     """Circuit ansatz containing all the parametrized gates"""
-    qml.QubitStateVector(unitary_group.rvs(2**4, random_state=0)[0], wires=wires)
+    qml.QubitStateVector(unitary_group.rvs(2 ** 4, random_state=0)[0], wires=wires)
     qml.RX(params[0], wires=wires[0])
     qml.RY(params[1], wires=wires[1])
     qml.RX(params[2], wires=wires[2]).inv()

tests/test_apply.py

@@ -626,7 +626,7 @@ def test_sample_values(self, qubit_device_2_wires, tol):
 
         # s1 should only contain 1 and -1, which is guaranteed if
         # they square to 1
-        assert np.allclose(s1**2, 1, atol=tol, rtol=0)
+        assert np.allclose(s1 ** 2, 1, atol=tol, rtol=0)
 
 
 class TestLightningQubitIntegration:
@@ -747,7 +747,7 @@ def circuit(x):
     def test_nonzero_shots(self, tol):
         """Test that the default qubit plugin provides correct result for high shot number"""
 
-        shots = 10**4
+        shots = 10 ** 4
         dev = qml.device("lightning.qubit", wires=1, shots=shots)
 
         p = 0.543
@@ -1330,13 +1330,13 @@ def test_paulix_pauliy(self, theta, phi, varphi, shots, monkeypatch, tol):
         p = dev.probability(wires=obs.wires)
 
         # s1 should only contain 1 and -1
-        assert np.allclose(s1**2, 1, atol=tolerance, rtol=0)
+        assert np.allclose(s1 ** 2, 1, atol=tolerance, rtol=0)
 
         mean = s1 @ p
         expected = np.sin(theta) * np.sin(phi) * np.sin(varphi)
         assert np.allclose(mean, expected, atol=tolerance, rtol=0)
 
-        var = (s1**2) @ p - (s1 @ p).real ** 2
+        var = (s1 ** 2) @ p - (s1 @ p).real ** 2
         expected = (
             8 * np.sin(theta) ** 2 * np.cos(2 * varphi) * np.sin(phi) ** 2
             - np.cos(2 * (theta - phi))
@@ -1372,13 +1372,13 @@ def test_pauliz_hadamard(
         p = dev.marginal_prob(dev.probability(), wires=obs.wires)
 
         # s1 should only contain 1 and -1
-        assert np.allclose(s1**2, 1, atol=tol, rtol=0)
+        assert np.allclose(s1 ** 2, 1, atol=tol, rtol=0)
 
         mean = s1 @ p
         expected = -(np.cos(varphi) * np.sin(phi) + np.sin(varphi) * np.cos(theta)) / np.sqrt(2)
         assert np.allclose(mean, expected, atol=tol, rtol=0)
 
-        var = (s1**2) @ p - (s1 @ p).real ** 2
+        var = (s1 ** 2) @ p - (s1 @ p).real ** 2
         expected = (
             3
             + np.cos(2 * phi) * np.cos(varphi) ** 2
@@ -1418,13 +1418,13 @@ def test_qubitunitary_rotation_hadamard(
         p = dev.marginal_prob(dev.probability(), wires=obs.wires)
 
         # s1 should only contain 1 and -1
-        assert np.allclose(s1**2, 1, atol=tol, rtol=0)
+        assert np.allclose(s1 ** 2, 1, atol=tol, rtol=0)
 
         mean = s1 @ p
         expected = -(np.cos(varphi) * np.sin(phi) + np.sin(varphi) * np.cos(theta)) / np.sqrt(2)
         assert np.allclose(mean, expected, atol=tol, rtol=0)
 
-        var = (s1**2) @ p - (s1 @ p).real ** 2
+        var = (s1 ** 2) @ p - (s1 @ p).real ** 2
         expected = (
             3
             + np.cos(2 * phi) * np.cos(varphi) ** 2