Add generate_samples() to lightning

pennylane_lightning/_version.py

@@ -15,5 +15,4 @@
 """Version information.
    Version number (major.minor.patch[-label])
 """
-
-__version__ = "0.23.0-dev5"
+__version__ = "0.23.0-dev6"

pennylane_lightning/lightning_qubit.py

@@ -576,6 +576,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[int]: array of samples in binary representation with shape ``(dev.shots, dev.num_wires)``
+        """
+
+        # 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(len(self.wires), self.shots).astype(int)
+
     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,23 @@ 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_wires, size_t num_shots) {
+                 auto &&result = M.generate_samples(num_shots);
+                 const size_t ndim = 2;
+                 const std::vector<size_t> shape{num_shots, num_wires};
+                 constexpr auto sz = sizeof(size_t);
+                 const std::vector<size_t> strides{sz * num_wires, sz};
+                 // return 2-D NumPy array
+                 return py::array(py::buffer_info(
+                     result.data(), /* data as contiguous array  */
+                     sz,            /* size of one scalar        */
+                     py::format_descriptor<size_t>::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/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) {
+        const 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;
+
+        const 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.empty() && !overfull_bucket_ids.empty()) {
+            // get an overfull bucket
+            size_t i = overfull_bucket_ids.top();
+
+            // get an underfull bucket
+            size_t 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 = pct;
+            if (pct - idx > bucket[idx]) {
+                idx = bucket_partner[idx];
+            }
+            // If cached, retrieve sample from 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/tests/Test_Measures.cpp

@@ -18,20 +18,21 @@ using std::string;
 using std::vector;
 }; // namespace
 
-StateVectorManaged<double> Initializing_StateVector() {
+template <typename T = double>
+StateVectorManaged<T> Initializing_StateVector() {
     // Defining a StateVector in a non-trivial configuration:
     size_t num_qubits = 3;
     size_t data_size = std::pow(2, num_qubits);
 
-    std::vector<std::complex<double>> arr(data_size, 0);
+    std::vector<std::complex<T>> arr(data_size, 0);
     arr[0] = 1;
-    StateVectorManaged<double> Measured_StateVector(arr.data(), data_size);
+    StateVectorManaged<T> Measured_StateVector(arr.data(), data_size);
 
     std::vector<size_t> wires;
 
-    double alpha = 0.7;
-    double beta = 0.5;
-    double gamma = 0.2;
+    T alpha = 0.7;
+    T beta = 0.5;
+    T gamma = 0.2;
     Measured_StateVector.applyOperations(
         {"RX", "RY", "RX", "RY", "RX", "RY"}, {{0}, {0}, {1}, {1}, {2}, {2}},
         {false, false, false, false, false, false},
@@ -163,39 +164,93 @@ TEST_CASE("Expected Values", "[Measures]") {
     }
 }
 
-TEST_CASE("Variances", "[Measures]") {
+TEMPLATE_TEST_CASE("Sample", "[Measures]", float, double) {
+    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();
+    StateVectorManaged<TestType> Measured_StateVector =
+        Initializing_StateVector<TestType>();
 
     // Initializing the measures class.
     // It will attach to the StateVector, allowing measures to keep been taken.
-    Measures<double, StateVectorManaged<double>> Measurer(Measured_StateVector);
+    Measures<TestType, StateVectorManaged<TestType>> Measurer(
+        Measured_StateVector);
+    vector<TestType> 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 = 100000;
+    auto &&samples = Measurer.generate_samples(num_samples);
+
+    std::vector<size_t> counts(N, 0);
+    std::vector<size_t> 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<TestType> probabilities(counts.size());
+    for (size_t i = 0; i < counts.size(); i++) {
+        probabilities[i] = counts[i] / (TestType)num_samples;
+    }
+
+    // compare estimated probabilities to real probabilities
+    SECTION("No wires provided:") {
+        REQUIRE_THAT(probabilities,
+                     Catch::Approx(expected_probabilities).margin(.05));
+    }
+}
+
+TEMPLATE_TEST_CASE("Variances", "[Measures]", float, double) {
+    // Defining the State Vector that will be measured.
+    StateVectorManaged<TestType> Measured_StateVector =
+        Initializing_StateVector<TestType>();
+
+    // Initializing the measures class.
+    // It will attach to the StateVector, allowing measures to keep been taken.
+    Measures<TestType, StateVectorManaged<TestType>> Measurer(
+        Measured_StateVector);
 
     SECTION("Testing single operation defined by a matrix:") {
-        vector<std::complex<double>> PauliX = {0, 1, 1, 0};
+        vector<std::complex<TestType>> PauliX = {0, 1, 1, 0};
         vector<size_t> wires_single = {0};
-        double variance = Measurer.var(PauliX, wires_single);
-        double variances_ref = 0.757222;
+        TestType variance = Measurer.var(PauliX, wires_single);
+        TestType variances_ref = 0.757222;
         REQUIRE(variance == Approx(variances_ref).margin(1e-6));
     }
 
     SECTION("Testing single operation defined by its name:") {
         vector<size_t> wires_single = {0};
-        double variance = Measurer.var("PauliX", wires_single);
-        double variances_ref = 0.757222;
+        TestType variance = Measurer.var("PauliX", wires_single);
+        TestType variances_ref = 0.757222;
         REQUIRE(variance == Approx(variances_ref).margin(1e-6));
     }
 
     SECTION("Testing list of operators defined by a matrix:") {
-        vector<std::complex<double>> PauliX = {0, 1, 1, 0};
-        vector<std::complex<double>> PauliY = {0, {0, -1}, {0, 1}, 0};
-        vector<std::complex<double>> PauliZ = {1, 0, 0, -1};
-
-        vector<double> variances;
-        vector<double> variances_ref;
+        vector<std::complex<TestType>> PauliX = {
+            {0, 0}, {1, 0}, {1, 0}, {0, 0}};
+        vector<std::complex<TestType>> PauliY = {
+            {0, 0}, {0, -1}, {0, 1}, {0, 0}};
+        vector<std::complex<TestType>> PauliZ = {
+            {1, 0}, {0, 0}, {0, 0}, {-1, 0}};
+
+        vector<TestType> variances;
+        vector<TestType> variances_ref;
         vector<vector<size_t>> wires_list = {{0}, {1}, {2}};
-        vector<vector<std::complex<double>>> operations_list;
+        vector<vector<std::complex<TestType>>> operations_list;
 
         operations_list = {PauliX, PauliX, PauliX};
         variances = Measurer.var(operations_list, wires_list);
@@ -214,8 +269,8 @@ TEST_CASE("Variances", "[Measures]") {
     }
 
     SECTION("Testing list of operators defined by its name:") {
-        vector<double> variances;
-        vector<double> variances_ref;
+        vector<TestType> variances;
+        vector<TestType> variances_ref;
         vector<vector<size_t>> wires_list = {{0}, {1}, {2}};
         vector<string> operations_list;