diff --git a/runcards/rb_correction_qua.py b/runcards/rb_correction_qua.py new file mode 100644 index 000000000..222969b67 --- /dev/null +++ b/runcards/rb_correction_qua.py @@ -0,0 +1,202 @@ +from argparse import ArgumentParser +from dataclasses import dataclass, field +from pathlib import Path +from typing import Optional + +import numpy as np +import numpy.typing as npt +import plotly.graph_objects as go +from qibolab import create_platform +from qibolab.qubits import QubitId + +from qibocal.auto.execute import Executor +from qibocal.auto.operation import Data +from qibocal.cli.report import report +from qibocal.protocols.randomized_benchmarking import utils + +RBCorrectionType = np.dtype( + [ + ("biases", np.float64), + ("qubit_frequency", np.float64), + ("pulse_fidelity_uncorrected", np.float64), + ("pulse_fidelity_corrected", np.float64), + ] +) +"""Custom dtype for RBCorrection routines.""" + + +biases = np.arange(-0.02, 0.02, 0.01) +"bias points to sweep" + +# Flipping +nflips_max = 200 +"""Maximum number of flips ([RX(pi) - RX(pi)] sequences). """ +nflips_step = 10 +"""Flip step.""" + +# Ramsey signal +detuning = 3_000_000 +"""Frequency detuning [Hz].""" +delay_between_pulses_start = 16 +"""Initial delay between RX(pi/2) pulses in ns.""" +delay_between_pulses_end = 5_000 +"""Final delay between RX(pi/2) pulses in ns.""" +delay_between_pulses_step = 200 +"""Step delay between RX(pi/2) pulses in ns.""" + +# Std Rb Ondevice Parameters +num_of_sequences: int = 100 +max_circuit_depth: int = 250 +"Maximum circuit depth" +delta_clifford: int = 50 +"Play each sequence with a depth step equals to delta_clifford" +seed: Optional[int] = 1234 +"Pseudo-random number generator seed" +n_avg: int = 128 +"Number of averaging loops for each random sequence" +save_sequences: bool = True +apply_inverse: bool = True +state_discrimination: bool = True +"Flag to enable state discrimination if the readout has been calibrated (rotated blobs and threshold)" + + +@dataclass +class RBCorrectionSignalData(Data): + """Coherence acquisition outputs.""" + + data: dict[QubitId, npt.NDArray] = field(default_factory=dict) + """Raw data acquired.""" + + @property + def average(self): + if len(next(iter(self.data.values())).shape) > 1: + return utils.average_single_shots(self.__class__, self.data) + return self + + +parser = ArgumentParser() +parser.add_argument("--target", default="D2", help="Target qubit index") +parser.add_argument("--platform", type=str, default="qw11q", help="Platform name") +parser.add_argument( + "--path", type=str, default="TESTRBCorrection", help="Path for the output" +) +args = parser.parse_args() + +target = args.target +path = args.path + +data = RBCorrectionSignalData() + +platform = create_platform(args.platform) +for target in [args.target]: + + # NOTE: Center around the sweetspot [Optional] + centered_biases = biases + platform.qubits[target].sweetspot + + for i, bias in enumerate(centered_biases): + with Executor.open( + f"myexec_{i}", + path=path / Path(f"flux_{i}"), + platform=platform, + targets=[target], + update=True, + force=True, + ) as e: + + discrimination_output = e.single_shot_classification(nshots=5000) + + rb_output_uncorrected = e.rb_ondevice( + num_of_sequences=num_of_sequences, + max_circuit_depth=max_circuit_depth, + delta_clifford=delta_clifford, + seed=seed, + n_avg=n_avg, + save_sequences=save_sequences, + apply_inverse=apply_inverse, + state_discrimination=state_discrimination, + ) + + ramsey_output = e.ramsey_signal( + delay_between_pulses_start=delay_between_pulses_start, + delay_between_pulses_end=delay_between_pulses_end, + delay_between_pulses_step=delay_between_pulses_step, + detuning=detuning, + ) + + flipping_output = e.flipping_signal( + nflips_max=nflips_max, + nflips_step=nflips_step, + ) + + discrimination_output = e.single_shot_classification(nshots=5000) + + rb_output_corrected = e.rb_ondevice( + num_of_sequences=num_of_sequences, + max_circuit_depth=max_circuit_depth, + delta_clifford=delta_clifford, + seed=seed, + n_avg=n_avg, + save_sequences=save_sequences, + apply_inverse=apply_inverse, + state_discrimination=state_discrimination, + ) + + data.register_qubit( + RBCorrectionType, + (target), + dict( + biases=[bias], + qubit_frequency=[platform.qubits[target].native_gates.RX.frequency], + pulse_fidelity_uncorrected=[ + rb_output_uncorrected.results.pars[target][ + 2 + ] # Get error covs[2] + ], + pulse_fidelity_corrected=[ + rb_output_corrected.results.pars[target][2] + ], + ), + ) + + report(e.path, e.history) + + +def plot(data: RBCorrectionSignalData, target: QubitId, path): + """Plotting function for Coherence experiment.""" + + figure = go.Figure() + + figure.add_trace( + go.Scatter( + x=data[target].qubit_frequency, + y=data[target].pulse_fidelity_uncorrected, + opacity=1, + name="Pulse Fidelity Uncorrected", + showlegend=True, + legendgroup="Pulse Fidelity Uncorrected", + ) + ) + + figure.add_trace( + go.Scatter( + x=data[target].qubit_frequency, + y=data[target].pulse_fidelity_corrected, + opacity=1, + name="Pulse Fidelity Corrected", + showlegend=True, + legendgroup="Pulse Fidelity Corrected", + ) + ) + + # last part + figure.update_layout( + showlegend=True, + xaxis_title="Frequency [GHZ]", + yaxis_title="Pulse Fidelity", + ) + + if path is not None: + figure.write_html(path / Path("plot.html")) + + +plot(data, target, path=path) diff --git a/runcards/rx_calibration.py b/runcards/rx_calibration.py index ebb65e10b..f71fd20f5 100644 --- a/runcards/rx_calibration.py +++ b/runcards/rx_calibration.py @@ -1,77 +1,75 @@ from qibocal.auto.execute import Executor from qibocal.cli.report import report -executor = Executor.create(name="myexec", platform="dummy") - -from myexec import close, drag_tuning, init, rabi_amplitude, ramsey - target = 0 -platform = executor.platform -platform.settings.nshots = 2048 -init("test_rx_calibration", force=True, targets=[target]) +with Executor.open( + "myexec", + path="test_rx_calibration", + platform="dummy", + targets=[target], + update=True, + force=True, +) as e: + e.platform.settings.nshots = 2000 -rabi_output = rabi_amplitude( - min_amp_factor=0.5, - max_amp_factor=1.5, - step_amp_factor=0.01, - pulse_length=platform.qubits[target].native_gates.RX.duration, -) -# update only if chi2 is satisfied -if rabi_output.results.chi2[target][0] > 2: - raise RuntimeError( - f"Rabi fit has chi2 {rabi_output.results.chi2[target][0]} greater than 2. Stopping." + rabi_output = e.rabi_amplitude( + min_amp_factor=0.5, + max_amp_factor=1.5, + step_amp_factor=0.01, + pulse_length=e.platform.qubits[target].native_gates.RX.duration, ) -rabi_output.update_platform(platform) + # update only if chi2 is satisfied + if rabi_output.results.chi2[target][0] > 2: + raise RuntimeError( + f"Rabi fit has chi2 {rabi_output.results.chi2[target][0]} greater than 2. Stopping." + ) -ramsey_output = ramsey( - delay_between_pulses_start=10, - delay_between_pulses_end=5000, - delay_between_pulses_step=100, - detuning=1_000_000, -) -if ramsey_output.results.chi2[target][0] > 2: - raise RuntimeError( - f"Ramsey fit has chi2 {ramsey_output.results.chi2[target][0]} greater than 2. Stopping." + ramsey_output = e.ramsey( + delay_between_pulses_start=10, + delay_between_pulses_end=5000, + delay_between_pulses_step=100, + detuning=1_000_000, + update=False, ) -if ramsey_output.results.delta_phys[target][0] < 1e4: - print( - f"Ramsey frequency not updated, correction too small { ramsey_output.results.delta_phys[target][0]}" - ) -else: - ramsey_output.update_platform(platform) + if ramsey_output.results.chi2[target][0] > 2: + raise RuntimeError( + f"Ramsey fit has chi2 {ramsey_output.results.chi2[target][0]} greater than 2. Stopping." + ) + if ramsey_output.results.delta_phys[target][0] < 1e4: + print( + f"Ramsey frequency not updated, correction too small { ramsey_output.results.delta_phys[target][0]}" + ) + else: + ramsey_output.update_platform(e.platform) -rabi_output_2 = rabi_amplitude( - min_amp_factor=0.5, - max_amp_factor=1.5, - step_amp_factor=0.01, - pulse_length=platform.qubits[target].native_gates.RX.duration, -) -# update only if chi2 is satisfied -if rabi_output_2.results.chi2[target][0] > 2: - raise RuntimeError( - f"Rabi fit has chi2 {rabi_output_2.results.chi2[target][0]} greater than 2. Stopping." + rabi_output_2 = e.rabi_amplitude( + min_amp_factor=0.5, + max_amp_factor=1.5, + step_amp_factor=0.01, + pulse_length=e.platform.qubits[target].native_gates.RX.duration, ) -rabi_output_2.update_platform(platform) + # update only if chi2 is satisfied + if rabi_output_2.results.chi2[target][0] > 2: + raise RuntimeError( + f"Rabi fit has chi2 {rabi_output_2.results.chi2[target][0]} greater than 2. Stopping." + ) -drag_output = drag_tuning(beta_start=-4, beta_end=4, beta_step=0.5) -if drag_output.results.chi2[target][0] > 2: - raise RuntimeError( - f"Drag fit has chi2 {drag_output.results.chi2[target][0]} greater than 2. Stopping." - ) -drag_output.update_platform(platform) + drag_output = e.drag_tuning(beta_start=-4, beta_end=4, beta_step=0.5) + if drag_output.results.chi2[target][0] > 2: + raise RuntimeError( + f"Drag fit has chi2 {drag_output.results.chi2[target][0]} greater than 2. Stopping." + ) -rabi_output_3 = rabi_amplitude( - min_amp_factor=0.5, - max_amp_factor=1.5, - step_amp_factor=0.01, - pulse_length=platform.qubits[target].native_gates.RX.duration, -) -# update only if chi2 is satisfied -if rabi_output_3.results.chi2[target][0] > 2: - raise RuntimeError( - f"Rabi fit has chi2 {rabi_output_3.results.chi2[target][0]} greater than 2. Stopping." + rabi_output_3 = e.rabi_amplitude( + min_amp_factor=0.5, + max_amp_factor=1.5, + step_amp_factor=0.01, + pulse_length=e.platform.qubits[target].native_gates.RX.duration, ) -rabi_output_3.update_platform(platform) + # update only if chi2 is satisfied + if rabi_output_3.results.chi2[target][0] > 2: + raise RuntimeError( + f"Rabi fit has chi2 {rabi_output_3.results.chi2[target][0]} greater than 2. Stopping." + ) -close() -report(executor.path, executor.history) +report(e.path, e.history) diff --git a/runcards/single_shot.py b/runcards/single_shot.py index 0025995fd..714b909d2 100644 --- a/runcards/single_shot.py +++ b/runcards/single_shot.py @@ -1,13 +1,30 @@ -from qibocal.auto.execute import Executor -from qibocal.cli.report import report +"""Minimal Qibocal script example. + +In this example, the default Qibocal executor is used. Additional +configurations can still be passed through the `init()` call, but there +is no explicit API to access the execution details. + +If more fine grained control is needed, refer to the `rx_calibration.py` example. -executor = Executor.create(name="myexec", platform="dummy") +.. note:: -from myexec import close, init, single_shot_classification + though simple, this example is not limited to single protocol execution, but + multiple protocols can be added as well, essentially in the same fashion of a plain + runcard - still with the advantage of handling execution and results + programmatically +""" + +from pathlib import Path + +from qibocal.cli.report import report +from qibocal.routines import close, init, single_shot_classification -init("test_x", force=True) +path = Path("test_x") +init(platform="dummy", path=path, force=True) -completed = single_shot_classification(nshots=1000) +# Here +ssc = single_shot_classification(nshots=1000) +print("\nfidelities:\n", ssc.results.fidelity, "\n") close() -report(executor.path, executor.history) +report(path) diff --git a/src/qibocal/auto/execute.py b/src/qibocal/auto/execute.py index a33b34ced..ec4db1659 100644 --- a/src/qibocal/auto/execute.py +++ b/src/qibocal/auto/execute.py @@ -4,6 +4,7 @@ import importlib.util import os import sys +from contextlib import contextmanager from copy import deepcopy from dataclasses import dataclass, fields from pathlib import Path @@ -205,6 +206,8 @@ def wrapper( parameters: Optional[dict] = None, id: str = operation, mode: ExecutionMode = AUTOCALIBRATION, + update: bool = True, + targets: Optional[Targets] = None, **kwargs, ): positional = dict( @@ -215,6 +218,8 @@ def wrapper( source={ "id": id, "operation": operation, + "targets": targets, + "update": update, "parameters": params | positional | kwargs, } ) @@ -248,6 +253,7 @@ def init( path: os.PathLike, force: bool = False, platform: Union[Platform, str, None] = None, + update: Optional[bool] = None, targets: Optional[Targets] = None, ): """Initialize execution.""" @@ -258,6 +264,8 @@ def init( platform = self.platform = backend.platform assert isinstance(platform, Platform) + if update is not None: + self.update = update if targets is not None: self.targets = targets @@ -295,3 +303,42 @@ def close(self, path: Optional[os.PathLike] = None): # attempt unloading self.__del__() + + @classmethod + @contextmanager + def open( + cls, + name: str, + path: os.PathLike, + force: bool = False, + platform: Union[Platform, str, None] = None, + update: Optional[bool] = None, + targets: Optional[Targets] = None, + ): + """Enter the execution context.""" + ex = cls.create(name, platform) + ex.init(path, force, platform, update, targets) + try: + yield ex + finally: + ex.close() + + def __enter__(self): + """Reenter the execution context. + + This method its here to reuse an already existing (and + initialized) executor, in a new context. + + It should not be used with new executors. In which case, cf. :meth:`__open__`. + """ + # connect and initialize platform + self.platform.connect() + return self + + def __exit__(self, exc_type, exc_value, traceback): + """Exit execution context. + + This pairs with :meth:`__enter__`. + """ + self.close() + return False diff --git a/src/qibocal/auto/operation.py b/src/qibocal/auto/operation.py index 8b1f602eb..1d7527ae8 100644 --- a/src/qibocal/auto/operation.py +++ b/src/qibocal/auto/operation.py @@ -181,7 +181,7 @@ def load_data(path: Path, filename: str): """Load data stored in a npz file.""" file = path / f"{filename}.npz" if file.is_file(): - raw_data_dict = dict(np.load(file)) + raw_data_dict = dict(np.load(file, allow_pickle=True)) data_dict = {} for data_key, array in raw_data_dict.items(): diff --git a/src/qibocal/cli/report.py b/src/qibocal/cli/report.py index cb6d89e42..c4576b100 100644 --- a/src/qibocal/cli/report.py +++ b/src/qibocal/cli/report.py @@ -44,14 +44,17 @@ def plotter( completed node on specific target. """ figures, fitting_report = generate_figures_and_report(node, target) - buffer = io.StringIO() - html_list = [] - for figure in figures: - figure.write_html(buffer, include_plotlyjs=False, full_html=False) - buffer.seek(0) - html_list.append(buffer.read()) - buffer.close() - all_html = "".join(html_list) + if isinstance(figures[0], str): + all_html = "".join(figures) + else: + buffer = io.StringIO() + html_list = [] + for figure in figures: + figure.write_html(buffer, include_plotlyjs=False, full_html=False) + buffer.seek(0) + html_list.append(buffer.read()) + buffer.close() + all_html = "".join(html_list) return all_html, fitting_report diff --git a/src/qibocal/protocols/__init__.py b/src/qibocal/protocols/__init__.py index 19345fa05..eef55a1c8 100644 --- a/src/qibocal/protocols/__init__.py +++ b/src/qibocal/protocols/__init__.py @@ -27,6 +27,7 @@ from .flux_dependence.qubit_flux_tracking import qubit_flux_tracking from .flux_dependence.resonator_crosstalk import resonator_crosstalk from .flux_dependence.resonator_flux_dependence import resonator_flux +from .qua import rb_ondevice from .qubit_power_spectroscopy import qubit_power_spectroscopy from .qubit_spectroscopy import qubit_spectroscopy from .qubit_spectroscopy_ef import qubit_spectroscopy_ef @@ -144,4 +145,5 @@ "rabi_length_frequency", "rabi_length_frequency_signal", "standard_rb_2q", + "rb_ondevice", ] diff --git a/src/qibocal/protocols/qua/__init__.py b/src/qibocal/protocols/qua/__init__.py new file mode 100644 index 000000000..deeb29b40 --- /dev/null +++ b/src/qibocal/protocols/qua/__init__.py @@ -0,0 +1 @@ +from .rb_ondevice import rb_ondevice diff --git a/src/qibocal/protocols/qua/configuration.py b/src/qibocal/protocols/qua/configuration.py new file mode 100644 index 000000000..d457e3958 --- /dev/null +++ b/src/qibocal/protocols/qua/configuration.py @@ -0,0 +1,160 @@ +from dataclasses import asdict + +from qibolab.instruments.qm import QMController +from qibolab.instruments.qm.config import QMConfig + +NATIVE_OPS = { + "x180": lambda q: (f"plus_i_{q}", f"plus_q_{q}"), + "y180": lambda q: (f"minus_q_{q}", f"plus_i_{q}"), + "x90": lambda q: (f"plus_i_half_{q}", f"plus_q_half_{q}"), + "y90": lambda q: (f"minus_q_half_{q}", f"plus_i_half_{q}"), + "-x90": lambda q: (f"minus_i_half_{q}", f"minus_q_half_{q}"), + "-y90": lambda q: (f"plus_q_half_{q}", f"minus_i_half_{q}"), +} + + +def native_operations(qubit): + return {op: f"{op}_{qubit}" for op in NATIVE_OPS.keys()} + + +def drive_waveform_components(qubit, mode, samples): + return { + f"plus_{mode}_{qubit}": { + "type": "arbitrary", + "samples": samples, + }, + f"minus_{mode}_{qubit}": { + "type": "arbitrary", + "samples": -samples, + }, + f"plus_{mode}_half_{qubit}": { + "type": "arbitrary", + "samples": samples / 2, + }, + f"minus_{mode}_half_{qubit}": { + "type": "arbitrary", + "samples": -samples / 2, + }, + } + + +def drive_waveforms(platform, qubit): + pulse = platform.qubits[qubit].native_gates.RX.pulse(start=0) + envelope_i, envelope_q = pulse.envelope_waveforms(sampling_rate=1) + return drive_waveform_components( + qubit, "i", envelope_i.data + ) | drive_waveform_components(qubit, "q", envelope_q.data) + + +def waveforms(platform, qubits): + _waveforms = { + "zero": { + "type": "constant", + "sample": 0.0, + }, + } + _waveforms.update( + { + f"mz_{q}": { + "type": "constant", + "sample": platform.qubits[q].native_gates.MZ.amplitude, + } + for q in qubits + } + ) + for q in qubits: + _waveforms.update(drive_waveforms(platform, q)) + return _waveforms + + +def drive_pulses(platform, qubit): + _pulses = {} + for op, wf in NATIVE_OPS.items(): + i, q = wf(qubit) + _pulses[f"{op}_{qubit}"] = { + "operation": "control", + "length": platform.qubits[qubit].native_gates.RX.duration, + "waveforms": { + "I": i, + "Q": q, + }, + "digital_marker": "ON", + } + return _pulses + + +def pulses(platform, qubits): + _pulses = { + f"mz_{q}": { + "operation": "measurement", + "length": platform.qubits[q].native_gates.MZ.duration, + "waveforms": { + "I": f"mz_{q}", + "Q": "zero", + }, + "integration_weights": { + "cos": f"cosine_weights{q}", + "sin": f"sine_weights{q}", + "minus_sin": f"minus_sine_weights{q}", + }, + "digital_marker": "ON", + } + for q in qubits + } + for q in qubits: + _pulses.update(drive_pulses(platform, q)) + return _pulses + + +def integration_weights(platform, qubits): + _integration_weights = {} + for q in qubits: + _duration = platform.qubits[q].native_gates.MZ.duration + _integration_weights.update( + { + f"cosine_weights{q}": { + "cosine": [(1.0, _duration)], + "sine": [(-0.0, _duration)], + }, + f"sine_weights{q}": { + "cosine": [(0.0, _duration)], + "sine": [(1.0, _duration)], + }, + f"minus_sine_weights{q}": { + "cosine": [(-0.0, _duration)], + "sine": [(-1.0, _duration)], + }, + } + ) + return _integration_weights + + +def register_element(config, qubit, time_of_flight, smearing): + config.register_port(qubit.readout.port) + config.register_port(qubit.feedback.port) + mz_frequency = qubit.native_gates.MZ.frequency - qubit.readout.lo_frequency + config.register_readout_element(qubit, mz_frequency, time_of_flight, smearing) + config.register_port(qubit.drive.port) + rx_frequency = qubit.native_gates.RX.frequency - qubit.drive.lo_frequency + config.register_drive_element(qubit, rx_frequency) + if qubit.flux is not None: + config.register_port(qubit.flux.port) + config.register_flux_element(qubit) + + +def generate_config(platform, qubits): + con = [ + instr + for instr in platform.instruments.values() + if isinstance(instr, QMController) + ][0] + config = QMConfig() + for q in qubits: + qubit = platform.qubits[q] + register_element(config, qubit, con.time_of_flight, con.smearing) + config.elements[f"readout{q}"]["operations"]["measure"] = f"mz_{q}" + config.elements[f"drive{q}"]["operations"] = native_operations(q) + config.pulses = pulses(platform, qubits) + config.waveforms = waveforms(platform, qubits) + config.integration_weights = integration_weights(platform, qubits) + return asdict(config) diff --git a/src/qibocal/protocols/qua/qua_clifford_group.npz b/src/qibocal/protocols/qua/qua_clifford_group.npz new file mode 100644 index 000000000..4cc9f4d1c Binary files /dev/null and b/src/qibocal/protocols/qua/qua_clifford_group.npz differ diff --git a/src/qibocal/protocols/qua/rb_ondevice.py b/src/qibocal/protocols/qua/rb_ondevice.py new file mode 100644 index 000000000..5d11473d2 --- /dev/null +++ b/src/qibocal/protocols/qua/rb_ondevice.py @@ -0,0 +1,544 @@ +import time +from dataclasses import dataclass, field +from typing import Optional + +import matplotlib.pyplot as plt +import mpld3 +import numpy as np +import numpy.typing as npt +from qibolab.platform import Platform +from qibolab.qubits import QubitId +from qm import generate_qua_script +from qm.qua import * # nopycln: import +from qualang_tools.bakery.randomized_benchmark_c1 import c1_table +from qualang_tools.results import fetching_tool +from scipy.optimize import curve_fit + +from qibocal.auto.operation import Data, Parameters, Results, Routine +from qibocal.protocols.utils import table_dict, table_html + +from .configuration import generate_config +from .utils import RBType, generate_depths, power_law, process_data + +# parser.add_argument("--simulation-duration", type=int, default=None) +# parser.add_argument("--relaxation-time", type=int, default=None) + + +@dataclass +class RbOnDeviceParameters(Parameters): + num_of_sequences: int + max_circuit_depth: int + "Maximum circuit depth" + delta_clifford: int + "Play each sequence with a depth step equals to delta_clifford" + seed: Optional[int] = None + "Pseudo-random number generator seed" + n_avg: int = 1 + "Number of averaging loops for each random sequence" + save_sequences: bool = True + apply_inverse: bool = False + state_discrimination: bool = True + "Flag to enable state discrimination if the readout has been calibrated (rotated blobs and threshold)" + + def __post_init__(self): + if self.seed is None: + self.seed = np.random.randint(0, int(1e6)) + + +def generate_sequence(max_circuit_depth, seed): + sequence = declare(int, size=max_circuit_depth + 1) + i = declare(int) + rand = Random(seed=seed) + + with for_(i, 0, i < max_circuit_depth, i + 1): + assign(sequence[i], rand.rand_int(24)) + + return sequence + + +def generate_sequence_with_inverse(max_circuit_depth, seed, c1_table, inv_gates): + cayley = declare(int, value=c1_table.flatten().tolist()) + inv_list = declare(int, value=inv_gates) + current_state = declare(int) + step = declare(int) + sequence = declare(int, size=max_circuit_depth + 1) + inv_gate = declare(int, size=max_circuit_depth + 1) + i = declare(int) + rand = Random(seed=seed) + + assign(current_state, 0) + with for_(i, 0, i < max_circuit_depth, i + 1): + assign(step, rand.rand_int(24)) + assign(current_state, cayley[current_state * 24 + step]) + assign(sequence[i], step) + assign(inv_gate[i], inv_list[current_state]) + + return sequence, inv_gate + + +def play_sequence(qubit, sequence_list, apply_inverse, depth, rx_duration): + i = declare(int) + condition = (i <= depth) if apply_inverse else (i < depth) + with for_(i, 0, condition, i + 1): + with switch_(sequence_list[i], unsafe=True): + with case_(0): + wait(rx_duration // 4, qubit) + with case_(1): + play("x180", qubit) + with case_(2): + play("y180", qubit) + with case_(3): + play("y180", qubit) + play("x180", qubit) + with case_(4): + play("x90", qubit) + play("y90", qubit) + with case_(5): + play("x90", qubit) + play("-y90", qubit) + with case_(6): + play("-x90", qubit) + play("y90", qubit) + with case_(7): + play("-x90", qubit) + play("-y90", qubit) + with case_(8): + play("y90", qubit) + play("x90", qubit) + with case_(9): + play("y90", qubit) + play("-x90", qubit) + with case_(10): + play("-y90", qubit) + play("x90", qubit) + with case_(11): + play("-y90", qubit) + play("-x90", qubit) + with case_(12): + play("x90", qubit) + with case_(13): + play("-x90", qubit) + with case_(14): + play("y90", qubit) + with case_(15): + play("-y90", qubit) + with case_(16): + play("-x90", qubit) + play("y90", qubit) + play("x90", qubit) + with case_(17): + play("-x90", qubit) + play("-y90", qubit) + play("x90", qubit) + with case_(18): + play("x180", qubit) + play("y90", qubit) + with case_(19): + play("x180", qubit) + play("-y90", qubit) + with case_(20): + play("y180", qubit) + play("x90", qubit) + with case_(21): + play("y180", qubit) + play("-x90", qubit) + with case_(22): + play("x90", qubit) + play("y90", qubit) + play("x90", qubit) + with case_(23): + play("-x90", qubit) + play("y90", qubit) + play("-x90", qubit) + + +RbOnDeviceType = np.dtype( + [ + ("state", np.int32), + ("sequences", np.int32), + ] +) + + +@dataclass +class RbOnDeviceData(Data): + rb_type: str + relaxation_time: int + depths: list[int] + data: dict[QubitId, dict[str, npt.NDArray[np.int32]]] + + def _get_data(self, key: str): + qubit = self.qubits[0] + try: + arrays = self.data[qubit].item(0) + return arrays[key] + except AttributeError: + return self.data[qubit][key] + + @property + def state(self): + return self._get_data("state") + + @property + def sequences(self): + return self._get_data("sequences") + + +def _acquisition( + params: RbOnDeviceParameters, platform: Platform, targets: list[QubitId] +) -> RbOnDeviceData: + assert len(targets) == 1 + target = targets[0] + qubit = f"drive{target}" + resonator = f"readout{target}" if target is not None else f"readout{target}" + save_sequences = params.save_sequences + apply_inverse = params.apply_inverse + relaxation_time = params.relaxation_time + + num_of_sequences = params.num_of_sequences + n_avg = params.n_avg + max_circuit_depth = params.max_circuit_depth + delta_clifford = params.delta_clifford + assert ( + max_circuit_depth / delta_clifford + ).is_integer(), "max_circuit_depth / delta_clifford must be an integer." + seed = params.seed + state_discrimination = params.state_discrimination + # List of recovery gates from the lookup table + inv_gates = [int(np.where(c1_table[i, :] == 0)[0][0]) for i in range(24)] + + ################### + # The QUA program # + ################### + with program() as rb: + depth = declare(int) # QUA variable for the varying depth + depth_target = declare( + int + ) # QUA variable for the current depth (changes in steps of delta_clifford) + # QUA variable to store the last Clifford gate of the current sequence which is replaced by the recovery gate + saved_gate = declare(int) + m = declare(int) # QUA variable for the loop over random sequences + n = declare(int) # QUA variable for the averaging loop + I = declare(fixed) # QUA variable for the 'I' quadrature + Q = declare(fixed) # QUA variable for the 'Q' quadrature + state = declare(bool) # QUA variable for state discrimination + # The relevant streams + m_st = declare_stream() + if state_discrimination: + state_st = declare_stream() + else: + I_st = declare_stream() + Q_st = declare_stream() + # save random sequences to return to host + if save_sequences: + sequence_st = declare_stream() + + with for_( + m, 0, m < num_of_sequences, m + 1 + ): # QUA for_ loop over the random sequences + if apply_inverse: + sequence_list, inv_gate_list = generate_sequence_with_inverse( + max_circuit_depth, seed, c1_table, inv_gates + ) # Generate the random sequence of length max_circuit_depth + else: + sequence_list = generate_sequence(max_circuit_depth, seed) + + if delta_clifford == 1: + assign(depth_target, 1) + else: + assign(depth_target, 0) # Initialize the current depth to 0 + + if save_sequences: + save(sequence_list[0], sequence_st) # save sequence indices in stream + with for_( + depth, 1, depth <= max_circuit_depth, depth + 1 + ): # Loop over the depths + if save_sequences: + save(sequence_list[depth], sequence_st) + # Replacing the last gate in the sequence with the sequence's inverse gate + # The original gate is saved in 'saved_gate' and is being restored at the end + assign(saved_gate, sequence_list[depth]) + if apply_inverse: + assign(sequence_list[depth], inv_gate_list[depth - 1]) + # Only played the depth corresponding to target_depth + with if_((depth == 1) | (depth == depth_target)): + with for_(n, 0, n < n_avg, n + 1): # Averaging loop + # Can be replaced by active reset + if relaxation_time > 0: + wait(relaxation_time // 4, resonator) + # Align the two elements to play the sequence after qubit initialization + align(resonator, qubit) + # The strict_timing ensures that the sequence will be played without gaps + with strict_timing_(): + # Play the random sequence of desired depth + rx_duration = platform.qubits[ + target + ].native_gates.RX.duration + play_sequence( + qubit, sequence_list, apply_inverse, depth, rx_duration + ) + # Align the two elements to measure after playing the circuit. + align(qubit, resonator) + # wait(3 * (depth + 1) * (DRIVE_DURATION // 4), resonator) + # Make sure you updated the ge_threshold and angle if you want to use state discrimination + # state, I, Q = readout_macro(threshold=ge_threshold, state=state, I=I, Q=Q) + measure( + "measure", + resonator, + None, + dual_demod.full("cos", "out1", "sin", "out2", I), + dual_demod.full("minus_sin", "out1", "cos", "out2", Q), + ) + # Save the results to their respective streams + if state_discrimination: + threshold = platform.qubits[target].threshold + iq_angle = platform.qubits[target].iq_angle + cos = np.cos(iq_angle) + sin = np.sin(iq_angle) + assign(state, I * cos - Q * sin > threshold) + save(state, state_st) + else: + save(I, I_st) + save(Q, Q_st) + # Go to the next depth + assign(depth_target, depth_target + delta_clifford) + # Reset the last gate of the sequence back to the original Clifford gate + # (that was replaced by the recovery gate at the beginning) + assign(sequence_list[depth], saved_gate) + # Save the counter for the progress bar + save(m, m_st) + + with stream_processing(): + ndepth = max_circuit_depth / delta_clifford + int(delta_clifford > 1) + m_st.save("iteration") + if save_sequences: + sequence_st.buffer(max_circuit_depth + 1).buffer(num_of_sequences).save( + "sequences" + ) + + if state_discrimination: + if n_avg > 1: + # saves a 2D array of depth and random pulse sequences in order to get error bars along the random sequences + state_st.boolean_to_int().buffer(n_avg).map( + FUNCTIONS.average() + ).buffer(ndepth).buffer(num_of_sequences).save("state") + # returns a 1D array of averaged random pulse sequences vs depth of circuit for live plotting + state_st.boolean_to_int().buffer(n_avg).map( + FUNCTIONS.average() + ).buffer(ndepth).average().save("state_avg") + else: + # saves a 2D array of depth and random pulse sequences in order to get error bars along the random sequences + state_st.boolean_to_int().buffer(ndepth).buffer( + num_of_sequences + ).save("state") + # returns a 1D array of averaged random pulse sequences vs depth of circuit for live plotting + state_st.boolean_to_int().buffer(ndepth).average().save("state_avg") + else: + I_st.buffer(n_avg).map(FUNCTIONS.average()).buffer(ndepth).buffer( + num_of_sequences + ).save("I") + Q_st.buffer(n_avg).map(FUNCTIONS.average()).buffer(ndepth).buffer( + num_of_sequences + ).save("Q") + I_st.buffer(n_avg).map(FUNCTIONS.average()).buffer( + ndepth + ).average().save("I_avg") + Q_st.buffer(n_avg).map(FUNCTIONS.average()).buffer( + ndepth + ).average().save("Q_avg") + + # Print total relaxation time (estimate of execution time) + total_relaxation = ( + relaxation_time + * num_of_sequences + * n_avg + * (max_circuit_depth // delta_clifford) + * 1e-9 + ) + print("\nTotal relaxation time: %.2f sec\n" % total_relaxation) + + ##################################### + # Open Communication with the QOP # + ##################################### + controller = platform._controller + qmm = controller.manager + + ########################### + # Run or Simulate Program # + ########################### + # Open the quantum machine + config = generate_config(platform, list(platform.qubits.keys())) + qm = qmm.open_qm(config) + # Send the QUA program to the OPX, which compiles and executes it + job = qm.execute(rb) + + with open("qua_script.py", "w") as file: + file.write(generate_qua_script(rb, config)) + + # Get results from QUA program + 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" + ) + + start_time = time.time() + while results.is_processing(): + # data analysis + if state_discrimination: + state_avg, iteration = results.fetch_all() + value_avg = state_avg + else: + I, Q, iteration = results.fetch_all() + value_avg = I + final_time = time.time() + + # At the end of the program, fetch the non-averaged results to get the error-bars + rb_type = RBType.infer(apply_inverse, relaxation_time).value + if state_discrimination: + if save_sequences: + results = fetching_tool(job, data_list=["state", "sequences"]) + state, sequences = results.fetch_all() + data = {target: {"state": state, "sequences": sequences}} + else: + results = fetching_tool(job, data_list=["state"]) + state = results.fetch_all()[0] + data = {target: {"state": state}} + else: + raise NotImplementedError + # if save_sequences: + # results = fetching_tool(job, data_list=["I", "Q", "sequences"]) + # I, Q, sequences = results.fetch_all() + # data.voltage_i[target] = I + # data.voltage_q[target] = Q + # data.sequences[target] = sequences + # else: + # results = fetching_tool(job, data_list=["I", "Q"]) + # I, Q = results.fetch_all() + # data.voltage_i[target] = I + # data.voltage_q[target] = Q + return RbOnDeviceData( + rb_type=rb_type, + relaxation_time=relaxation_time, + depths=[int(x) for x in generate_depths(max_circuit_depth, delta_clifford)], + data=data, + ) + + +@dataclass +class RbOnDeviceResults(Results): + pars: dict[QubitId, list[float]] = field(default_factory=dict) + cov: dict[QubitId, list[float]] = field(default_factory=dict) + + +def _fit(data: RbOnDeviceData) -> RbOnDeviceResults: + qubit = data.qubits[0] + rb_type = RBType(data.rb_type) + depths = data.depths + state = data.state + sequences = data.sequences + + value_avg = np.mean(state, axis=0) + error_avg = np.std(state, axis=0) + + results = RbOnDeviceResults() + ydata, ysigma = process_data(rb_type, state, depths, sequences) + + pars, cov = curve_fit( + f=power_law, + xdata=depths, + ydata=ydata, + sigma=ysigma, + p0=[0.5, 0.0, 0.9], + bounds=(-np.inf, np.inf), + maxfev=2000, + ) + results.pars[qubit] = list(pars) + results.cov[qubit] = list(cov.flatten()) + + return results + + +def _plot(data: RbOnDeviceData, target: QubitId, fit: RbOnDeviceResults): + pars = fit.pars.get(target) + rb_type = RBType(data.rb_type) + relaxation_time = data.relaxation_time + depths = data.depths + state = data.state + sequences = data.sequences + + fitting_report = table_html( + table_dict( + target, + [ + "RB type", + "Number of sequences", + "Relaxation time (us)", + ], + [ + (rb_type.value.capitalize(),), + (len(state),), + (data.relaxation_time // 1000,), + ], + ) + ) + if pars is not None: + stdevs = np.sqrt(np.diag(np.reshape(fit.cov[target], (3, 3)))) + 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 + fitting_report = "\n".join( + [ + fitting_report, + table_html( + table_dict( + target, + [ + "A", + "B", + "p", + "Error rate (1-p)", + "Clifford set infidelity", + "Gate infidelity", + ], + [ + (np.round(pars[0], 3), np.round(stdevs[0], 3)), + (np.round(pars[1], 3), np.round(stdevs[1], 3)), + (np.round(pars[2], 3), np.round(stdevs[2], 3)), + (one_minus_p, stdevs[2]), + (r_c, r_c_std), + (r_g, r_g_std), + ], + display_error=True, + ) + ), + ] + ) + + ydata, ysigma = process_data(rb_type, state, depths, sequences) + + fig = plt.figure(figsize=(16, 6)) + title = f"{rb_type.value.capitalize()} RB" + plt.errorbar( + depths, ydata, ysigma, marker="o", linestyle="-", markersize=4, label="data" + ) + if pars is not None: + max_circuit_depth = depths[-1] + x = np.linspace(0, max_circuit_depth + 0.1, 1000) + plt.plot(x, power_law(x, *pars), linestyle="--", label="fit") + plt.xlabel("Depth") + plt.ylabel("Survival probability") + plt.legend() + + figures = [mpld3.fig_to_html(fig)] + return figures, fitting_report + + +def _update(results: RbOnDeviceResults, platform: Platform, target: QubitId): + pass + + +rb_ondevice = Routine(_acquisition, _fit, _plot, _update) diff --git a/src/qibocal/protocols/qua/utils.py b/src/qibocal/protocols/qua/utils.py new file mode 100644 index 000000000..8d4e76c3e --- /dev/null +++ b/src/qibocal/protocols/qua/utils.py @@ -0,0 +1,110 @@ +from enum import Enum +from pathlib import Path + +import numpy as np + +CLIFFORD = np.load(Path(__file__).parent / "qua_clifford_group.npz") + + +RESTLESS_RELAXATION_CUTOFF = int(20e3) + + +class RBType(Enum): + STANDARD = "standard" + FILTERED = "filtered" + RESTLESS = "restless" + + @classmethod + def infer(cls, apply_inverse, relaxation_time): + if apply_inverse: + return cls.STANDARD + if relaxation_time < RESTLESS_RELAXATION_CUTOFF: + assert not apply_inverse + return cls.RESTLESS + return cls.FILTERED + + +def process_data(rb_type, state, depths=None, sequences=None): + if rb_type is RBType.STANDARD: + return 1 - np.mean(state, axis=0), np.std(state, axis=0) / np.sqrt( + state.shape[0] + ) + + is_restless = rb_type is RBType.RESTLESS + term = filter_term(depths, state, sequences, is_restless=is_restless) + ff = filter_function(term) + return np.mean(ff, axis=1), np.std(ff, axis=1) / np.sqrt(ff.shape[1]) + + +def power_law(power, a, b, p): + """Function to fit to survival probability vs circuit depths.""" + return a * (p**power) + b + + +def clifford_mul(sequences, interleave=None): + """Quick multiplication over single-qubit Clifford group (using integers). + + Requires the `CLIFFORD` file (loaded on top) which contains + the group multiplication matrix and the mapping from integers to + unitary matrices. + The indexing is the one defined by the QUA script used for data acquisition. + """ + group = CLIFFORD["algebra"] + result = np.zeros_like(sequences[:, 0]) + for i in range(sequences.shape[1]): + if interleave is not None: + matrix = group[sequences[:, i], interleave] + else: + matrix = sequences[:, i] + result = group[matrix, result] + return result + + +def generate_depths(max_circuit_depth, delta_clifford): + """Generate vector of depths compatible with the QUA acquisition script.""" + if delta_clifford == 1: + return np.arange(1, max_circuit_depth + 0.1, delta_clifford).astype(int) + depths = np.arange(0, max_circuit_depth + 0.1, delta_clifford).astype(int) + depths[0] = 1 + return depths + + +def filter_term(depths, state, sequences, is_restless=False, interleave=None): + """Calculate Clifford sandwich term that appears in the filter function. + + Args: + depths (np.ndarray): Vector of depths. size: (ndepths,) + state (np.ndarray): Matrix with acquired shots. size (num_of_sequences, ndepths) + sequences (np.ndarray): Matrix with Clifford indices used. size (num_of_sequences, max_circuit_depth) + interleave (int): Optional integer from 0 to 23, corresponding to the Clifford matrix to interleave. + is_restless (bool): If `True` the restless filter function is used. + """ + state = state.astype(int) + seqids = np.arange(len(sequences)) + terms = [] + for i, depth in enumerate(depths): + clifford_indices = clifford_mul(sequences[:, :depth], interleave=interleave) + # `clifford_indices`: (num_of_sequences,) + clifford_matrices = CLIFFORD["matrices"][clifford_indices] + # `clifford_matrices`: (num_of_sequences, 2, 2) + if is_restless: + if i > 0: + state_before = state[:, i - 1] + else: + state_before = np.concatenate(([0], state[:-1, i - 1])) + else: + state_before = 0 + terms.append(clifford_matrices[seqids, state[:, i], state_before]) + return terms + + +def filter_function(x): + """Calculate filter function using Eq. (7) from notes. + + Args: + x: Term calculated by the ``filter_term`` function above. + + Returns: + Filter function output for each data point of shape (ndepths, num_of_sequences). + """ + return 3 * (np.abs(x) ** 2 - 0.5) diff --git a/tests/test_executor.py b/tests/test_executor.py index 84417978f..bbecfca1a 100644 --- a/tests/test_executor.py +++ b/tests/test_executor.py @@ -151,3 +151,13 @@ def test_close(tmp_path: Path, executor: Executor): assert executor.meta is not None assert executor.meta.start is not None assert executor.meta.end is not None + + +def test_context_manager(tmp_path: Path, executor: Executor): + path = tmp_path / "my-open-folder" + + executor.init(path) + + with executor: + assert executor.meta is not None + assert executor.meta.start is not None