First shot memory (#61)

Change the API of `qec::sample_memory_circuit` to return nRounds of data
as opposed to nRounds - 1 of data. Along with returning nRounds - 1 of
xor'd syndrome data, the initial unxor'd shot is also returned first.

docs/sphinx/components/qec/introduction.rst

@@ -486,14 +486,14 @@ The decoder base class defines the core interface for syndrome decoding:
 
     class decoder {
     protected:
-        std::size_t block_size;     // For [n,k] code, this is n
-        std::size_t syndrome_size;   // For [n,k] code, this is n-k
-        tensor<uint8_t> H;          // Parity check matrix
+        std::size_t block_size;       // For [n,k] code, this is n
+        std::size_t syndrome_size;    // For [n,k] code, this is n-k
+        tensor<uint8_t> H;            // Parity check matrix
 
     public:
         struct decoder_result {
-            bool converged;                  // Decoder convergence status
-            std::vector<float_t> result;     // Soft error probabilities
+            bool converged;                 // Decoder convergence status
+            std::vector<float_t> result;    // Soft error probabilities
         };
 
         virtual decoder_result decode(
@@ -658,20 +658,22 @@ Usage Example
         import cudaq_qec as qec
 
         # Get a code instance
-        code = qec.get_code('steane')
+        steane = qec.get_code("steane")
 
         # Create decoder with code's parity matrix
-        decoder = qec.get_decoder('single_error_lut',
-                                H=code.get_parity())
+        decoder = qec.get_decoder('single_error_lut', steane.get_parity())
 
         # Run stabilizer measurements
-        syndromes, dataQubitResults = qec.sample_memory_circuit(steane, numShots, numRounds)
+        syndromes, dataQubitResults = qec.sample_memory_circuit(steane, numShots=1, numRounds=1)
 
-        # Decode syndrome
+        # Decode a syndrome
         result = decoder.decode(syndromes[0])
         if result.converged:
             print("Error locations:",
                 [i for i,p in enumerate(result.result) if p > 0.5])
+            # No errors as we did not include a noise model and
+            # thus prints:
+            # Error locations: []
 
 .. tab:: C++
 
@@ -783,31 +785,39 @@ Function Variants
 
     .. code-block:: python
 
+        import cudaq
         import cudaq_qec as qec
 
+        # Use the stim backend for performance in QEC settings
+        cudaq.set_target("stim")
+
+        # Get a code instance
+        code = qec.get_code("steane")
+
         # Basic memory circuit with |0⟩ state
         syndromes, measurements = qec.sample_memory_circuit(
             code,           # QEC code instance
-            numShots=1000, # Number of circuit executions
-            numRounds=1    # Number of stabilizer rounds
+            numShots=1000,  # Number of circuit executions
+            numRounds=1     # Number of stabilizer rounds
         )
 
         # Memory circuit with custom initial state
         syndromes, measurements = qec.sample_memory_circuit(
-            code,                       # QEC code instance
-            state_prep=qec.operation.prep1,  # Initial state
+            code,                     # QEC code instance
+            op=qec.operation.prep1,   # Initial state
             numShots=1000,            # Number of shots
             numRounds=1               # Number of rounds
         )
 
         # Memory circuit with noise model
-        noise = cudaq.noise_model()
-        noise.add_channel(...)  # Configure noise
+        noise = cudaq.NoiseModel()
+        # Configure noise
+        noise.add_all_qubit_channel("x", qec.TwoQubitDepolarization(0.01), 1)
         syndromes, measurements = qec.sample_memory_circuit(
-            code,           # QEC code instance
-            numShots=1000, # Number of shots
-            numRounds=1,   # Number of rounds
-            noise=noise     # Noise model
+            code,             # QEC code instance
+            numShots=1000,    # Number of shots
+            numRounds=1,      # Number of rounds
+            noise=noise       # Noise model
         )
 
 .. tab:: C++
@@ -816,27 +826,27 @@ Function Variants
 
         // Basic memory circuit with |0⟩ state
         auto [syndromes, measurements] = qec::sample_memory_circuit(
-            code,           // QEC code instance
-            numShots,       // Number of circuit executions
-            numRounds       // Number of stabilizer rounds
+            code,       // QEC code instance
+            numShots,   // Number of circuit executions
+            numRounds   // Number of stabilizer rounds
         );
 
         // Memory circuit with custom initial state
         auto [syndromes, measurements] = qec::sample_memory_circuit(
-            code,                    // QEC code instance
-            operation::prep1,  // Initial state preparation
-            numShots,                // Number of circuit executions
-            numRounds               // Number of stabilizer rounds
+            code,               // QEC code instance
+            operation::prep1,   // Initial state preparation
+            numShots,           // Number of circuit executions
+            numRounds           // Number of stabilizer rounds
         );
 
         // Memory circuit with noise model
         auto noise_model = cudaq::noise_model();
         noise_model.add_channel(...);  // Configure noise
         auto [syndromes, measurements] = qec::sample_memory_circuit(
-            code,           // QEC code instance
-            numShots,       // Number of circuit executions
-            numRounds,      // Number of stabilizer rounds
-            noise_model     // Noise model to apply
+            code,         // QEC code instance
+            numShots,     // Number of circuit executions
+            numRounds,    // Number of stabilizer rounds
+            noise_model   // Noise model to apply
         );
 
 Return Values
@@ -846,7 +856,7 @@ The functions return a tuple containing:
 
 1. **Syndrome Measurements** (:code:`tensor<uint8_t>`):
 
-   * Shape: :code:`(num_shots, (num_rounds-1) * syndrome_size)`
+   * Shape: :code:`(num_shots, num_rounds * syndrome_size)`
    * Contains stabilizer measurement results
    * Values are 0 or 1 representing measurement outcomes
 
@@ -859,7 +869,7 @@ The functions return a tuple containing:
 Example Usage
 ~~~~~~~~~~~~~
 
-Here's a complete example of running a memory experiment:
+Example of running a memory experiment:
 
 .. tab:: Python
 
@@ -868,21 +878,24 @@ Here's a complete example of running a memory experiment:
         import cudaq
         import cudaq_qec as qec
 
+        # Use the stim backend for performance in QEC settings
+        cudaq.set_target("stim")
+
         # Create code and decoder
         code = qec.get_code('steane')
         decoder = qec.get_decoder('single_error_lut',
                                   code.get_parity())
 
         # Configure noise
-        noise = cudaq.noise_model()
-        noise.add_channel('x', depolarizing=0.001)
+        noise = cudaq.NoiseModel()
+        noise.add_all_qubit_channel("x", qec.TwoQubitDepolarization(0.01), 1)
 
         # Run memory experiment
         syndromes, measurements = qec.sample_memory_circuit(
             code,
-            state_prep=qec.operation.prep0,
-            num_shots=1000,
-            num_rounds=10,
+            op=qec.operation.prep0,
+            numShots=1000,
+            numRounds=10,
             noise=noise
         )
 
@@ -901,38 +914,46 @@ Here's a complete example of running a memory experiment:
 
     .. code-block:: cpp
 
-        // Top of file
-        #include "cudaq/qec/experiments.h"
-
-        // Create a Steane code instance
-        auto code = cudaq::qec::get_code("steane");
-
-        // Configure noise model
-        auto noise = cudaq::noise_model();
-        noise.add_all_qubit_channel("x", cudaq::qec::two_qubit_depolarization(0.1),
-                            /*num_controls=*/1);
+        // Compile and run with:
+        // nvq++ --enable-mlir --target=stim -lcudaq-qec example.cpp
+        // ./a.out
 
-        // Run memory experiment
-        auto [syndromes, measurements] = qec::sample_memory_circuit(
-            code,                    // Code instance
-            operation::prep0,  // Prepare |0⟩ state
-            1000,                    // 1000 shots
-            10,                      // 10 rounds
-            noise                    // Apply noise
-        );
-
-        // Analyze results
-        auto decoder = qec::get_decoder("single_error_lut", code->get_parity());
-        for (std::size_t shot = 0; shot < 1000; shot++) {
+        #include "cudaq.h"
+        #include "cudaq/qec/decoder.h"
+        #include "cudaq/qec/experiments.h"
+        #include "cudaq/qec/noise_model.h"
+
+        int main(){
+          // Create a Steane code instance
+          auto code = cudaq::qec::get_code("steane");
+
+          // Configure noise model
+          cudaq::noise_model noise;
+          noise.add_all_qubit_channel("x", cudaq::qec::two_qubit_depolarization(0.1),
+                              /*num_controls=*/1);
+
+          // Run memory experiment
+          auto [syndromes, data] = cudaq::qec::sample_memory_circuit(
+              *code,                          // Code instance
+              cudaq::qec::operation::prep0,   // Prepare |0⟩ state
+              1000,                           // 1000 shots
+              1,                              // 1 rounds
+              noise                           // Apply noise
+          );
+
+          // Analyze results
+          auto decoder = cudaq::qec::get_decoder("single_error_lut", code->get_parity());
+          for (std::size_t shot = 0; shot < 1000; shot++) {
             // Get syndrome for this shot
-            std::vector<float> syndrome(code->get_syndrome_size());
+            std::vector<cudaq::qec::float_t> syndrome(syndromes.shape()[1]);
             for (std::size_t i = 0; i < syndrome.size(); i++)
-                syndrome[i] = syndromes.at({shot, i});
+              syndrome[i] = syndromes.at({shot, i});
 
             // Decode syndrome
-            auto result = decoder->decode(syndrome);
+            auto [converged, v_result] = decoder->decode(syndrome);
             // Process correction
             // ...
+          }
         }
 
 Additional Noise Models:

docs/sphinx/examples/qec/cpp/circuit_level_noise.cpp

@@ -6,9 +6,9 @@
  * the terms of the Apache License 2.0 which accompanies this distribution.    *
  ******************************************************************************/
 
-// Compile and run with
-// nvq++ --enable-mlir -lcudaq-qec circuit_level_noise.cpp -o circuit_level
-// ./circuit_level
+// Compile and run with:
+// nvq++ --enable-mlir --target=stim -lcudaq-qec circuit_level_noise.cpp
+// ./a.out
 
 #include "cudaq.h"
 #include "cudaq/qec/decoder.h"
@@ -80,11 +80,11 @@ int main() {
   for (size_t shot = 0; shot < nShots; ++shot) {
     std::cout << "shot: " << shot << "\n";
 
-    for (size_t round = 0; round < nRounds - 1; ++round) {
+    for (size_t round = 0; round < nRounds; ++round) {
       std::cout << "round: " << round << "\n";
 
       // Access one row of the syndrome tensor
-      size_t count = shot * (nRounds - 1) + round;
+      size_t count = shot * nRounds + round;
       size_t stride = syndromes.shape()[1];
       cudaqx::tensor<uint8_t> syndrome({stride});
       syndrome.borrow(syndromes.data() + stride * count);

docs/sphinx/examples/qec/cpp/code_capacity_noise.cpp

@@ -10,9 +10,8 @@
 // decoder, code
 //
 // Compile and run with
-// nvq++ --enable-mlir -lcudaq-qec code_capacity_noise.cpp -o
-// code_capacity_noise
-// ./code_capacity_noise
+// nvq++ --enable-mlir --target=stim -lcudaq-qec code_capacity_noise.cpp
+// ./a.out
 
 #include <algorithm>
 #include <cmath>

docs/sphinx/examples_rst/qec/circuit_level_noise.rst

@@ -36,7 +36,7 @@ Here's how to use CUDA-Q QEC to perform a circuit-level noise model experiment i
 
    .. code-block:: bash
 
-      nvq++ --enable-mlir -lcudaq-qec circuit_level_noise.cpp -o circuit_level_noise
+      nvq++ --enable-mlir --target=stim -lcudaq-qec circuit_level_noise.cpp -o circuit_level_noise
       ./circuit_level_noise
 
 
@@ -58,8 +58,8 @@ Here's how to use CUDA-Q QEC to perform a circuit-level noise model experiment i
     - Each memory circuit runs for an input number of `nRounds`, which specifies how many `stabilizer_round` kernels are ran.
     - After `nRounds` the data qubits are measured and the run is over.
     - This is performed `nShots` number of times.
-    - During a shot, each syndrome is `xor`'d against the preceding syndrome, so that we can track a sparser flow of data showing which round each parity check was violated.
-    - This means we will get a total of `nShots * (nRounds - 1)` syndromes to decode and analyze.
+    - During a shot, each stabilizer round's syndrome is `xor`'d against the preceding syndrome, so that we can track a sparser flow of data showing which round each parity check was violated.
+    - The first round returns the syndrome as is, as there is nothing preceding to `xor` against.
 
 5. Data qubit measurements:
     - The data qubits are only read out after the end of each shot, so there are `nShots` worth of data readouts.

docs/sphinx/examples_rst/qec/code_capacity_noise.rst

@@ -46,7 +46,7 @@ Here's how to use CUDA-Q QEC to perform a code capacity noise model experiment i
 
    .. code-block:: bash
 
-      nvq++ --enable-mlir -lcudaq-qec code_capacity_noise.cpp -o code_capacity_noise
+      nvq++ --enable-mlir --target=stim -lcudaq-qec code_capacity_noise.cpp -o code_capacity_noise
       ./code_capacity_noise
 
 

libs/qec/lib/experiments.cpp

@@ -119,7 +119,7 @@ sample_memory_circuit(const code &code, operation statePrep,
   std::size_t numCols = numAncx + numAncz;
 
   // Allocate the tensor data for the syndromes and data.
-  cudaqx::tensor<uint8_t> syndromeTensor({numShots * (numRounds - 1), numCols});
+  cudaqx::tensor<uint8_t> syndromeTensor({numShots * numRounds, numCols});
   cudaqx::tensor<uint8_t> dataResults({numShots, numData});
 
   // Run the memory circuit experiment
@@ -151,20 +151,25 @@ sample_memory_circuit(const code &code, operation statePrep,
   auto &measurements = getMemoryCircuitAncillaMeasurements();
 
   std::size_t numMeasRows = numShots * numRounds;
-  std::size_t numSyndRows = numShots * (numRounds - 1);
   cudaqx::tensor<uint8_t> measuresTensor({numMeasRows, numCols});
   measuresTensor.borrow(measurements.data());
 
-  // Convert to Syndromes
+  // First round, store bare syndrome measurement
+  for (std::size_t col = 0; col < numCols; ++col) {
+    std::size_t shot = 0;
+    std::size_t round = 0;
+    std::size_t measIdx = shot * numRounds + round;
+    syndromeTensor.at({measIdx, col}) = measuresTensor.at({measIdx, col});
+  }
 
+  // After first round, store syndrome flips
   // #pragma omp parallel for collapse(2)
   for (std::size_t shot = 0; shot < numShots; ++shot)
     for (std::size_t round = 1; round < numRounds; ++round)
       for (std::size_t col = 0; col < numCols; ++col) {
         std::size_t measIdx = shot * numRounds + round;
         std::size_t prevMeasIdx = shot * numRounds + (round - 1);
-        std::size_t syndIdx = shot * (numRounds - 1) + (round - 1);
-        syndromeTensor.at({syndIdx, col}) =
+        syndromeTensor.at({measIdx, col}) =
             measuresTensor.at({measIdx, col}) ^
             measuresTensor.at({prevMeasIdx, col});
       }

libs/qec/python/tests/test_code.py

@@ -58,14 +58,14 @@ def test_sample_memory_circuit():
                                                        numShots=10,
                                                        numRounds=4)
     assert isinstance(syndromes, np.ndarray)
-    assert syndromes.shape == (30, 6)
+    assert syndromes.shape == (40, 6)
     print(syndromes)
 
     syndromes_with_op, dataResults = qec.sample_memory_circuit(
         steane, qec.operation.prep1, 10, 4)
     assert isinstance(syndromes_with_op, np.ndarray)
     print(syndromes_with_op)
-    assert syndromes_with_op.shape == (30, 6)
+    assert syndromes_with_op.shape == (40, 6)
 
 
 def test_custom_steane_code():
@@ -99,7 +99,7 @@ def test_noisy_simulation():
                                                        numRounds=4,
                                                        noise=noise)
     assert isinstance(syndromes, np.ndarray)
-    assert syndromes.shape == (30, 6)
+    assert syndromes.shape == (40, 6)
     print(syndromes)
     assert np.any(syndromes)
     cudaq.reset_target()
@@ -111,7 +111,7 @@ def test_python_code():
                                                        numShots=10,
                                                        numRounds=4)
     assert isinstance(syndromes, np.ndarray)
-    assert syndromes.shape == (30, 6)
+    assert syndromes.shape == (40, 6)
     print(syndromes)
     assert not np.any(syndromes)