Fix/2q rb

CHANGELOG.md

@@ -4,6 +4,11 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/)
 
 ## [Unreleased]
+### Added
+- two-qubit rb - Added feature to plot the two qubit state distribution.
+
+### Fixed
+- two-qubit rb - Swapped the order of the circuit_depth and repeat axis for better performance.
 
 ## [0.19.1] - 2025-02-17
 ### Changed

qualang_tools/characterization/two_qubit_rb/two_qubit_rb_example.py → examples/two_qubit_rb.py

@@ -138,6 +138,9 @@ def meas():
 res.plot_fidelity()
 plt.show()
 
+res.plot_two_qubit_state_distribution()
+plt.show()
+
 # verify/save the random sequences created during the experiment
 rb.save_sequences_to_file("sequences.txt")  # saves the gates used in each random sequence
 rb.save_command_mapping_to_file("commands.txt")  # saves mapping from "command id" to sequence

qualang_tools/characterization/two_qubit_rb/README.md

@@ -37,7 +37,7 @@ res = rb.run(
 )
 ```
 
-[Here](two_qubit_rb_example.py) is an example for flux-tunable transmon qubits.
+[Here](/examples/two_qubit_rb.py) is an example for flux-tunable transmon qubits.
 
 ## Introduction
 Two-Qubit Randomized Benchmarking (RB) has become a popular protocol that allows to experimentally quantify the performance of a quantum processor by applying sequences of randomly sampled Clifford gates and measuring the average error rate. Due to its universality it has been implemented in various qubit platforms such as trapped-ions [^1], NMR [^2], spin [^3] and superconducting qubits [^4]. Two-Qubit RB can be challenging to implement with state-of-the-art control electronics because of the necessity to sample from a large Clifford gate set. The Clifford group consists of 11520 operations [^4] and contains the single qubit Clifford operations (576), the CNOT-like class (5184), the iSWAP-like class (5184) and the SWAP-like class (576). In the provided example we introduce an implementation on the OPX+ using the current version of the generic `TwoQubitRb` class. The implementation exploits the [baking](https://github.com/qua-platform/py-qua-tools/blob/main/qualang_tools/bakery/README.md) tool to generate the individual Clifford operations. The class then uses the [Input Stream](https://docs.quantum-machines.co/latest/qm-qua-sdk/docs/Guides/features/?h=declare_input_stream#input-streams) feature to send a string of Clifford indices to the OPX that represent the executed gate sequence which is terminated with the inverse operation. The execution is based on the [Switch Case](https://docs.quantum-machines.co/latest/qm-qua-sdk/docs/Guides/features/?h=switch#switch-case) flow control of QUA, which sets the current minimal gate duration limit to 40 ns.  <!--The inverse is calculated in Python using Clifford tableaus.

qualang_tools/characterization/two_qubit_rb/__init__.py

@@ -1,3 +1,4 @@
 from .two_qubit_rb import *
+from .two_qubit_rb.RBResult import RBResult
 
-__all__ = ["TwoQubitRb", "TwoQubitRbDebugger"]
+__all__ = ["TwoQubitRb", "TwoQubitRbDebugger", "RBResult"]

qualang_tools/characterization/two_qubit_rb/two_qubit_rb/RBResult.py

@@ -27,10 +27,10 @@ def __post_init__(self):
         Initializes the xarray Dataset to store the RB experiment data.
         """
         self.data = xr.Dataset(
-            data_vars={"state": (["circuit_depth", "repeat", "average"], self.state)},
+            data_vars={"state": (["repeat", "circuit_depth", "average"], self.state)},
             coords={
-                "circuit_depth": self.circuit_depths,
                 "repeat": range(self.num_repeats),
+                "circuit_depth": self.circuit_depths,
                 "average": range(self.num_averages),
             },
         )
@@ -64,23 +64,93 @@ def plot(self):
     def plot_with_fidelity(self):
         """
         Plots the RB fidelity as a function of circuit depth, including a fit to an exponential decay model.
-        The fitted curve is overlaid with the raw data points.
+        The fitted curve is overlaid with the raw data points, and error bars are included.
         """
         A, alpha, B = self.fit_exponential()
         fidelity = self.get_fidelity(alpha)
 
+        # Compute error bars
+        error_bars = (self.data == 0).mean(dim="average").std(dim="repeat").state.data
+
         plt.figure()
-        plt.plot(self.circuit_depths, self.get_decay_curve(), "o", label="Data")
+        plt.errorbar(
+            self.circuit_depths,
+            self.get_decay_curve(),
+            yerr=error_bars,
+            fmt=".",
+            capsize=2,
+            elinewidth=0.5,
+            color="blue",
+            label="Experimental Data",
+        )
+
+        circuit_depths_smooth_axis = np.linspace(self.circuit_depths[0], self.circuit_depths[-1], 100)
+        plt.plot(
+            circuit_depths_smooth_axis,
+            rb_decay_curve(np.array(circuit_depths_smooth_axis), A, alpha, B),
+            color="red",
+            linestyle="--",
+            label="Exponential Fit",
+        )
+
+        plt.text(
+            0.5,
+            0.95,
+            f"2Q Clifford Fidelity = {fidelity * 100:.2f}%",
+            horizontalalignment="center",
+            verticalalignment="top",
+            fontdict={"fontsize": "large", "fontweight": "bold"},
+            transform=plt.gca().transAxes,
+        )
+
+        plt.xlabel("Circuit Depth")
+        plt.ylabel(r"Probability to recover to $|00\rangle$")
+        plt.title("2Q Randomized Benchmarking")
+        plt.legend(framealpha=0)
+        plt.show()
+
+    def plot_two_qubit_state_distribution(self):
+        """
+        Plot how the two-qubit state is distributed as a function of circuit-depth on average.
+        """
         plt.plot(
             self.circuit_depths,
-            rb_decay_curve(np.array(self.circuit_depths), A, alpha, B),
-            "-",
-            label=f"Fidelity={fidelity * 100:.3f}%\nalpha={alpha:.4f}",
+            (self.data.state == 0).mean(dim="average").mean(dim="repeat").data,
+            label=r"$|00\rangle$",
+            marker=".",
+            color="c",
+            linewidth=3,
         )
+        plt.plot(
+            self.circuit_depths,
+            (self.data.state == 1).mean(dim="average").mean(dim="repeat").data,
+            label=r"$|01\rangle$",
+            marker=".",
+            color="b",
+            linewidth=1,
+        )
+        plt.plot(
+            self.circuit_depths,
+            (self.data.state == 2).mean(dim="average").mean(dim="repeat").data,
+            label=r"$|10\rangle$",
+            marker=".",
+            color="y",
+            linewidth=1,
+        )
+        plt.plot(
+            self.circuit_depths,
+            (self.data.state == 3).mean(dim="average").mean(dim="repeat").data,
+            label=r"$|11\rangle$",
+            marker=".",
+            color="r",
+            linewidth=1,
+        )
+        plt.axhline(0.25, color="grey", linestyle="--", linewidth=2, label="2Q mixed-state")
+
         plt.xlabel("Circuit Depth")
-        plt.ylabel("Fidelity")
-        plt.title("2Q Randomized Benchmarking Fidelity")
-        plt.legend()
+        plt.ylabel(r"Probability to recover to a given state")
+        plt.title("2Q State Distribution vs. Circuit Depth")
+        plt.legend(framealpha=0, title=r"2Q State $\mathbf{|q_cq_t\rangle}$", title_fontproperties={"weight": "bold"})
         plt.show()
 
     def fit_exponential(self):
@@ -95,7 +165,7 @@ def fit_exponential(self):
         """
         decay_curve = self.get_decay_curve()
 
-        popt, _ = curve_fit(rb_decay_curve, self.circuit_depths, decay_curve, p0=[0.75, -0.1, 0.25], maxfev=10000)
+        popt, _ = curve_fit(rb_decay_curve, self.circuit_depths, decay_curve, p0=[0.75, 0.9, 0.25], maxfev=10000)
         A, alpha, B = popt
 
         return A, alpha, B

qualang_tools/characterization/two_qubit_rb/two_qubit_rb/TwoQubitRB.py

@@ -171,8 +171,8 @@ def _gen_qua_program(self, sequence_depths: list[int], num_repeats: int, num_ave
             }
 
             assign(progress, 0)
-            with for_each_(sequence_depth, sequence_depths):
-                with for_(repeat, 0, repeat < num_repeats, repeat + 1):
+            with for_(repeat, 0, repeat < num_repeats, repeat + 1):
+                with for_each_(sequence_depth, sequence_depths):
                     assign(progress, progress + 1)
                     save(progress, progress_os)
                     advance_input_stream(gates_len_is)
@@ -187,8 +187,9 @@ def _gen_qua_program(self, sequence_depths: list[int], num_repeats: int, num_ave
                         save(state, state_os)
 
             with stream_processing():
-                state_os.buffer(len(sequence_depths), num_repeats, num_averages).save("state")
+                state_os.buffer(num_repeats, len(sequence_depths), num_averages).save("state")
                 progress_os.save("progress")
+
         return prog
 
     def _input_stream_name(self, element: str):
@@ -208,14 +209,16 @@ def _insert_all_input_stream(
         num_repeats: int,
         callback: Optional[Callable[[List[int]], None]] = None,
     ):
-        for sequence_depth in sequence_depths:
-            for repeat in range(num_repeats):
+        for repeat in range(num_repeats):
+            for sequence_depth in sequence_depths:
                 sequence = self._gen_rb_sequence(sequence_depth)
                 if self._sequence_tracker is not None:
                     self._sequence_tracker.make_sequence(sequence)
                 job.insert_input_stream("__gates_len_is__", len(sequence))
                 for qe in self._rb_baker.all_elements:
-                    job.insert_input_stream(f"{qe}_is", self._decode_sequence_for_element(qe, sequence))
+                    job.insert_input_stream(
+                        f"{self._input_stream_name(qe)}_is", self._decode_sequence_for_element(qe, sequence)
+                    )
 
                 if callback is not None:
                     callback(sequence)