Skip to content

Randomized Benchmarking

Welcome, reader! This chapter will shine a light on:

  1. QUA embedded in Python
    • QUA Macros
  2. Real-Time Processing
    • Random Number Generation
  3. Real-Time Control Flow:
    • Switch/Case

Shall we commence?

The QUA program

QUA macros improve code readability and promote reuse. In our real-time RB implementation, macros condense numerous lines related to random circuit generation, play, and qubit measurement, offering a clear sequence overview.

rb.py
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
# Go to line 35, 72, 92, 187, 195, and 197
import qm.qua as qua
from qm import QuantumMachinesManager
from qm import SimulationConfig, LoopbackInterface
from scipy.optimize import curve_fit
from OPX1000configuration import config, qop_ip, cluster_name, cloud_username, cloud_password # or OPXplusconfiguration, depending on your hardware    import matplotlib.pyplot as plt
import numpy as np
from qualang_tools.results import progress_counter, fetching_tool
from qualang_tools.plot import interrupt_on_close
from qualang_tools.bakery.randomized_benchmark_c1 import c1_table
from macros import readout_macro
import warnings

warnings.filterwarnings("ignore")

##############################
# Program-specific variables #
##############################

num_of_sequences = 3  # Number of random sequences
n_avg = 1  # Number of averaging loops for each random sequence
max_circuit_depth = 28  # Maximum circuit depth
delta_clifford = 2  #  Play each sequence with a depth step equals to 'delta_clifford - Must be > 1
assert (max_circuit_depth / delta_clifford).is_integer(), "max_circuit_depth / delta_clifford must be an integer."
seed = 345324  # Pseudo-random number generator seed
# Flag to enable state discrimination if the readout has been calibrated (rotated blobs and threshold)
state_discrimination = False
# List of recovery gates from the lookup table
inv_gates = [int(np.where(c1_table[i, :] == 0)[0][0]) for i in range(24)]


###################################
# Helper functions and QUA macros #
###################################
def readout_macro(threshold=None, state=None, I=None, Q=None): # (1)!
    """
    A macro for performing the readout, with the ability to perform state discrimination.
    If `threshold` is given, the information in the `I` quadrature will be compared against the threshold and `state`
    would be `True` if `I > threshold`.
    Note that it is assumed that the results are rotated such that all the information is in the `I` quadrature.

    :param threshold: Optional. The threshold to compare `I` against.
    :param state: A QUA variable for the state information, only used when a threshold is given.
        Should be of type `bool`. If not given, a new variable will be created
    :param I: A QUA variable for the information in the `I` quadrature. Should be of type `Fixed`. If not given, a new
        variable will be created
    :param Q: A QUA variable for the information in the `Q` quadrature. Should be of type `Fixed`. If not given, a new
        variable will be created
    :return: Three QUA variables populated with the results of the readout: (`state`, `I`, `Q`)
    """
    if I is None:
        I = qua.declare(qua.fixed)
    if Q is None:
        Q = qua.declare(qua.fixed)
    if threshold is not None and state is None:
        state = qua.declare(bool)
    qua.measure(
        "readout",
        "resonator",
        None,
        qua.dual_demod.full("rotated_cos", "rotated_sin", I),
        qua.dual_demod.full("rotated_minus_sin", "rotated_cos", Q),
    )
    if threshold is not None:
        qua.assign(state, I > threshold)
    return state, I, Q

def power_law(power, a, b, p):
    return a * (p**power) + b


def generate_sequence():  # (1)!
    cayley = qua.declare(int, value=c1_table.flatten().tolist())
    inv_list = qua.declare(int, value=inv_gates)
    current_state = qua.declare(int)
    step = qua.declare(int)
    sequence = qua.declare(int, size=max_circuit_depth + 1)
    inv_gate = qua.declare(int, size=max_circuit_depth + 1)
    i = qua.declare(int)
    rand = qua.Random(seed=seed)

    qua.assign(current_state, 0)
    with qua.for_(i, 0, i < max_circuit_depth, i + 1):
        qua.assign(step, rand.rand_int(24))
        qua.assign(current_state, cayley[current_state * 24 + step])
        qua.assign(sequence[i], step)
        qua.assign(inv_gate[i], inv_list[current_state])

    return sequence, inv_gate


def play_sequence(sequence_list, depth): # (1)!
    i = qua.declare(int)
    with qua.for_(i, 0, i <= depth, i + 1):
        with qua.switch_(sequence_list[i], unsafe=True):
            with qua.case_(0):
                qua.wait(x180_len // 4, "qubit")
            with qua.case_(1):
                qua.play("x180", "qubit")
            with qua.case_(2):
                qua.play("y180", "qubit")
            with qua.case_(3):
                qua.play("y180", "qubit")
                qua.play("x180", "qubit")
            with qua.case_(4):
                qua.play("x90", "qubit")
                qua.play("y90", "qubit")
            with qua.case_(5):
                qua.play("x90", "qubit")
                qua.play("-y90", "qubit")
            with qua.case_(6):
                qua.play("-x90", "qubit")
                qua.play("y90", "qubit")
            with qua.case_(7):
                qua.play("-x90", "qubit")
                qua.play("-y90", "qubit")
            with qua.case_(8):
                qua.play("y90", "qubit")
                qua.play("x90", "qubit")
            with qua.case_(9):
                qua.play("y90", "qubit")
                qua.play("-x90", "qubit")
            with qua.case_(10):
                qua.play("-y90", "qubit")
                qua.play("x90", "qubit")
            with qua.case_(11):
                qua.play("-y90", "qubit")
                qua.play("-x90", "qubit")
            with qua.case_(12):
                qua.play("x90", "qubit")
            with qua.case_(13):
                qua.play("-x90", "qubit")
            with qua.case_(14):
                qua.play("y90", "qubit")
            with qua.case_(15):
                qua.play("-y90", "qubit")
            with qua.case_(16):
                qua.play("-x90", "qubit")
                qua.play("y90", "qubit")
                qua.play("x90", "qubit")
            with qua.case_(17):
                qua.play("-x90", "qubit")
                qua.play("-y90", "qubit")
                qua.play("x90", "qubit")
            with qua.case_(18):
                qua.play("x180", "qubit")
                qua.play("y90", "qubit")
            with qua.case_(19):
                qua.play("x180", "qubit")
                qua.play("-y90", "qubit")
            with qua.case_(20):
                qua.play("y180", "qubit")
                qua.play("x90", "qubit")
            with qua.case_(21):
                qua.play("y180", "qubit")
                qua.play("-x90", "qubit")
            with qua.case_(22):
                qua.play("x90", "qubit")
                qua.play("y90", "qubit")
                qua.play("x90", "qubit")
            with qua.case_(23):
                qua.play("-x90", "qubit")
                qua.play("y90", "qubit")
                qua.play("-x90", "qubit")


###################
# The QUA program #
###################
with qua.program() as rb:
    depth = qua.declare(int)
    depth_target = qua.declare(int) 
    saved_gate = qua.declare(int)
    m = qua.declare(int)
    n = qua.declare(int)
    I = qua.declare(qua.fixed)
    Q = qua.declare(qua.fixed)
    state = qua.declare(bool)
    m_st = qua.declare_stream()
    if state_discrimination:
        state_st = qua.declare_stream()
    else:
        I_st = qua.declare_stream()
        Q_st = qua.declare_stream()

    with qua.for_(m, 0, m < num_of_sequences, m + 1):
        sequence_list, inv_gate_list = generate_sequence()

        qua.assign(depth_target, 0) 
        with qua.for_(depth, 1, depth <= max_circuit_depth, depth + 1):
            qua.assign(saved_gate, sequence_list[depth])
            qua.assign(sequence_list[depth], inv_gate_list[depth - 1])
            with qua.if_((depth == 1) | (depth == depth_target)):
                with qua.for_(n, 0, n < n_avg, n + 1):
                    play_sequence(sequence_list, depth)
                    qua.align("qubit", "resonator")
                    state, I, Q = readout_macro(threshold=ge_threshold, state=state, I=I, Q=Q)
                    if state_discrimination:
                        qua.save(state, state_st)
                    else:
                        qua.save(I, I_st)
                        qua.save(Q, Q_st)
                    qua.align("resonator", "qubit")
                    qua.wait(thermalization_time * u.ns, "resonator")
                qua.assign(depth_target, depth_target + delta_clifford)
            qua.assign(sequence_list[depth], saved_gate)
        qua.save(m, m_st)

    with qua.stream_processing():
        m_st.save("iteration")
        if state_discrimination:
            state_st.boolean_to_int()
                    .buffer(n_avg)
                    .map(FUNCTIONS.average())
                    .buffer(max_circuit_depth / delta_clifford + 1)
                    .buffer(num_of_sequences)
                    .save("state")
            state_st.boolean_to_int()
                    .buffer(n_avg)
                    .map(FUNCTIONS.average())
                    .buffer(max_circuit_depth / delta_clifford + 1)
                    .average()
                    .save("state_avg")
        else:
            I_st.buffer(n_avg)
                .map(FUNCTIONS.average())
                .buffer(max_circuit_depth / delta_clifford + 1)
                .buffer(num_of_sequences)
                .save("I")
            Q_st.buffer(n_avg)
                .map(FUNCTIONS.average())
                .buffer(max_circuit_depth / delta_clifford + 1)
                .buffer(num_of_sequences)
                .save("Q")
            I_st.buffer(n_avg)
                .map(FUNCTIONS.average())
                .buffer(max_circuit_depth / delta_clifford + 1)
                .average()
                .save("I_avg")
            Q_st.buffer(n_avg)
                .map(FUNCTIONS.average())
                .buffer(max_circuit_depth / delta_clifford + 1)
                .average()
                .save("Q_avg")


#####################################
#  Open Communication with the QOP  #
#####################################
simulate = True
simulation_in_cloud = True

if simulate: 

    if simulation_in_cloud:
        client = QmSaas(email=cloud_username, password=cloud_password)
        instance = client.simulator(QoPVersion.v3_2_0)
        instance.spawn()
        qmm = QuantumMachinesManager(host=instance.host,
                                    port=instance.port,
                                    connection_headers=instance.default_connection_headers)
    else:
        qmm = QuantumMachinesManager(host=qop_ip, 
                                    cluster_name=cluster_name)

    simulation_config = SimulationConfig(
            duration=50_000
        )
    job = qmm.simulate(config, rb, simulation_config)
    job.get_simulated_samples().con1.plot()
    plt.figure()
    if state_discrimination:
        results = fetching_tool(job, data_list=["state_avg", "iteration"])
        state_avg, iteration = results.fetch_all()
        value_avg = state_avg
    else:
        results = fetching_tool(job, data_list=["I_avg", "Q_avg", "iteration"])
        I, Q, iteration = results.fetch_all()
        value_avg = I
    x = np.arange(0, max_circuit_depth + 0.1, delta_clifford)
    x[0] = 1 
    plt.plot(x, value_avg, marker=".")
    plt.xlabel("Number of Clifford gates")
    plt.ylabel("Sequence Fidelity")
    plt.title("Single qubit RB")
    plt.tight_layout()
    plt.show()
else:
    qm = qmm.open_qm(config)
    job = qm.execute(rb)
    if state_discrimination:
        results = fetching_tool(job, data_list=["state_avg", "iteration"], mode="live")
    else:
        results = fetching_tool(job, data_list=["I_avg", "Q_avg", "iteration"], mode="live")
    fig = plt.figure()
    interrupt_on_close(fig, job) 
    x = np.arange(0, max_circuit_depth + 0.1, delta_clifford)
    x[0] = 1 
    while results.is_processing():
        if state_discrimination:
            state_avg, iteration = results.fetch_all()
            value_avg = state_avg
        else:
            I, Q, iteration = results.fetch_all()
            value_avg = I

        progress_counter(iteration, num_of_sequences, start_time=results.get_start_time())
        plt.cla()
        plt.plot(x, value_avg, marker=".")
        plt.xlabel("Number of Clifford gates")
        plt.ylabel("Sequence Fidelity")
        plt.title("Single qubit RB")
        plt.pause(0.1)
    plt.show()

    if state_discrimination:
        results = fetching_tool(job, data_list=["state"])
        state = results.fetch_all()[0]
        value_avg = np.mean(state, axis=0)
        error_avg = np.std(state, axis=0)
    else:
        results = fetching_tool(job, data_list=["I", "Q"])
        I, Q = results.fetch_all()
        value_avg = np.mean(I, axis=0)
        error_avg = np.std(I, axis=0)
    pars, cov = curve_fit(
        f=power_law,
        xdata=x,
        ydata=value_avg,
        p0=[0.5, 0.5, 0.9],
        bounds=(-np.inf, np.inf),
        maxfev=2000,
    )
    stdevs = np.sqrt(np.diag(cov))

    print("#########################")
    print("### Fitted Parameters ###")
    print("#########################")
    print(f"A = {pars[0]:.3} ({stdevs[0]:.1}), B = {pars[1]:.3} ({stdevs[1]:.1}), p = {pars[2]:.3} ({stdevs[2]:.1})")
    print("Covariance Matrix")
    print(cov)

    one_minus_p = 1 - pars[2]
    r_c = one_minus_p * (1 - 1 / 2**1)
    r_g = r_c / 1.875  # 1.875 is the average number of gates in clifford operation
    r_c_std = stdevs[2] * (1 - 1 / 2**1)
    r_g_std = r_c_std / 1.875

    print("#########################")
    print("### Useful Parameters ###")
    print("#########################")
    print(
        f"Error rate: 1-p = {np.format_float_scientific(one_minus_p, precision=2)} ({stdevs[2]:.1})\n"
        f"Clifford set infidelity: r_c = {np.format_float_scientific(r_c, precision=2)} ({r_c_std:.1})\n"
        f"Gate infidelity: r_g = {np.format_float_scientific(r_g, precision=2)}  ({r_g_std:.1})"
    )

    # Plots
    plt.figure()
    plt.errorbar(x, value_avg, yerr=error_avg, marker=".")
    plt.plot(x, power_law(x, *pars), linestyle="--", linewidth=2)
    plt.xlabel("Number of Clifford gates")
    plt.ylabel("Sequence Fidelity")
    plt.title("Single qubit RB")
  1. QUA macros have a similarity to Python functions but the indentend lines are indeed QUA code that is put into place in the process of translation and sending from your computer to the compiler. There are good practices to follow when working with QUA macros.

The PPU and QUA gives you the capability to generate random numbers and store them in QUA variables. This feature ensures that QUA is mathematically comprehensive, advanced, and practical. In RB it is used to generate the random Clifford-circuit in real-time.

rb.py
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
# Go to lines 80 and 84
import qm.qua as qua
from qm import QuantumMachinesManager
from qm import SimulationConfig, LoopbackInterface
from scipy.optimize import curve_fit
from OPX1000configuration import config, qop_ip, cluster_name, cloud_username, cloud_password # or OPXplusconfiguration, depending on your hardware    import matplotlib.pyplot as plt
import numpy as np
from qualang_tools.results import progress_counter, fetching_tool
from qualang_tools.plot import interrupt_on_close
from qualang_tools.bakery.randomized_benchmark_c1 import c1_table
from macros import readout_macro
import warnings

warnings.filterwarnings("ignore")

##############################
# Program-specific variables #
##############################

num_of_sequences = 3  # Number of random sequences
n_avg = 1  # Number of averaging loops for each random sequence
max_circuit_depth = 28  # Maximum circuit depth
delta_clifford = 2  #  Play each sequence with a depth step equals to 'delta_clifford - Must be > 1
assert (max_circuit_depth / delta_clifford).is_integer(), "max_circuit_depth / delta_clifford must be an integer."
seed = 345324  # Pseudo-random number generator seed
# Flag to enable state discrimination if the readout has been calibrated (rotated blobs and threshold)
state_discrimination = False
# List of recovery gates from the lookup table
inv_gates = [int(np.where(c1_table[i, :] == 0)[0][0]) for i in range(24)]


###################################
# Helper functions and QUA macros #
###################################
def readout_macro(threshold=None, state=None, I=None, Q=None):
    """
    A macro for performing the readout, with the ability to perform state discrimination.
    If `threshold` is given, the information in the `I` quadrature will be compared against the threshold and `state`
    would be `True` if `I > threshold`.
    Note that it is assumed that the results are rotated such that all the information is in the `I` quadrature.

    :param threshold: Optional. The threshold to compare `I` against.
    :param state: A QUA variable for the state information, only used when a threshold is given.
        Should be of type `bool`. If not given, a new variable will be created
    :param I: A QUA variable for the information in the `I` quadrature. Should be of type `Fixed`. If not given, a new
        variable will be created
    :param Q: A QUA variable for the information in the `Q` quadrature. Should be of type `Fixed`. If not given, a new
        variable will be created
    :return: Three QUA variables populated with the results of the readout: (`state`, `I`, `Q`)
    """
    if I is None:
        I = qua.declare(qua.fixed)
    if Q is None:
        Q = qua.declare(qua.fixed)
    if threshold is not None and state is None:
        state = qua.declare(bool)
    qua.measure(
        "readout",
        "resonator",
        None,
        qua.dual_demod.full("rotated_cos", "rotated_sin", I),
        qua.dual_demod.full("rotated_minus_sin", "rotated_cos", Q),
    )
    if threshold is not None:
        qua.assign(state, I > threshold)
    return state, I, Q

def power_law(power, a, b, p):
    return a * (p**power) + b


def generate_sequence(): 
    cayley = qua.declare(int, value=c1_table.flatten().tolist())
    inv_list = qua.declare(int, value=inv_gates)
    current_state = qua.declare(int)
    step = qua.declare(int)
    sequence = qua.declare(int, size=max_circuit_depth + 1)
    inv_gate = qua.declare(int, size=max_circuit_depth + 1)
    i = qua.declare(int)
    rand = qua.Random(seed=seed) # (1)!

    qua.assign(current_state, 0)
    with qua.for_(i, 0, i < max_circuit_depth, i + 1):
        qua.assign(step, rand.rand_int(24)) # (2)!
        qua.assign(current_state, cayley[current_state * 24 + step])
        qua.assign(sequence[i], step)
        qua.assign(inv_gate[i], inv_list[current_state])

    return sequence, inv_gate


def play_sequence(sequence_list, depth):
    i = qua.declare(int)
    with qua.for_(i, 0, i <= depth, i + 1):
        with qua.switch_(sequence_list[i], unsafe=True):
            with qua.case_(0):
                qua.wait(x180_len // 4, "qubit")
            with qua.case_(1):
                qua.play("x180", "qubit")
            with qua.case_(2):
                qua.play("y180", "qubit")
            with qua.case_(3):
                qua.play("y180", "qubit")
                qua.play("x180", "qubit")
            with qua.case_(4):
                qua.play("x90", "qubit")
                qua.play("y90", "qubit")
            with qua.case_(5):
                qua.play("x90", "qubit")
                qua.play("-y90", "qubit")
            with qua.case_(6):
                qua.play("-x90", "qubit")
                qua.play("y90", "qubit")
            with qua.case_(7):
                qua.play("-x90", "qubit")
                qua.play("-y90", "qubit")
            with qua.case_(8):
                qua.play("y90", "qubit")
                qua.play("x90", "qubit")
            with qua.case_(9):
                qua.play("y90", "qubit")
                qua.play("-x90", "qubit")
            with qua.case_(10):
                qua.play("-y90", "qubit")
                qua.play("x90", "qubit")
            with qua.case_(11):
                qua.play("-y90", "qubit")
                qua.play("-x90", "qubit")
            with qua.case_(12):
                qua.play("x90", "qubit")
            with qua.case_(13):
                qua.play("-x90", "qubit")
            with qua.case_(14):
                qua.play("y90", "qubit")
            with qua.case_(15):
                qua.play("-y90", "qubit")
            with qua.case_(16):
                qua.play("-x90", "qubit")
                qua.play("y90", "qubit")
                qua.play("x90", "qubit")
            with qua.case_(17):
                qua.play("-x90", "qubit")
                qua.play("-y90", "qubit")
                qua.play("x90", "qubit")
            with qua.case_(18):
                qua.play("x180", "qubit")
                qua.play("y90", "qubit")
            with qua.case_(19):
                qua.play("x180", "qubit")
                qua.play("-y90", "qubit")
            with qua.case_(20):
                qua.play("y180", "qubit")
                qua.play("x90", "qubit")
            with qua.case_(21):
                qua.play("y180", "qubit")
                qua.play("-x90", "qubit")
            with qua.case_(22):
                qua.play("x90", "qubit")
                qua.play("y90", "qubit")
                qua.play("x90", "qubit")
            with qua.case_(23):
                qua.play("-x90", "qubit")
                qua.play("y90", "qubit")
                qua.play("-x90", "qubit")


###################
# The QUA program #
###################
with qua.program() as rb:
    depth = qua.declare(int)
    depth_target = qua.declare(int) 
    saved_gate = qua.declare(int)
    m = qua.declare(int)
    n = qua.declare(int)
    I = qua.declare(qua.fixed)
    Q = qua.declare(qua.fixed)
    state = qua.declare(bool)
    m_st = qua.declare_stream()
    if state_discrimination:
        state_st = qua.declare_stream()
    else:
        I_st = qua.declare_stream()
        Q_st = qua.declare_stream()

    with qua.for_(m, 0, m < num_of_sequences, m + 1):
        sequence_list, inv_gate_list = generate_sequence()

        qua.assign(depth_target, 0) 
        with qua.for_(depth, 1, depth <= max_circuit_depth, depth + 1):
            qua.assign(saved_gate, sequence_list[depth])
            qua.assign(sequence_list[depth], inv_gate_list[depth - 1])
            with qua.if_((depth == 1) | (depth == depth_target)):
                with qua.for_(n, 0, n < n_avg, n + 1):
                    play_sequence(sequence_list, depth)
                    qua.align("qubit", "resonator")
                    state, I, Q = readout_macro(threshold=ge_threshold, state=state, I=I, Q=Q)
                    if state_discrimination:
                        qua.save(state, state_st)
                    else:
                        qua.save(I, I_st)
                        qua.save(Q, Q_st)
                    qua.align("resonator", "qubit")
                    qua.wait(thermalization_time * u.ns, "resonator")
                qua.assign(depth_target, depth_target + delta_clifford)
            qua.assign(sequence_list[depth], saved_gate)
        qua.save(m, m_st)

    with qua.stream_processing():
        m_st.save("iteration")
        if state_discrimination:
            state_st.boolean_to_int()
                    .buffer(n_avg)
                    .map(FUNCTIONS.average())
                    .buffer(max_circuit_depth / delta_clifford + 1)
                    .buffer(num_of_sequences)
                    .save("state")
            state_st.boolean_to_int()
                    .buffer(n_avg)
                    .map(FUNCTIONS.average())
                    .buffer(max_circuit_depth / delta_clifford + 1)
                    .average()
                    .save("state_avg")
        else:
            I_st.buffer(n_avg)
                .map(FUNCTIONS.average())
                .buffer(max_circuit_depth / delta_clifford + 1)
                .buffer(num_of_sequences)
                .save("I")
            Q_st.buffer(n_avg)
                .map(FUNCTIONS.average())
                .buffer(max_circuit_depth / delta_clifford + 1)
                .buffer(num_of_sequences)
                .save("Q")
            I_st.buffer(n_avg)
                .map(FUNCTIONS.average())
                .buffer(max_circuit_depth / delta_clifford + 1)
                .average()
                .save("I_avg")
            Q_st.buffer(n_avg)
                .map(FUNCTIONS.average())
                .buffer(max_circuit_depth / delta_clifford + 1)
                .average()
                .save("Q_avg")


#####################################
#  Open Communication with the QOP  #
#####################################
simulate = True
simulation_in_cloud = True

if simulate: 

    if simulation_in_cloud:
        client = QmSaas(email=cloud_username, password=cloud_password)
        instance = client.simulator(QoPVersion.v3_2_0)
        instance.spawn()
        qmm = QuantumMachinesManager(host=instance.host,
                                    port=instance.port,
                                    connection_headers=instance.default_connection_headers)
    else:
        qmm = QuantumMachinesManager(host=qop_ip, 
                                    cluster_name=cluster_name)

    simulation_config = SimulationConfig(
            duration=50_000
        )
    job = qmm.simulate(config, rb, simulation_config)
    job.get_simulated_samples().con1.plot()
    plt.figure()
    if state_discrimination:
        results = fetching_tool(job, data_list=["state_avg", "iteration"])
        state_avg, iteration = results.fetch_all()
        value_avg = state_avg
    else:
        results = fetching_tool(job, data_list=["I_avg", "Q_avg", "iteration"])
        I, Q, iteration = results.fetch_all()
        value_avg = I
    x = np.arange(0, max_circuit_depth + 0.1, delta_clifford)
    x[0] = 1 
    plt.plot(x, value_avg, marker=".")
    plt.xlabel("Number of Clifford gates")
    plt.ylabel("Sequence Fidelity")
    plt.title("Single qubit RB")
    plt.tight_layout()
    plt.show()
else:
    qm = qmm.open_qm(config)
    job = qm.execute(rb)
    if state_discrimination:
        results = fetching_tool(job, data_list=["state_avg", "iteration"], mode="live")
    else:
        results = fetching_tool(job, data_list=["I_avg", "Q_avg", "iteration"], mode="live")
    fig = plt.figure()
    interrupt_on_close(fig, job) 
    x = np.arange(0, max_circuit_depth + 0.1, delta_clifford)
    x[0] = 1 
    while results.is_processing():
        if state_discrimination:
            state_avg, iteration = results.fetch_all()
            value_avg = state_avg
        else:
            I, Q, iteration = results.fetch_all()
            value_avg = I

        progress_counter(iteration, num_of_sequences, start_time=results.get_start_time())
        plt.cla()
        plt.plot(x, value_avg, marker=".")
        plt.xlabel("Number of Clifford gates")
        plt.ylabel("Sequence Fidelity")
        plt.title("Single qubit RB")
        plt.pause(0.1)
    plt.show()

    if state_discrimination:
        results = fetching_tool(job, data_list=["state"])
        state = results.fetch_all()[0]
        value_avg = np.mean(state, axis=0)
        error_avg = np.std(state, axis=0)
    else:
        results = fetching_tool(job, data_list=["I", "Q"])
        I, Q = results.fetch_all()
        value_avg = np.mean(I, axis=0)
        error_avg = np.std(I, axis=0)
    pars, cov = curve_fit(
        f=power_law,
        xdata=x,
        ydata=value_avg,
        p0=[0.5, 0.5, 0.9],
        bounds=(-np.inf, np.inf),
        maxfev=2000,
    )
    stdevs = np.sqrt(np.diag(cov))

    print("#########################")
    print("### Fitted Parameters ###")
    print("#########################")
    print(f"A = {pars[0]:.3} ({stdevs[0]:.1}), B = {pars[1]:.3} ({stdevs[1]:.1}), p = {pars[2]:.3} ({stdevs[2]:.1})")
    print("Covariance Matrix")
    print(cov)

    one_minus_p = 1 - pars[2]
    r_c = one_minus_p * (1 - 1 / 2**1)
    r_g = r_c / 1.875  # 1.875 is the average number of gates in clifford operation
    r_c_std = stdevs[2] * (1 - 1 / 2**1)
    r_g_std = r_c_std / 1.875

    print("#########################")
    print("### Useful Parameters ###")
    print("#########################")
    print(
        f"Error rate: 1-p = {np.format_float_scientific(one_minus_p, precision=2)} ({stdevs[2]:.1})\n"
        f"Clifford set infidelity: r_c = {np.format_float_scientific(r_c, precision=2)} ({r_c_std:.1})\n"
        f"Gate infidelity: r_g = {np.format_float_scientific(r_g, precision=2)}  ({r_g_std:.1})"
    )

    # Plots
    plt.figure()
    plt.errorbar(x, value_avg, yerr=error_avg, marker=".")
    plt.plot(x, power_law(x, *pars), linestyle="--", linewidth=2)
    plt.xlabel("Number of Clifford gates")
    plt.ylabel("Sequence Fidelity")
    plt.title("Single qubit RB")
  1. This line of code initializes the PPU's random number generator. For a detailed overview, please consult our documentation.
  2. Each time the PPU processes this instruction, it selects a new random number, enabling the use of dynamic circuits.

For effective control and decision-making, the PPU supports switch/case functionality. In the RB code, we use this to randomly play the single qubit Clifford set.

rb.py
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
# Go to lines 92-164
import qm.qua as qua
from qm import QuantumMachinesManager
from qm import SimulationConfig, LoopbackInterface
from scipy.optimize import curve_fit
from OPX1000configuration import config, qop_ip, cluster_name, cloud_username, cloud_password # or OPXplusconfiguration, depending on your hardware    import matplotlib.pyplot as plt
import numpy as np
from qualang_tools.results import progress_counter, fetching_tool
from qualang_tools.plot import interrupt_on_close
from qualang_tools.bakery.randomized_benchmark_c1 import c1_table
from macros import readout_macro
import warnings

warnings.filterwarnings("ignore")

##############################
# Program-specific variables #
##############################

num_of_sequences = 3  # Number of random sequences
n_avg = 1  # Number of averaging loops for each random sequence
max_circuit_depth = 28  # Maximum circuit depth
delta_clifford = 2  #  Play each sequence with a depth step equals to 'delta_clifford - Must be > 1
assert (max_circuit_depth / delta_clifford).is_integer(), "max_circuit_depth / delta_clifford must be an integer."
seed = 345324  # Pseudo-random number generator seed
# Flag to enable state discrimination if the readout has been calibrated (rotated blobs and threshold)
state_discrimination = False
# List of recovery gates from the lookup table
inv_gates = [int(np.where(c1_table[i, :] == 0)[0][0]) for i in range(24)]


###################################
# Helper functions and QUA macros #
###################################
def readout_macro(threshold=None, state=None, I=None, Q=None):
    """
    A macro for performing the readout, with the ability to perform state discrimination.
    If `threshold` is given, the information in the `I` quadrature will be compared against the threshold and `state`
    would be `True` if `I > threshold`.
    Note that it is assumed that the results are rotated such that all the information is in the `I` quadrature.

    :param threshold: Optional. The threshold to compare `I` against.
    :param state: A QUA variable for the state information, only used when a threshold is given.
        Should be of type `bool`. If not given, a new variable will be created
    :param I: A QUA variable for the information in the `I` quadrature. Should be of type `Fixed`. If not given, a new
        variable will be created
    :param Q: A QUA variable for the information in the `Q` quadrature. Should be of type `Fixed`. If not given, a new
        variable will be created
    :return: Three QUA variables populated with the results of the readout: (`state`, `I`, `Q`)
    """
    if I is None:
        I = qua.declare(qua.fixed)
    if Q is None:
        Q = qua.declare(qua.fixed)
    if threshold is not None and state is None:
        state = qua.declare(bool)
    qua.measure(
        "readout",
        "resonator",
        None,
        qua.dual_demod.full("rotated_cos", "rotated_sin", I),
        qua.dual_demod.full("rotated_minus_sin", "rotated_cos", Q),
    )
    if threshold is not None:
        qua.assign(state, I > threshold)
    return state, I, Q

def power_law(power, a, b, p):
    return a * (p**power) + b


def generate_sequence(): 
    cayley = qua.declare(int, value=c1_table.flatten().tolist())
    inv_list = qua.declare(int, value=inv_gates)
    current_state = qua.declare(int)
    step = qua.declare(int)
    sequence = qua.declare(int, size=max_circuit_depth + 1)
    inv_gate = qua.declare(int, size=max_circuit_depth + 1)
    i = qua.declare(int)
    rand = qua.Random(seed=seed)

    qua.assign(current_state, 0)
    with qua.for_(i, 0, i < max_circuit_depth, i + 1):
        qua.assign(step, rand.rand_int(24))
        qua.assign(current_state, cayley[current_state * 24 + step])
        qua.assign(sequence[i], step)
        qua.assign(inv_gate[i], inv_list[current_state])

    return sequence, inv_gate


def play_sequence(sequence_list, depth): # (1)!
    i = qua.declare(int)
    with qua.for_(i, 0, i <= depth, i + 1):
        with qua.switch_(sequence_list[i], unsafe=True):
            with qua.case_(0):
                qua.wait(x180_len // 4, "qubit")
            with qua.case_(1):
                qua.play("x180", "qubit")
            with qua.case_(2):
                qua.play("y180", "qubit")
            with qua.case_(3):
                qua.play("y180", "qubit")
                qua.play("x180", "qubit")
            with qua.case_(4):
                qua.play("x90", "qubit")
                qua.play("y90", "qubit")
            with qua.case_(5):
                qua.play("x90", "qubit")
                qua.play("-y90", "qubit")
            with qua.case_(6):
                qua.play("-x90", "qubit")
                qua.play("y90", "qubit")
            with qua.case_(7):
                qua.play("-x90", "qubit")
                qua.play("-y90", "qubit")
            with qua.case_(8):
                qua.play("y90", "qubit")
                qua.play("x90", "qubit")
            with qua.case_(9):
                qua.play("y90", "qubit")
                qua.play("-x90", "qubit")
            with qua.case_(10):
                qua.play("-y90", "qubit")
                qua.play("x90", "qubit")
            with qua.case_(11):
                qua.play("-y90", "qubit")
                qua.play("-x90", "qubit")
            with qua.case_(12):
                qua.play("x90", "qubit")
            with qua.case_(13):
                qua.play("-x90", "qubit")
            with qua.case_(14):
                qua.play("y90", "qubit")
            with qua.case_(15):
                qua.play("-y90", "qubit")
            with qua.case_(16):
                qua.play("-x90", "qubit")
                qua.play("y90", "qubit")
                qua.play("x90", "qubit")
            with qua.case_(17):
                qua.play("-x90", "qubit")
                qua.play("-y90", "qubit")
                qua.play("x90", "qubit")
            with qua.case_(18):
                qua.play("x180", "qubit")
                qua.play("y90", "qubit")
            with qua.case_(19):
                qua.play("x180", "qubit")
                qua.play("-y90", "qubit")
            with qua.case_(20):
                qua.play("y180", "qubit")
                qua.play("x90", "qubit")
            with qua.case_(21):
                qua.play("y180", "qubit")
                qua.play("-x90", "qubit")
            with qua.case_(22):
                qua.play("x90", "qubit")
                qua.play("y90", "qubit")
                qua.play("x90", "qubit")
            with qua.case_(23):
                qua.play("-x90", "qubit")
                qua.play("y90", "qubit")
                qua.play("-x90", "qubit")


###################
# The QUA program #
###################
with qua.program() as rb:
    depth = qua.declare(int)
    depth_target = qua.declare(int) 
    saved_gate = qua.declare(int)
    m = qua.declare(int)
    n = qua.declare(int)
    I = qua.declare(qua.fixed)
    Q = qua.declare(qua.fixed)
    state = qua.declare(bool)
    m_st = qua.declare_stream()
    if state_discrimination:
        state_st = qua.declare_stream()
    else:
        I_st = qua.declare_stream()
        Q_st = qua.declare_stream()

    with qua.for_(m, 0, m < num_of_sequences, m + 1):
        sequence_list, inv_gate_list = generate_sequence()

        qua.assign(depth_target, 0) 
        with qua.for_(depth, 1, depth <= max_circuit_depth, depth + 1):
            qua.assign(saved_gate, sequence_list[depth])
            qua.assign(sequence_list[depth], inv_gate_list[depth - 1])
            with qua.if_((depth == 1) | (depth == depth_target)):
                with qua.for_(n, 0, n < n_avg, n + 1):
                    play_sequence(sequence_list, depth)
                    qua.align("qubit", "resonator")
                    state, I, Q = readout_macro(threshold=ge_threshold, state=state, I=I, Q=Q)
                    if state_discrimination:
                        qua.save(state, state_st)
                    else:
                        qua.save(I, I_st)
                        qua.save(Q, Q_st)
                    qua.align("resonator", "qubit")
                    qua.wait(thermalization_time * u.ns, "resonator")
                qua.assign(depth_target, depth_target + delta_clifford)
            qua.assign(sequence_list[depth], saved_gate)
        qua.save(m, m_st)

    with qua.stream_processing():
        m_st.save("iteration")
        if state_discrimination:
            state_st.boolean_to_int()
                    .buffer(n_avg)
                    .map(FUNCTIONS.average())
                    .buffer(max_circuit_depth / delta_clifford + 1)
                    .buffer(num_of_sequences)
                    .save("state")
            state_st.boolean_to_int()
                    .buffer(n_avg)
                    .map(FUNCTIONS.average())
                    .buffer(max_circuit_depth / delta_clifford + 1)
                    .average()
                    .save("state_avg")
        else:
            I_st.buffer(n_avg)
                .map(FUNCTIONS.average())
                .buffer(max_circuit_depth / delta_clifford + 1)
                .buffer(num_of_sequences)
                .save("I")
            Q_st.buffer(n_avg)
                .map(FUNCTIONS.average())
                .buffer(max_circuit_depth / delta_clifford + 1)
                .buffer(num_of_sequences)
                .save("Q")
            I_st.buffer(n_avg)
                .map(FUNCTIONS.average())
                .buffer(max_circuit_depth / delta_clifford + 1)
                .average()
                .save("I_avg")
            Q_st.buffer(n_avg)
                .map(FUNCTIONS.average())
                .buffer(max_circuit_depth / delta_clifford + 1)
                .average()
                .save("Q_avg")


#####################################
#  Open Communication with the QOP  #
#####################################
simulate = True
simulation_in_cloud = True

if simulate: 

    if simulation_in_cloud:
        client = QmSaas(email=cloud_username, password=cloud_password)
        instance = client.simulator(QoPVersion.v3_2_0)
        instance.spawn()
        qmm = QuantumMachinesManager(host=instance.host,
                                    port=instance.port,
                                    connection_headers=instance.default_connection_headers)
    else:
        qmm = QuantumMachinesManager(host=qop_ip, 
                                    cluster_name=cluster_name)

    simulation_config = SimulationConfig(
            duration=50_000
        )
    job = qmm.simulate(config, rb, simulation_config)
    job.get_simulated_samples().con1.plot()
    plt.figure()
    if state_discrimination:
        results = fetching_tool(job, data_list=["state_avg", "iteration"])
        state_avg, iteration = results.fetch_all()
        value_avg = state_avg
    else:
        results = fetching_tool(job, data_list=["I_avg", "Q_avg", "iteration"])
        I, Q, iteration = results.fetch_all()
        value_avg = I
    x = np.arange(0, max_circuit_depth + 0.1, delta_clifford)
    x[0] = 1 
    plt.plot(x, value_avg, marker=".")
    plt.xlabel("Number of Clifford gates")
    plt.ylabel("Sequence Fidelity")
    plt.title("Single qubit RB")
    plt.tight_layout()
    plt.show()
else:
    qm = qmm.open_qm(config)
    job = qm.execute(rb)
    if state_discrimination:
        results = fetching_tool(job, data_list=["state_avg", "iteration"], mode="live")
    else:
        results = fetching_tool(job, data_list=["I_avg", "Q_avg", "iteration"], mode="live")
    fig = plt.figure()
    interrupt_on_close(fig, job) 
    x = np.arange(0, max_circuit_depth + 0.1, delta_clifford)
    x[0] = 1 
    while results.is_processing():
        if state_discrimination:
            state_avg, iteration = results.fetch_all()
            value_avg = state_avg
        else:
            I, Q, iteration = results.fetch_all()
            value_avg = I

        progress_counter(iteration, num_of_sequences, start_time=results.get_start_time())
        plt.cla()
        plt.plot(x, value_avg, marker=".")
        plt.xlabel("Number of Clifford gates")
        plt.ylabel("Sequence Fidelity")
        plt.title("Single qubit RB")
        plt.pause(0.1)
    plt.show()

    if state_discrimination:
        results = fetching_tool(job, data_list=["state"])
        state = results.fetch_all()[0]
        value_avg = np.mean(state, axis=0)
        error_avg = np.std(state, axis=0)
    else:
        results = fetching_tool(job, data_list=["I", "Q"])
        I, Q = results.fetch_all()
        value_avg = np.mean(I, axis=0)
        error_avg = np.std(I, axis=0)
    pars, cov = curve_fit(
        f=power_law,
        xdata=x,
        ydata=value_avg,
        p0=[0.5, 0.5, 0.9],
        bounds=(-np.inf, np.inf),
        maxfev=2000,
    )
    stdevs = np.sqrt(np.diag(cov))

    print("#########################")
    print("### Fitted Parameters ###")
    print("#########################")
    print(f"A = {pars[0]:.3} ({stdevs[0]:.1}), B = {pars[1]:.3} ({stdevs[1]:.1}), p = {pars[2]:.3} ({stdevs[2]:.1})")
    print("Covariance Matrix")
    print(cov)

    one_minus_p = 1 - pars[2]
    r_c = one_minus_p * (1 - 1 / 2**1)
    r_g = r_c / 1.875  # 1.875 is the average number of gates in clifford operation
    r_c_std = stdevs[2] * (1 - 1 / 2**1)
    r_g_std = r_c_std / 1.875

    print("#########################")
    print("### Useful Parameters ###")
    print("#########################")
    print(
        f"Error rate: 1-p = {np.format_float_scientific(one_minus_p, precision=2)} ({stdevs[2]:.1})\n"
        f"Clifford set infidelity: r_c = {np.format_float_scientific(r_c, precision=2)} ({r_c_std:.1})\n"
        f"Gate infidelity: r_g = {np.format_float_scientific(r_g, precision=2)}  ({r_g_std:.1})"
    )

    # Plots
    plt.figure()
    plt.errorbar(x, value_avg, yerr=error_avg, marker=".")
    plt.plot(x, power_law(x, *pars), linestyle="--", linewidth=2)
    plt.xlabel("Number of Clifford gates")
    plt.ylabel("Sequence Fidelity")
    plt.title("Single qubit RB")
  1. Below you see an example of switch/case in action, and as you can see the input to the switch evaluation can be a QUA variable.

Test your knowledge

  1. Follow the steps in rb.py and using the c1_table, generate random Clifford circuits and its corresponding inverse gates with pen and paper for the following conditions:

    1. For a depth of 3, and the random number in each iteration to be 7.
    2. For a depth of 13, and the random number list for the sequential iterations to be [3, 17, 22, 8, 0, 14, 19, 6, 21, 11, 15, 5, 20]
    3. Given that we declared the switch/case to be unsafe=True, What would occur if you remove with case_(5) statement in the play_sequence() macro?

The End is the Beginning

After mastering both the basic and advanced concepts of the QM Platform, you're well-prepared to face your experimental challenges. Continue on your journey!