How to automate calibration of control hardware

Calibrate RF control channels for maximum pulse performance

The Q-CTRL Python Package outputs various control solutions that can be used to define operations in real control hardware. Here we will show you how to automate the process of calibrating the RF or microwave channels used to implement controls for real quantum hardware in order to achieve maximum-fidelity operations. The calibration process allows translation of the pulse that you define at the Hamiltonian level to a useful control applied to a device, and also provides a means to account for small hardware imperfections.

In this notebook we outline the general approach, demonstrate the calibration procedure using IBM hardware, and experimentally validate the procedure's performance using complex optimized gates formatted for Qiskit. For concreteness we focus on examples calibrating single-channel RF drives. However these approaches can be generalized for the calibration of any control signal used in any quantum hardware or control electronics.

IQ channel calibration formulation and objective

Most control solutions involve programming a complex time-varying control that can be written using a traditional $IQ$ decomposition familiar in RF engineering:

\begin{align*} \gamma(t) &= \Omega(t) e^{i\phi(t)} = I(t) + i Q(t) \end{align*}

using the $IQ$ decomposition

\begin{align*} I = \Omega \cos(\phi),\\ Q = \Omega \sin(\phi),\\ |\gamma|\le\Omega_\text{max}. \end{align*}

The central calibration task is to determine the relationship between these dimensionful units and the dimensionless variables employed in programming a hardware device, $A = A_I + i A_Q$. For instance on the IBM hardware these variables take values $A_{I,Q} \in [-1,1]$, although these may alternatively be e.g. integers representing a 16-bit number. This calibration must be performed with sufficient fidelity to allow the implementation of high-fidelity controls which faithfully represent the intended Hamiltonian terms. We break this into two steps below.

Summary workflow for IQ calibration

1. Coarse IQ channel calibration

Coarsely calibrate the drive amplitude by measuring any recognizable system response (such as Rabi oscillations) using square pulses with variable amplitude or length. Repeat for different square-pulse settings to produce a coarse map of $(A_I, A_Q) \longleftrightarrow (I, Q)$.

This procedure can identify nonlinear regimes that may required pre-distortion and also mitigates the $\mathrm{mod} (2\pi)$ calibration ambiguity that inhibits direct automation in systems exhibiting periodicity (e.g. Rabi oscillations). Performing an initial coarse calibration of system response to a single pulse stimulus permits appropriate bounding of parameter spaces for automation.

The workflow for this coarse calibration is as follows:

  • Run Rabi experiments for a range of pulse amplitudes $A_I$.
  • Fit the resulting oscillations to obtain the Rabi rates $I$ for each $A_I$.
  • Interpolate the results to establish the mapping $A_I \longleftrightarrow I$.

A detailed implementation of this procedure on the IBM Quantum hardware is presented in the Worked example section below.

2. Fine IQ channel calibration

Use automated closed-loop optimization (or scripting) to determine two additional scaling factors $(S_{amp}, S_{rel})$

\begin{align*} \gamma(t) = S_{amp} (S_{rel} A_I(t) + i A_Q(t)). \end{align*}

The $S_{rel}$ factor allows you to account for relative differences between the designed and the actual values of $I$ and $Q$ reaching the device, while $S_{amp}$ accounts for an absolute amplitude scaling. Starting values of these fine-tuning scaling factors should be $\approx1$ as $S_{amp} = S_{rel} = 1$ corresponds to the values obtained from the previous coarse calibration.

In the calibration procedure, error amplification may be employed in order to accurately discern small differences in the impact of $S_{amp}$ and $S_{rel}$. This is easily achieved by employing gate repetition - for simplicity you should choose multiples of the target gate repetition to return a net operation the same as a single operation.

Determining $(S_{amp}, S_{rel})$ is in general sufficient to produce high-fidelity implementation of an arbitrary control operation. We note that further calibration may require system identification to determine unknown signal distortion via transmission lines.

Worked example: Calibration of single-qubit gates on IBM hardware

In the following we will perform these steps using IBM hardware and demonstrate a series of verification experiments to illustrate the benefits of this approach.

Our target will be the efficient calibration of the control lines used to implement a single-qubit gate. In the first section of this example, you'll learn how to perform coarse calibration of Rabi rates on a hardware device - in this case from IBM. With this initial information we will then show how to use automated closed-loop optimization in order to perform fine calibration. As an alternative we will also present a simpler but less efficient approach based on scripting to perform the same fine calibration. We will conclude with pulse-repetition experiments that illustrate the efficacy of the calibration approach outlined above.

Some cells in this notebook require an account with IBM-Q to execute correctly. If you want to run them, please go to the IBM-Q experience to set up an account.

Coarse calibration of Rabi rates

We begin by identifying the relationship $(A_I, A_Q) \longleftrightarrow (I, Q)$ by performing experiments on IBM hardware, accessed using Qiskit pulse (analog-layer programming). In all circumstances Qiskit-specific code can be replaced with your own hardware-programming scripts.

# Imports and initialization
import jsonpickle
import matplotlib.gridspec as gridspec
import matplotlib.pyplot as plt

import numpy as np
from qctrlvisualizer import get_qctrl_style, plot_controls
from scipy import interpolate
from scipy.optimize import curve_fit

# Q-CTRL imports
from qctrl import Qctrl

# Starting a session with the API
qctrl = Qctrl()

# Choose to run experiments or to use saved data
use_saved_data = True

# Plotting parameters
plt.style.use(get_qctrl_style())
markers = {"x": "x", "y": "s", "z": "o"}
lines = {"x": "--", "y": "-.", "z": "-"}

# Definition of operators and functions to be calibrated
sigma_z = np.array([[1.0, 0.0], [0.0, -1.0]], dtype=complex)
sigma_x = np.array([[0.0, 1.0], [1.0, 0.0]], dtype=complex)
sigma_y = np.array([[0.0, -1.0j], [1.0j, 0.0]], dtype=complex)
X90_gate = np.array([[1.0, -1j], [-1j, 1.0]], dtype=complex) / np.sqrt(2)
bloch_basis = ["x", "y", "z"]


# Initialize Qiskit parameters and IBM hardware
import warnings

warnings.simplefilter("ignore")
use_IBM = False

if use_IBM == True:
    # IBM-Q imports
    import qiskit as qs
    import qiskit.pulse as pulse
    from qiskit.pulse.library import Waveform
    from qiskit import IBMQ, execute
    from qiskit.compiler import assemble
    from qiskit.pulse import Acquire, Play, Schedule
    from qiskit.tools.jupyter import *
    from qiskit.tools.monitor import job_monitor
    from qiskit.ignis.verification import marginal_counts
    from qiskit.ignis.mitigation.measurement import (
        complete_meas_cal,
        CompleteMeasFitter,
    )

    # IBM credentials and backend selection
    IBMQ.enable_account("your IBM token")
    provider = IBMQ.get_provider(
        hub="your hub", group="your group", project="your project"
    )
    backend = provider.get_backend("ibmq_jakarta")
    backend_defaults = backend.defaults()
    backend_config = backend.configuration()

    # Backend properties
    dt = backend_config.dt
    print(f"Hardware sampling time: {dt/1e-9} ns")

    qubit_freq_est = []
    for qubit in backend_config.meas_map[0]:
        qubit_freq_est.append(backend_defaults.qubit_freq_est[qubit])
        print(f"Qubit [{qubit}] frequency estimate: {qubit_freq_est[qubit]/1e9} GHz")
# define functions used in performance verification
def save_var(file_name, var):
    # saves a single var to a file using jsonpickle
    f = open(file_name, "w+")
    to_write = jsonpickle.encode(var)
    f.write(to_write)
    f.close()


def load_var(file_name):
    # returns a var from a json file
    f = open(file_name, "r+")
    encoded = f.read()
    decoded = jsonpickle.decode(encoded)
    f.close()
    return decoded


def fit_function_bounds(x_values, y_values, function, bound_values):
    fitparams, conv = curve_fit(function, x_values, y_values, bounds=bound_values)
    y_fit = function(x_values, *fitparams)
    return fitparams, y_fit


def movingaverage(x, w):
    return np.convolve(x, np.ones(w), "valid") / w
# Setting up calibration experiments
qubit = 0
dt = 2 / 9 * 1e-9
num_shots_per_point = 1024

# Variable square pulse amplitudes
pulse_amp_array = np.linspace(0.05, 0.2, 7)

# Variable square pulse durations
pulse_times = np.array(
    [4 + np.arange(0, int(3.6 / (amplitude)), 1) for amplitude in pulse_amp_array]
)
pulse_times = pulse_times * 16


# Experimental setup
if use_saved_data == False:
    backend.properties(refresh=True)
    qubit_frequency_updated = backend.properties().qubit_property(qubit, "frequency")[0]

    meas_map_idx = None
    for i, measure_group in enumerate(backend_config.meas_map):
        if qubit in measure_group:
            meas_map_idx = i
            break
    assert meas_map_idx is not None, f"Couldn't find qubit {qubit} in the meas_map!"

    inst_sched_map = backend_defaults.instruction_schedule_map
    measure_schedule = inst_sched_map.get("measure", qubits=[qubit])
    drive_chan = pulse.DriveChannel(qubit)

    rabi_programs_dic_I = {}
    for idx, pulse_amplitude in enumerate(pulse_amp_array):
        rabi_schedules_I = []
        for duration_pulse in pulse_times[idx]:
            drive_pulse = qs.pulse.library.gaussian_square(
                duration=duration_pulse,
                sigma=1,
                amp=pulse_amplitude,
                risefall=1,
                name=f"square_pulse_{duration_pulse}",
            )
            schedule = pulse.Schedule(name=str(duration_pulse))
            schedule |= (
                Play(drive_pulse, pulse.DriveChannel(qubit)) << schedule.duration
            )
            schedule += measure_schedule << schedule.duration
            rabi_schedules_I.append(schedule)
        rabi_experiment_program_I = assemble(
            rabi_schedules_I,
            backend=backend,
            meas_level=2,
            meas_return="single",
            shots=num_shots_per_point,
            schedule_los=[{drive_chan: qubit_frequency_updated}]
            * len(pulse_times[idx]),
        )
        rabi_programs_dic_I[pulse_amplitude] = rabi_experiment_program_I

    # Running calibration experiments
    rabi_calibration_exp_I = []
    rabi_oscillations_results = []
    for idx, pulse_amplitude in enumerate(pulse_amp_array):
        job = backend.run(rabi_programs_dic_I[pulse_amplitude])
        job_monitor(job)
        rabi_results = job.result(timeout=120)
        rabi_values = []
        time_array = pulse_times[idx] * dt
        for time_idx in pulse_times[idx]:
            counts = rabi_results.get_counts(str(time_idx))
            excited_pop = 0
            for bits, count in counts.items():
                excited_pop += count if bits[::-1][qubit] == "1" else 0
            rabi_values.append(excited_pop / num_shots_per_point)

        rabi_oscillations_results.append(rabi_values)
        fit_parameters, y_fit = fit_function_bounds(
            time_array,
            rabi_values,
            lambda x, A, rabi_freq, phi: A
            * np.cos(2 * np.pi * rabi_freq * x + phi) ** 2,
            (
                [0.8, np.abs(pulse_amplitude * 8 * 1e7), -4],
                [1, np.abs(pulse_amplitude * 11 * 1e7), 4],
            ),
        )

        rabi_calibration_exp_I.append(fit_parameters[1])

    save_var(
        "resources/superconducting-qubits-pulse-calibration/rabi_calibration_jakarta_qubit_0",
        rabi_calibration_exp_I,
    )
    save_var(
        "resources/superconducting-qubits-pulse-calibration/fit_parameters",
        fit_parameters,
    )
    save_var(
        "resources/superconducting-qubits-pulse-calibration/rabi_values", rabi_values
    )
else:
    rabi_calibration_exp_I = load_var(
        "resources/superconducting-qubits-pulse-calibration/rabi_calibration_jakarta_qubit_0"
    )
    fit_parameters = load_var(
        "resources/superconducting-qubits-pulse-calibration/fit_parameters"
    )
    rabi_values = load_var(
        "resources/superconducting-qubits-pulse-calibration/rabi_values"
    )

time_array = pulse_times[-1] * dt
print("Drive amplitude:", pulse_amp_array[-1])
print("Fitted Rabi frequency [Hz]:", fit_parameters[1])

fig = plt.figure(figsize=(10, 5))
fig.suptitle("Exemplary Rabi oscillation data with fitting")
plt.xlabel("Time (s)")
plt.ylabel("Measured signal (a.u.)")
plt.scatter(time_array, np.real(rabi_values), color="black")
plt.xlim(0, time_array[-1])
plt.ylim(0, 1)
plot_times = np.linspace(0, time_array[-1], 100)
plt.plot(
    plot_times,
    fit_parameters[0]
    * np.cos(2 * np.pi * fit_parameters[1] * plot_times + fit_parameters[2]) ** 2,
)
plt.show()
Drive amplitude: 0.2
Fitted Rabi frequency [Hz]: 16000000.0000049

This figure shows the result of a Rabi experiment for an amplitude $A_I = 0.2$ (circles) and the correponding cosine square fit (line).

You can now collect all the $I$ values from the previous experiment and establish the correspondence between $I$ and $A_I$. This is performed in the cell below.

amplitude_interpolated_list = np.linspace(-0.2, 0.2, 100)
pulse_amp_array = np.concatenate((-pulse_amp_array[::-1], pulse_amp_array))
rabi_calibration_exp_I = np.concatenate(
    (-np.asarray(rabi_calibration_exp_I[::-1]), np.asarray(rabi_calibration_exp_I))
)
f_amp_to_rabi = interpolate.interp1d(pulse_amp_array, rabi_calibration_exp_I)
rabi_interpolated_exp_I = f_amp_to_rabi(amplitude_interpolated_list)

f_rabi_to_amp = interpolate.interp1d(
    rabi_interpolated_exp_I, amplitude_interpolated_list
)

fig = plt.figure(figsize=(10, 5))
fig.suptitle("IBMQ Jakarta: Rabi rates calibration")
plt.xlabel("Hardware input amplitude")
plt.ylabel("Rabi rate (MHz)")
plt.scatter(pulse_amp_array, rabi_calibration_exp_I * 1e-6)
plt.plot(amplitude_interpolated_list, rabi_interpolated_exp_I * 1e-6)
plt.axvline(0, color="black", linestyle="dashed")
plt.axhline(0, color="black", linestyle="dashed")
plt.show()

The figure shows the mapping between the Rabi rates and the dimensionless hardware input. The calibration was performed in the interval $[-0.2, 0.2]$, and the observed relationship is approximately linear. With this information we have completed the coarse calibration of Rabi rates for this specific complex drive channel.

Fine IQ calibration using an autonomous agent

We use automated closed-loop optimization to find the calibration parameters $S_{amp}, S_{rel}$ in the pulse parametrization $\gamma(t) = S_{amp} (S_{rel} A_I(t) + i A_Q(t))$ automatically. In each iteration of the procedure, the optimizer will pick a set of calibration parameters, transform the pulse, send it to the device, and receive in response measurements of the system (here measured infidelity under gate repetition) to inform the next actions of the autonomous agent. As more iterations are performed, the optimizable parameters will converge to the values that minimize the infidelity.

In what follows, the optimization bounds are set to be $S_{amp} \in [0.8,1.2]$ and $S_{rel} \in [0.8,1.2]$, that is, within $20 \%$ of the values of the initial coarse calibration step.

In each iteration of the calibration loop, you'll need to extract a measurement to be returned to the optimizer. The cost function may be created in a manner that best suits your hardware; we recommend measuring gate infidelity under different numbers of repetition $N$ and then averaging the fidelities for each value of gate repetition. Using repetitions amplifies the error in each gate, and averaging over different repetition numbers mitigates pathologies arising for specific experimental conditions. For the $X_{\pi/2}$ pulses, you can use $N = (1 + 4n)$ repetitions with $n \in \{0, 1 \cdots \}$ so that the Bloch vector always ends up in the same point if the pulse is perfect.

The functions in the cell below show how to do this using Qiskit. If you want to deploy this calibration in your own experiment, you'll need to modify these functions accordingly.

error_mitigation = True


def run_ibm_experiments(control, calibration_pairs, repetitions, qubit):

    I_values = np.array([segment["value"] for segment in control["I"]])
    Q_values = np.array([segment["value"] for segment in control["Q"]])

    schedules = []
    for scal in calibration_pairs:
        scaled_waveform = (I_values * scal[0] + 1j * Q_values) * scal[1]
        control_pulse = Play(Waveform(scaled_waveform), DriveChannel(qubit))
        for meas_basis in bloch_basis:
            for repetition in repetitions:
                schedule = pulse.Schedule(
                    name=f"pulse_scal_{scal}_reps_"
                    + str(repetition)
                    + "_basis_"
                    + str(meas_basis)
                )
                for reps in np.arange(repetition):
                    schedule += control_pulse
                if meas_basis == "x":
                    schedule += inst_sched_map.get("u2", [qubit], P0=0.0, P1=np.pi)
                    schedule += measure_schedule << schedule.duration
                if meas_basis == "y":
                    schedule += inst_sched_map.get("u2", [qubit], P0=0.0, P1=np.pi / 2)
                    schedule += measure_schedule << schedule.duration
                if meas_basis == "z":
                    schedule |= measure_schedule << schedule.duration
                schedules.append(schedule)

    # Running the jobs
    job = execute(schedules, backend, shots=num_shots_per_point)
    print("Job id:", job.job_id())
    job_monitor(job)
    exp_results = job.result(timeout=120)

    return exp_results


def extract_results(results, calibration_pairs):

    idx = 0
    pairs_results = []
    if mitigation:
        results = measurementFilter.apply(results)

    for scaling_pair in calibration_pairs:
        basis_results = []
        for meas_basis in bloch_basis:
            rep_results = []
            for number_of_repetitions in repetitions:
                count = marginal_counts(results.get_counts()[idx], meas_qubits=[0])
                excited_pop = 0
                idx += 1
                excited_pop += count["1"]
                rep_results.append(excited_pop / num_shots_per_point)
            basis_results.append(rep_results)
        pairs_results.append(1 - 2 * np.asarray(basis_results))
    calibration_data = np.asarray(pairs_results)

    # Infidelities for repetition experiments
    target_bloch = np.array([0, -1, 0])
    bloch_vectors = np.array(
        [
            np.array([calibration_data[idx, :, rep] for rep in range(len(repetitions))])
            for idx in range(len(calibration_pairs))
        ]
    )
    bloch_infidelity = 1 - ((1 + np.dot(bloch_vectors, target_bloch)) / 2) ** 2

    cost_infidelity = np.average(bloch_infidelity, axis=1)
    print("Best cost in this step:", np.min(cost_infidelity))

    return cost_infidelity


def conductErrorMitigation(shots):
    """Measure the state preparation and measurement error"""

    cal_circuits, state_labels = qs.ignis.mitigation.measurement.complete_meas_cal(
        qubit_list=np.arange(7), circlabel="MeasErrMitCal"
    )

    cal_job = qs.execute(cal_circuits, backend=backend, shots=shots)
    qs.tools.monitor.job_monitor(cal_job)
    cal_results = cal_job.result()

    meas_fitter = qs.ignis.mitigation.measurement.CompleteMeasFitter(
        cal_results, state_labels
    )

    # Set the object's measurementFilter here
    meas_filter = meas_fitter.filter

    ErrorMitigationComplete = True
    return meas_filter


if use_saved_data == False and error_mitigation:
    SPAMshots = 1024 * 4
    measurementFilter = conductErrorMitigation(shots=SPAMshots)

Here, you'll calibrate a specific control pulse using the coarse calibration information performed above. The cell below imports a specific pulse to calibrate and plots the waveform. In the Application Note "Designing noise-robust single-qubit gates for IBM Qiskit" you'll find a detailed explanation on how to construct robust pulses for superconducting qubit systems.

# Import pulses
gate_schemes = load_var(
    "resources/superconducting-qubits-pulse-calibration/pulse_amplitude_robust"
)

# Extract pulse properties
number_of_segments = {}
duration_int = {}
duration = {}
segment_scale = {}
waveform = {}
for scheme_name, control in gate_schemes.items():
    number_of_segments[scheme_name] = len(control["I"])
    duration[scheme_name] = sum(values["duration"] for values in control["I"])
    duration_int[scheme_name] = int(round(duration[scheme_name] / dt))
    segment_scale[scheme_name] = (
        duration_int[scheme_name] / number_of_segments[scheme_name]
    )
    I_values = [segment["value"] / 2 / np.pi for segment in control["I"]]
    Q_values = [segment["value"] / 2 / np.pi for segment in control["Q"]]
    A_I_values = np.repeat(f_rabi_to_amp(I_values), segment_scale[scheme_name])
    A_Q_values = np.repeat(f_rabi_to_amp(Q_values), segment_scale[scheme_name])
    waveform[scheme_name] = A_I_values + 1j * A_Q_values

# Plot pulses
for scheme_name, gate in gate_schemes.items():
    fig = plt.figure()
    fig.suptitle(scheme_name)
    plot_controls(fig, gate)
    plt.show()

The above plots display I and Q components of the microwave drive on the qubit for an amplitude robust pulse that performs a $\pi/2$ rotation around the $x$ axis (an $X_{\pi/2}$ gate). Note that the values of the pulse waveform are given in natural units of MHz.

# Invoke the pulse described above in the calibration procedure
pulse_to_calibrate = {
    "I": [
        {"duration": d, "value": v}
        for d, v in zip(
            dt * np.ones(waveform["Amplitude robust"].shape[0]),
            np.real(waveform["Amplitude robust"]),
        )
    ],
    "Q": [
        {"duration": d, "value": v}
        for d, v in zip(
            dt * np.ones(waveform["Amplitude robust"].shape[0]),
            np.imag(waveform["Amplitude robust"]),
        )
    ],
}

The next cell sets up the automated optimization and runs the calibration loop. Here you will directly set the optimization cost as the average infidelity for multiple repetitions, streamlining the calibration process by requiring less human intervention.

To minimize the number of API calls to the IBM Quantum device, the calibration procedure will make use of the batching feature of the Q-CTRL Python package (see "How to automate closed-loop hardware optimization" for more information). By setting test_point_count to 5, you'll send five independent jobs under a single request in the queue. In this way, the optimizer will use the result from five different sets of calibration parameters in each iteration.

if use_saved_data:
    best_cost_array = load_var(
        "resources/superconducting-qubits-pulse-calibration/best_cost_array_jakarta_June_8"
    )
    best_param_array = load_var(
        "resources/superconducting-qubits-pulse-calibration/best_control_array_jakarta_June_8"
    )
    print(f"Infidelity: {best_cost_array[-1]}")
    print(f"Calibration parameters: {best_param_array[-1]}")
else:
    repetitions = np.arange(5, 42, 12)
    relative_bounds = (0.8, 1.2)
    amplitude_bounds = (0.8, 1.2)
    test_point_count = 5

    calibration_result = {}
    best_cost_array = []
    best_param_array = []

    # Create random pairs for initial batch
    rng = np.random.default_rng()
    parameter_set = [
        [
            rng.uniform(relative_bounds[0], relative_bounds[1]),
            rng.uniform(amplitude_bounds[0], amplitude_bounds[1]),
        ]
        for idx in np.arange(test_point_count)
    ]

    # Run the experiment with the initial random calibration parameters
    exp_results = run_ibm_experiments(
        pulse_to_calibrate, parameter_set, repetitions, qubit
    )
    cost_results = extract_results(exp_results, parameter_set)

    # OPTIMIZATION SETUP
    # Define initialization object for the automated closed-loop optimization
    bound = qctrl.types.closed_loop_optimization_step.BoxConstraint(
        lower_bound=relative_bounds[0],
        upper_bound=relative_bounds[1],
    )
    bound_2 = qctrl.types.closed_loop_optimization_step.BoxConstraint(
        lower_bound=amplitude_bounds[0],
        upper_bound=amplitude_bounds[1],
    )
    initializer_GP = (
        qctrl.types.closed_loop_optimization_step.GaussianProcessInitializer(
            bounds=[bound, bound_2],
            rng_seed=0,
        )
    )

    # Define state object for the closed-loop optimization
    optimizer = qctrl.types.closed_loop_optimization_step.Optimizer(
        gaussian_process_initializer=initializer_GP,
    )

    best_cost, best_controls = min(
        zip(cost_results, parameter_set), key=lambda params: params[0]
    )
    optimization_count = 0

    best_cost_array.append(best_cost)
    best_param_array.append(best_controls)

    # Run the optimization loop until the cost is sufficiently small or when it reaches 10 steps
    idx = 0
    while best_cost > 0.005 and idx <= 10:
        # Print the current best cost.
        optimization_steps = (
            "optimization step" if optimization_count == 1 else "optimization steps"
        )
        print(
            f"Best infidelity after {optimization_count} Boulder Opal {optimization_steps}: {best_cost}"
        )

        # Organize the experiment results into the proper input format
        results = [
            qctrl.types.closed_loop_optimization_step.CostFunctionResult(
                parameters=list(parameters),
                cost=cost,
                cost_uncertainty=0.1,
            )
            for parameters, cost in zip(parameter_set, cost_results)
        ]

        # Call the automated closed-loop optimizer and obtain the next set of test points
        optimization_result = qctrl.functions.calculate_closed_loop_optimization_step(
            optimizer=optimizer,
            results=results,
            test_point_count=test_point_count,
        )
        optimization_count += 1

        # Organize the data returned by the automated closed-loop optimizer
        parameter_set = np.array(
            [test_point.parameters for test_point in optimization_result.test_points]
        )
        optimizer = qctrl.types.closed_loop_optimization_step.Optimizer(
            state=optimization_result.state
        )

        # Obtain experiment results that the automated closed-loop optimizer requested
        exp_results = run_ibm_experiments(
            pulse_to_calibrate, parameter_set, repetitions, qubit
        )
        cost_results = extract_results(exp_results, parameter_set)

        params_array.append(parameter_set)

        # Record the best results after this round of experiments
        cost, controls = min(
            zip(cost_results, parameter_set), key=lambda params: params[0]
        )

        if cost < best_cost:
            best_cost = cost
            best_controls = controls
        idx += 1
        best_cost_array.append(best_cost)
        best_param_array.append(best_controls)
        print("Best parameters:", best_controls)

    # Print final best cost.
    print(f"Infidelity: {best_cost}")
    print(f"Calibration parameters: {best_controls}")
    calibration_result[qubit] = best_controls
Infidelity: 0.00827659524514815
Calibration parameters: [1.02161937 1.04297263]
# Plotting
fig = plt.figure(figsize=(10, 5))
fig.suptitle("Convergence of automated pulse calibration")
plt.scatter(np.arange(12), best_cost_array)
plt.plot(best_cost_array)
plt.xlabel("Optimization step")
plt.ylabel("Cost (average infidelity)")

plt.show()

Gate infidelity averaged for multiple repetitions of the pulses (5, 17, 29, and 41, in this case) as a function of the number of optimization steps. The calibration converges after a couple of steps, with the cost below $1\%$, after readout error mitigation. The final calibration parameters are $(S_{amp}, S_{rel}) = (1.043, 1.022)$.

Fine IQ calibration using scripted parametric scans over $S_{rel}$ and $S_{amp}$

An alternative approach to performing fine calibration involves the execution of scripted scans over the parameters $(S_{amp}, S_{rel})$. This approach may be used as a manual diagnostic or in cases where automated optimization is not accessible. The procedure is described below.

  1. Using two sequences with different numbers of pulse repetitions ($N_1$ and $N_2$), perform experiments where you measure gate infidelity as a function of $S_{rel}\approx 1$, scanning over a narrow range.

  2. Find the value of the scaling factor, $S_{rel}$, that minimizes the infidelity for both repetition numbers. You could be tempted to use a single repetition number $N$ to speed up your calibration but this is not advisable. Note that, if you pick $N$ that is too small, there will be little variation of the infidelity with $S_{rel}$ and your calibration would be inaccurate. If you select $N$ that is too large, you may see multiple infidelity minima as small changes in $S_{rel}$ can lead to full rotations in the Bloch sphere. In this case, you may be selecting the wrong scaling factor. Using two values allows you to find a single scaling that best calibrates your pulse.

  3. Repeat step 2, now fixing the value of $S_{rel}$ found to minimize infidelity in the previous step and scaning $S_{amp}$. Again, find the scaling factor that minimizes the infidelity for both pulse repetition numbers.

The following cell executes steps 1 and 2 above for the IBM device using an automated routine to identify the appropriate value of $S_{rel}$. Here we select $N_1 = 49$ and $N_2 =101$.

repetitions = np.array([49, 101])
scalings = np.linspace(0.95, 1.05, 50)

if use_saved_data == False:
    pulse_program = {}
    for scheme_name in gate_schemes.keys():
        schedules = []
        for meas_basis in bloch_basis:
            for repetition in repetitions:
                for scaling_factor_I in scalings:
                    Ival = np.real(waveform[scheme_name]) * scaling_factor_I
                    Qval = np.imag(waveform[scheme_name])
                    scaled_waveform = Ival + 1j * Qval
                    control_pulse = Play(Waveform(scaled_waveform), drive_chan)
                    schedule = pulse.Schedule(
                        name=f"pi_pulse_scal_{scaling_factor_I}_reps_"
                        + str(repetition)
                        + "_basis_"
                        + str(meas_basis)
                    )
                    for reps in np.arange(repetition):
                        schedule += control_pulse
                    if meas_basis == "x":
                        schedule += inst_sched_map.get("u2", [0], P0=0.0, P1=np.pi)
                        schedule += measure_schedule << schedule.duration
                    if meas_basis == "y":
                        schedule += inst_sched_map.get("u2", [0], P0=0.0, P1=np.pi / 2)
                        schedule += measure_schedule << schedule.duration
                    if meas_basis == "z":
                        schedule |= measure_schedule << schedule.duration
                    schedules.append(schedule)
        pulse_program[scheme_name] = assemble(
            schedules,
            backend=backend,
            meas_level=2,
            meas_return="single",
            shots=num_shots_per_point,
            schedule_los=[{drive_chan: qubit_frequency_updated}]
            * (repetitions.shape[0])
            * (scalings.shape[0])
            * 3,
        )

    # Running the jobs
    exp_results = {}
    for scheme_name in gate_schemes.keys():
        print("Control Scheme:", scheme_name)
        job = backend.run(pulse_program[scheme_name])
        print("Job id:", job.job_id())
        job_monitor(job)
        exp_results[scheme_name] = job.result(timeout=120)

    # Extract results
    I_calibration_data = {}
    for scheme_name, results in exp_results.items():
        print(scheme_name)
        idx = 0
        result_bloch = []
        for meas_basis in bloch_basis:
            full_res = []
            for idx_rep, number_of_repetitions in enumerate(repetitions):
                rep_res = []
                for scaling_factor in scalings:
                    counts = results.get_counts()[idx]
                    excited_pop = 0
                    idx += 1
                    for bits, count in counts.items():
                        excited_pop += count if bits[::-1][qubit] == "1" else 0
                    rep_res.append(excited_pop / num_shots_per_point)
                full_res.append(rep_res)
            result_bloch.append(1 - 2 * np.asarray(full_res))
        I_calibration_data[scheme_name] = np.asarray(result_bloch)

    save_var(
        "resources/superconducting-qubits-pulse-calibration/relative_calibration_data_jakarta_qubit_0",
        I_calibration_data,
    )

else:
    I_calibration_data = load_var(
        "resources/superconducting-qubits-pulse-calibration/relative_calibration_data_jakarta_qubit_0"
    )

# Infidelities for repetition experiments
target_bloch = np.array([0, -1, 0])
bloch_infidelity = {}
for scheme_name in gate_schemes.keys():
    bloch_vectors = np.array(
        [
            np.array(
                [
                    I_calibration_data[scheme_name][:, rep, idx]
                    for idx in range(scalings.shape[0])
                ]
            )
            for rep in range(repetitions.shape[0])
        ]
    )
    bloch_infidelity[scheme_name] = (
        1 - ((1 + np.dot(bloch_vectors, target_bloch)) / 2) ** 2
    )

# Plotting
fig, axs = plt.subplots(1, len(gate_schemes.keys()), figsize=(10, 5), squeeze=False)
fig.suptitle("Relative scaling calibration from pulse repetitions", y=1)

# Find the optimal value of S_rel by averaging data for N_1 and N_2 together and finding the minimum infidelity.
best_scaling_I = {}
average_infidelity = {}
for idx, scheme_name in enumerate(gate_schemes.keys()):
    ax = axs[0][idx]
    infidelity_smoothed = []
    for idx, pulse_repetition in enumerate(repetitions):
        infidelity_smoothed.append(
            movingaverage(bloch_infidelity[scheme_name][idx, :], 3)
        )
        ax.set_title(f"Control: {scheme_name}")
        ax.scatter(
            scalings,
            bloch_infidelity[scheme_name][idx, :],
            label="Infidelity data for %s pulse repetitions" % (pulse_repetition),
        )
        ax.plot(
            scalings[1:-1],
            infidelity_smoothed[-1],
            label="Smoothed infidelity for %s pulse repetitions" % (pulse_repetition),
        )
        ax.set_xlabel("Relative scaling")
        ax.set_ylabel("Infidelity")

    average_infidelity[scheme_name] = np.average(infidelity_smoothed, axis=0)
    best_scaling_I[scheme_name] = scalings[
        np.argmin(average_infidelity[scheme_name]) + 1
    ]
    ax.plot(
        scalings[1:-1],
        average_infidelity[scheme_name],
        label="Average infidelity",
        c="k",
    )
    ax.axvline(best_scaling_I[scheme_name], ls="--")
    ax.text(
        scalings[np.argmin(infidelity_smoothed) + 1],
        np.max(infidelity_smoothed),
        f"Best $S_$: {scalings[np.argmin(infidelity_smoothed) + 1]:.3f}",
        fontsize=14,
    )
    ax.set_ylim(0, 1.05)
ax.legend(loc="upper left", bbox_to_anchor=(1, 1))
plt.show()

Infidelity as a function of $S_{rel}$ for 49 and 101 repetitions of the pulse. The solid lines are a moving average of the experimental results (symbols) using the values from 3 points. The black solid line is the average between the 49 and 101 repetitions experiments. The scaling factor that you should use when calibrating your pulse is the one that minimizes the average infidelity.

The following cell now performs step 3 for $S_{amp}$.

repetitions = np.array([49, 101])
scalings = np.linspace(0.95, 1.05, 50)

if use_saved_data == False:
    pulse_program = {}
    for scheme_name in gate_schemes.keys():
        schedules = []
        for meas_basis in bloch_basis:
            for repetition in repetitions:
                for scaling_factor in scalings:
                    scaled_waveform = (
                        np.real(waveform[scheme_name]) * best_scaling_I[scheme_name]
                        + 1j * np.imag(waveform[scheme_name])
                    ) * scaling_factor
                    control_pulse = Play(Waveform(scaled_waveform), drive_chan)
                    schedule = pulse.Schedule(
                        name=f"pi_pulse_scal_{scaling_factor}_reps_"
                        + str(repetition)
                        + "_basis_"
                        + str(meas_basis)
                    )
                    for reps in np.arange(repetition):
                        schedule += control_pulse
                    if meas_basis == "x":
                        schedule += inst_sched_map.get("u2", [0], P0=0.0, P1=np.pi)
                        schedule += measure_schedule << schedule.duration
                    if meas_basis == "y":
                        schedule += inst_sched_map.get("u2", [0], P0=0.0, P1=np.pi / 2)
                        schedule += measure_schedule << schedule.duration
                    if meas_basis == "z":
                        schedule |= measure_schedule << schedule.duration
                    schedules.append(schedule)
        pulse_program[scheme_name] = assemble(
            schedules,
            backend=backend,
            meas_level=2,
            meas_return="single",
            shots=num_shots_per_point,
            schedule_los=[{drive_chan: qubit_frequency_updated}]
            * (repetitions.shape[0])
            * (scalings.shape[0])
            * 3,
        )

    # Running the jobs
    exp_results = {}
    for scheme_name in gate_schemes.keys():
        print("Control Scheme:", scheme_name)
        job = backend.run(pulse_program[scheme_name])
        print("Job id:", job.job_id())
        job_monitor(job)
        exp_results[scheme_name] = job.result(timeout=120)

    # Extract results
    amplitude_calibration_data = {}
    for scheme_name, results in exp_results.items():
        print(scheme_name)
        idx = 0
        result_bloch = []
        for meas_basis in bloch_basis:
            full_res = []
            for idx_rep, number_of_repetitions in enumerate(repetitions):
                rep_res = []
                for scaling_factor in scalings:
                    counts = results.get_counts()[idx]
                    excited_pop = 0
                    idx += 1
                    for bits, count in counts.items():
                        excited_pop += count if bits[::-1][qubit] == "1" else 0
                    rep_res.append(excited_pop / num_shots_per_point)
                full_res.append(rep_res)
            result_bloch.append(1 - 2 * np.asarray(full_res))
        amplitude_calibration_data[scheme_name] = np.asarray(result_bloch)

    save_var(
        "resources/superconducting-qubits-pulse-calibration/amplitude_calibration_data_jakarta_qubit_0",
        amplitude_calibration_data,
    )
else:
    amplitude_calibration_data = load_var(
        "resources/superconducting-qubits-pulse-calibration/amplitude_calibration_data_jakarta_qubit_0"
    )

# Infidelities for repetition experiments
bloch_infidelity = {}
for scheme_name in gate_schemes.keys():
    bloch_vectors = np.array(
        [
            np.array(
                [
                    amplitude_calibration_data[scheme_name][:, rep, idx]
                    for idx in range(scalings.shape[0])
                ]
            )
            for rep in range(repetitions.shape[0])
        ]
    )
    bloch_infidelity[scheme_name] = (
        1 - ((1 + np.dot(bloch_vectors, target_bloch)) / 2) ** 2
    )

fig, axs = plt.subplots(1, len(gate_schemes.keys()), figsize=(10, 5), squeeze=False)
fig.suptitle(f"Amplitude scaling calibration from pulse repetitions", y=1)

best_scaling_amplitude = {}
average_infidelity = {}
for idx, scheme_name in enumerate(gate_schemes.keys()):
    ax = axs[0][idx]
    infidelity_smoothed = []
    for idx, pulse_repetition in enumerate(repetitions):
        infidelity_smoothed.append(
            movingaverage(bloch_infidelity[scheme_name][idx, :], 5)
        )
        ax.set_title(f"Control: {scheme_name}")
        ax.scatter(
            scalings,
            bloch_infidelity[scheme_name][idx, :],
            label=f"Infidelity data for {pulse_repetition} pulse repetitions",
        )
        ax.plot(
            scalings[2:-2],
            infidelity_smoothed[-1],
            label=f"Smoothed infidelity for {pulse_repetition} pulse repetitions",
        )
        ax.set_xlabel("Amplitude scaling")
        ax.set_ylabel("Infidelity")

    average_infidelity[scheme_name] = np.average(infidelity_smoothed, axis=0)
    best_scaling_amplitude[scheme_name] = scalings[
        np.argmin(average_infidelity[scheme_name]) + 2
    ]
    ax.plot(
        scalings[2:-2],
        average_infidelity[scheme_name],
        label="Average infidelity",
        c="k",
    )
    ax.text(
        best_scaling_amplitude[scheme_name],
        0.35,
        f"Best $S_$: {best_scaling_amplitude[scheme_name]:.3f}",
        fontsize=14,
    )
    ax.axvline(best_scaling_amplitude[scheme_name], ls="--")
    ax.set_ylim(0, 0.5)
ax.legend(loc="upper left", bbox_to_anchor=(1, 1))
plt.show()

Infidelity as a function of the amplitude scaling factor $S_{amp}$. As in the previous step, the solid lines are a moving average of the experimental results (symbols), and the black solid line is the average between the two smoothed curves. In this figure, you can see that this particular pulse is mostly insensitive to the overall scaling factor as it has been designed to be robust against amplitude noise (check the Application Note "Designing noise-robust single-qubit gates for IBM Qiskit" for more information on the design of robust pulses).

Verifying calibration using a pulse-repetition experiment

Here we perform a simple experiment repeating the calibrated pulses (using the different methods above) and compare the performance. The selected pulse-repetition experiment amplifies coherent calibration errors.

repetitions = np.arange(1, 151, 4)
num_shots_per_point = 2048
calibration_results = {
    "Basic Rabi calibration": [1, 1],
    "Fine-tuned calibration": [
        best_scaling_I["Amplitude robust"],
        best_scaling_amplitude["Amplitude robust"],
    ],
    "Automated calibration": best_param_array[-1],
}

if use_saved_data == False:
    # Full pulse calibration
    calibrated_waveform = {}
    for scheme_name, scalings in calibration_results.items():
        calibrated_waveform[scheme_name] = (
            scalings[0] * np.real(waveform["Amplitude robust"])
            + 1j * np.imag(waveform["Amplitude robust"])
        ) * scalings[1]

    pulse_program = {}
    for scheme_name, calibrated_pulse in calibrated_waveform.items():
        schedules = []
        for meas_basis in bloch_basis:
            for repetition in repetitions:
                control_pulse = Play(Waveform(calibrated_pulse), drive_chan)
                schedule = pulse.Schedule(
                    name=f"pulse_reps_" + str(repetition) + "_basis_" + str(meas_basis)
                )
                for reps in np.arange(repetition):
                    schedule += control_pulse
                if meas_basis == "x":
                    schedule += inst_sched_map.get("u2", [0], P0=0.0, P1=np.pi)
                    schedule += measure_schedule << schedule.duration
                if meas_basis == "y":
                    schedule += inst_sched_map.get("u2", [0], P0=0.0, P1=np.pi / 2)
                    schedule += measure_schedule << schedule.duration
                if meas_basis == "z":
                    schedule |= measure_schedule << schedule.duration
                schedules.append(schedule)
        pulse_program[scheme_name] = assemble(
            schedules,
            backend=backend,
            meas_level=2,
            meas_return="single",
            shots=num_shots_per_point,
            schedule_los=[{drive_chan: qubit_frequency_updated}]
            * (repetitions.shape[0])
            * 3,
        )

    # Run the jobs
    exp_results = {}
    for scheme_name, program in pulse_program.items():
        print("Control Scheme:", scheme_name)
        job = backend.run(program)
        print("Job id:", job.job_id())
        job_monitor(job)
        exp_results[scheme_name] = job.result(timeout=120)

    # Extract results
    calibrated_pulse_performance = {}
    for scheme_name, results in exp_results.items():
        print(scheme_name)
        idx = 0
        result_bloch = []
        for meas_basis in bloch_basis:
            full_res = []
            for idx_rep, number_of_repetitions in enumerate(repetitions):
                rep_res = []
                counts = results.get_counts()[idx]
                excited_pop = 0
                for bits, count in counts.items():
                    excited_pop += count if bits[::-1][qubit] == "1" else 0
                rep_res.append(excited_pop / num_shots_per_point)
                full_res.append(rep_res)
                idx += 1
            result_bloch.append(1 - 2 * np.asarray(full_res))
        calibrated_pulse_performance[scheme_name] = np.asarray(result_bloch)

    save_var(
        "resources/superconducting-qubits-pulse-calibration/calibrated_pulse_performance_jakarta_qubit_0",
        calibrated_pulse_performance,
    )
else:
    calibrated_pulse_performance = load_var(
        "resources/superconducting-qubits-pulse-calibration/calibrated_pulse_performance_jakarta_qubit_0"
    )

# Infidelities for repetition experiments
bloch_infidelity = {}
for scheme_name, performance in calibrated_pulse_performance.items():
    bloch_vectors = performance[:, :, 0].T
    bloch_infidelity[scheme_name] = (
        1 - ((1 + np.dot(bloch_vectors, target_bloch)) / 2) ** 2
    )

# Plotting

fig = plt.figure(figsize=(10, 5))
fig.suptitle(f"Pulse performance after calibration")
for scheme_name, infidelity in bloch_infidelity.items():
    plt.plot(repetitions, infidelity, marker="o", label=scheme_name)
    plt.xlabel("Pulse repetition number")
    plt.ylabel("Infidelity")
    plt.ylim(0, 1.05)
plt.legend()
plt.show()

These data reveal that using only the basic coarse calibration we come close to an ideal implementation for a single pulse but the infidelity quickly worsens as the gate is repeated. The behavior observed is consistent with the onset of oscillations due to rotation errors. By contrast, both fine-calibration procedures remove this behavior and result in a decay which is approximately consistent with incoherent decay from the measured $T_{1}$ in the hardware.