[sample] Detect measurement operations, disallow on hardware targets, warning on simulators

docs/sphinx/applications/cpp/amplitude_estimation.cpp

@@ -107,8 +107,6 @@ struct run_circuit {
       S_0{}(q);
       statePrep_A{}(q, bmax);
     }
-    // Measure the last auxiliary qubit
-    mz(last_qubit);
   }
 };
 
@@ -124,7 +122,8 @@ int main() {
 
   for (size_t i = 0; i < schedule.size(); i++) {
     auto counts = cudaq::sample(run_circuit{}, n, schedule[i], bmax);
-    hits[i] = counts.count("1");
+    auto last_qubit_counts = counts.get_marginal({n});
+    hits[i] = last_qubit_counts.count("1");
   }
 
   // Print the number of hits for the good state for each circuit

docs/sphinx/applications/cpp/bernstein_vazirani.cpp

@@ -12,6 +12,7 @@
 
 #include <cudaq.h>
 #include <iostream>
+#include <numeric>
 #include <random>
 #include <vector>
 
@@ -57,7 +58,6 @@ struct bernstein_vazirani {
 
     oracle<nrOfBits>{}(bitvector, qs, aux);
     h(qs);
-    mz(qs);
   }
 };
 
@@ -80,12 +80,15 @@ int main(int argc, char *argv[]) {
   auto bitvector = random_bits<nr_qubits>(seed);
   auto kernel = bernstein_vazirani<nr_qubits>{};
   auto counts = cudaq::sample(nr_shots, kernel, bitvector);
+  std::vector<std::size_t> indices(nr_qubits);
+  std::iota(indices.begin(), indices.end(), 0);
+  auto qs_counts = counts.get_marginal(indices);
 
   if (!cudaq::mpi::is_initialized() || cudaq::mpi::rank() == 0) {
     printf("Encoded bitstring:  %s\n", asString(bitvector).c_str());
-    printf("Measured bitstring: %s\n\n", counts.most_probable().c_str());
+    printf("Measured bitstring: %s\n\n", qs_counts.most_probable().c_str());
 
-    for (auto &[bits, count] : counts) {
+    for (auto &[bits, count] : qs_counts) {
       printf("observed %s with %.0f%% probability\n", bits.data(),
              100.0 * count / nr_shots);
     }

docs/sphinx/applications/cpp/grover.cpp

@@ -32,7 +32,6 @@ struct run_grover {
       oracle(qs);
       reflect_about_uniform(qs);
     }
-    mz(qs);
   }
 };
 

docs/sphinx/applications/cpp/iterative_qpe.cpp

@@ -12,10 +12,11 @@
 // nvq++ iterative_qpe.cpp -o qpe.x && ./qpe.x
 // ```
 
+#include <algorithm>
 #include <cudaq.h>
 
 struct iqpe {
-  void operator()() __qpu__ {
+  std::vector<bool> operator()() __qpu__ {
     cudaq::qarray<2> q;
     h(q[0]);
     x(q[1]);
@@ -63,14 +64,22 @@ struct iqpe {
       rz(-M_PI_2, q[0]);
 
     h(q[0]);
-    mz(q[0]);
+    return {cr0, cr1, cr2, mz(q[0])};
   }
 };
 
 int main() {
-  auto counts = cudaq::sample(/*shots*/ 10, iqpe{});
-  counts.dump();
-
+  auto results = cudaq::run(/*shots*/ 10, iqpe{});
+  // Get the counts for cr0, cr1, cr2 and the final measurement
+  auto count_bit = [&](std::size_t idx) {
+    return std::count_if(results.begin(), results.end(),
+                         [idx](auto &r) { return r[idx]; });
+  };
+  printf("Iterative QPE Results:\n");
+  printf("cr0 : { 1:%zu }\n", count_bit(0));
+  printf("cr1 : { 1:%zu }\n", count_bit(1));
+  printf("cr2 : { 0:%zu }\n", 10 - count_bit(2));
+  printf("final: { 1:%zu }\n", count_bit(3));
   return 0;
 }
 // [End Documentation]

docs/sphinx/applications/cpp/phase_estimation.cpp

@@ -3,12 +3,12 @@
 // nvq++ phase_estimation.cpp -o qpe.x && ./qpe.x
 // ```
 
+#include <cmath>
 #include <cudaq.h>
 #include <cudaq/algorithm.h>
+#include <numeric>
 #include <stdio.h>
 
-#include <cmath>
-
 // A pure device quantum kernel defined as a free function
 // (cannot be called from host code).
 __qpu__ void iqft(cudaq::qview<> q) {
@@ -67,9 +67,6 @@ struct qpe {
     // Apply inverse quantum Fourier transform
     iqft(counting_qubits);
 
-    // Measure to gather sampling statistics
-    mz(counting_qubits);
-
     return;
   }
 };
@@ -83,7 +80,10 @@ int main() {
   for (auto nQubits : std::vector<int>{2, 4, 6, 8}) {
     auto counts = cudaq::sample(
         qpe{}, nQubits, [](cudaq::qubit &q) __qpu__ { x(q); }, r1PiGate{});
-    auto mostProbable = counts.most_probable();
+    std::vector<std::size_t> indices(nQubits);
+    std::iota(indices.begin(), indices.end(), 0);
+    auto counting_qubits = counts.get_marginal(indices);
+    auto mostProbable = counting_qubits.most_probable();
     double theta = cudaq::to_integer(mostProbable) / (double)(1UL << nQubits);
     auto piEstimate = 1. / (2 * theta);
     printf("Pi Estimate(nQubits == %d) = %lf \n", nQubits, piEstimate);

docs/sphinx/examples/cpp/basics/bernstein_vazirani.cpp

@@ -1,5 +1,6 @@
 #include <cudaq.h>
 #include <iostream>
+#include <numeric>
 
 // If you have a NVIDIA GPU you can use this example to see
 // that the GPU-accelerated backends can easily handle a
@@ -53,10 +54,6 @@ __qpu__ void bernstein_vazirani(std::vector<int> &hidden_bitstring) {
 
   // Apply another set of Hadamards to the qubits.
   h(qvector);
-
-  // Apply measurement gates to just the `qubits`
-  // (excludes the auxillary qubit).
-  mz(qvector);
 }
 
 int main() {
@@ -69,11 +66,17 @@ int main() {
   auto result = cudaq::sample(bernstein_vazirani, hidden_bitstring);
   result.dump();
 
+  // Extract the sample counts for the qubits, excluding the auxiliary qubit.
+  std::vector<std::size_t> indices(qubit_count);
+  std::iota(indices.begin(), indices.end(), 0);
+  auto counts = result.get_marginal(indices);
+  counts.dump();
+
   std::cout << "Encoded bitstring = ";
   for (auto bit : hidden_bitstring)
     std::cout << bit;
   std::cout << "\n";
-  std::cout << "Measured bitstring = " << result.most_probable() << "\n";
+  std::cout << "Measured bitstring = " << counts.most_probable() << "\n";
 
   return 0;
 }

docs/sphinx/examples/cpp/basics/cuquantum_backends.cpp

@@ -22,7 +22,6 @@ struct ghz {
     for (int i = 0; i < N - 1; i++) {
       x<cudaq::ctrl>(q[i], q[i + 1]);
     }
-    mz(q);
   }
 };
 

docs/sphinx/examples/cpp/basics/cutensornet_backends.cpp

@@ -24,7 +24,6 @@ struct ghz {
     for (int i = 0; i < N - 1; i++) {
       x<cudaq::ctrl>(q[i], q[i + 1]);
     }
-    mz(q);
   }
 };
 

docs/sphinx/examples/cpp/basics/mid_circuit_measurement.cpp

@@ -6,7 +6,7 @@
 #include <cudaq.h>
 
 struct kernel {
-  void operator()() __qpu__ {
+  auto operator()() __qpu__ {
     cudaq::qarray<3> q;
     // Initial state preparation
     x(q[0]);
@@ -26,24 +26,23 @@ struct kernel {
     if (b0)
       z(q[2]);
 
-    mz(q[2]);
+    return mz(q[2]);
   }
 };
 
 int main() {
 
   int nShots = 100;
-  // Sample
-  auto counts = cudaq::sample(/*shots*/ nShots, kernel{});
-  counts.dump();
-
-  // Get the marginal counts on the second qubit
-  auto resultsOnZero = counts.get_marginal({0});
-  resultsOnZero.dump();
-
-  // Count the "1"
-  auto nOnes = resultsOnZero.count("1");
-
+  auto results = cudaq::run(/*shots*/ nShots, kernel{});
+  std::size_t nOnes = 0;
+  // Count the number of times we measured "1"
+  for (auto r : results) {
+    if (r)
+      nOnes++;
+  }
+  // Print out the results
+  printf("Measured '1' on target qubit %zu times out of %d shots.\n", nOnes,
+         nShots);
   // Will fail if not equal to number of shots
   assert(nOnes == nShots && "Failure to teleport qubit in |1> state.");
 }

docs/sphinx/examples/cpp/basics/noise_amplitude_damping.cpp

@@ -34,7 +34,6 @@ int main() {
   auto kernel = []() __qpu__ {
     cudaq::qubit q;
     h(q);
-    mz(q);
   };
 
   // Now let's set the noise and we're ready to run the simulation!

docs/sphinx/examples/cpp/basics/noise_bit_flip.cpp

@@ -32,7 +32,6 @@ int main() {
   auto kernel = []() __qpu__ {
     cudaq::qubit q;
     x(q);
-    mz(q);
   };
 
   // Now let's set the noise and we're ready to run the simulation!

docs/sphinx/examples/cpp/basics/noise_callback.cpp

@@ -44,7 +44,6 @@ int main() {
   auto kernel = [](double angle) __qpu__ {
     cudaq::qubit q;
     rx(angle, q);
-    mz(q);
   };
 
   // Now let's set the noise and we're ready to run the simulation!

docs/sphinx/examples/cpp/basics/noise_depolarization.cpp

@@ -38,7 +38,6 @@ int main() {
   auto kernel = []() __qpu__ {
     cudaq::qubit q;
     y(q);
-    mz(q);
   };
 
   // Now let's set the noise and we're ready to run the simulation!

docs/sphinx/examples/cpp/basics/noise_modeling.cpp

@@ -10,7 +10,6 @@ int main() {
   auto xgate = []() __qpu__ {
     cudaq::qubit q;
     x(q);
-    mz(q);
   };
 
   // Run noise-less simulation

docs/sphinx/examples/cpp/basics/noise_phase_flip.cpp

@@ -37,7 +37,6 @@ int main() {
     z(q);
     // Apply another Hadamard.
     h(q);
-    mz(q);
   };
 
   // Now let's set the noise and we're ready to run the simulation!

docs/sphinx/examples/cpp/basics/static_kernel.cpp

@@ -18,7 +18,6 @@ struct ghz {
     for (int i = 0; i < N - 1; i++) {
       x<cudaq::ctrl>(q[i], q[i + 1]);
     }
-    mz(q);
   }
 };
 

docs/sphinx/examples/cpp/building_kernels.cpp

@@ -105,7 +105,6 @@ __qpu__ void kernel7() {
   x(qvector);
   x(qvector[1]);
   x<cudaq::ctrl>(qvector[0], qvector[1]);
-  mz(qvector);
 }
 
 int main() {

docs/sphinx/examples/cpp/executing_kernels_sample.cpp

@@ -19,9 +19,7 @@ __qpu__ void kernel(int qubit_count) {
   for (auto qubit : cudaq::range(qubit_count - 1)) {
     x<cudaq::ctrl>(qvector[qubit], qvector[qubit + 1]);
   }
-  // If we do not specify measurements, all qubits are measured in
-  // the Z-basis by default or we can manually specify it also
-  mz(qvector);
+  // All qubits are measured in the Z-basis by default
 }
 
 int main() {

docs/sphinx/examples/cpp/executing_kernels_sample_async.cpp

@@ -19,9 +19,7 @@ __qpu__ void kernel(int qubit_count) {
   for (auto qubit : cudaq::range(qubit_count - 1)) {
     x<cudaq::ctrl>(qvector[qubit], qvector[qubit + 1]);
   }
-  // If we do not specify measurements, all qubits are measured in
-  // the Z-basis by default or we can manually specify it also
-  mz(qvector);
+  // All qubits are measured in the Z-basis by default
 }
 
 int main() {

docs/sphinx/examples/cpp/measuring_kernels.cpp

@@ -8,26 +8,26 @@
 
 // [Begin Docs]
 #include <cudaq.h>
+#include <map>
 // [End Docs]
 
 // [Begin Sample1]
-__qpu__ void kernel0() {
+__qpu__ auto kernel0() {
   cudaq::qvector qubits(2);
-  mz(qubits[0]);
+  return mz(qubits[0]);
 }
 // [End Sample1]
 
 // [Begin Sample2]
-__qpu__ void kernel1() {
+__qpu__ auto kernel1() {
   cudaq::qvector qubits_a(2);
   cudaq::qubit qubits_b;
-  mz(qubits_a);
-  mx(qubits_b);
+  return mx(qubits_b);
 }
 // [End Sample2]
 
 // [Begin Sample3]
-__qpu__ void kernel2() {
+__qpu__ auto kernel2() {
   cudaq::qvector q(2);
   h(q[0]);
   auto b0 = mz(q[0]);
@@ -37,18 +37,32 @@ __qpu__ void kernel2() {
   if (b0) {
     h(q[1]);
   }
+  return mz(q);
 }
 
 int main() {
-  auto result = cudaq::sample(kernel2);
-  result.dump();
+  auto results = cudaq::run(1000, kernel2);
+  // Count occurrences of each bitstring
+  std::map<std::string, std::size_t> bitstring_counts;
+  for (const auto &result : results) {
+    std::string bits = std::to_string(result[0]) + std::to_string(result[1]);
+    bitstring_counts[bits]++;
+  }
+
+  printf("Bitstring counts:\n{\n");
+  for (const auto &[bits, count] : bitstring_counts) {
+    printf("  %s: %zu\n", bits.c_str(), count);
+  }
+  printf("}\n");
+
   return 0;
 }
 // [End Sample3]
 
 /* [Begin Sample4]
+Bitstring counts:
 {
-  __global__ : { 10:728 11:272 }
-   b0 : { 0:505 1:495 }
+  10: 771
+  11: 229
 }
  [End Sample4] */