diff --git a/doc/source/protocols/cross_resonance/amplitude.png b/doc/source/protocols/cross_resonance/amplitude.png new file mode 100644 index 0000000000..cca5dce4c9 Binary files /dev/null and b/doc/source/protocols/cross_resonance/amplitude.png differ diff --git a/doc/source/protocols/cross_resonance/cross_resonance.rst b/doc/source/protocols/cross_resonance/cross_resonance.rst new file mode 100644 index 0000000000..10638025e1 --- /dev/null +++ b/doc/source/protocols/cross_resonance/cross_resonance.rst @@ -0,0 +1,188 @@ +Calibration of CNOT gate using Cross-Resonance +=============================================== + +It is possible to generate an interaction between two superconducting qubits without requiring +flux tunability, through a mechanism known as Cross Resonance (CR). This mechanism relies only +on microwave drive pulses. Moreover, not using flux lines, results in a reduction of the number +of fridge lines and allows to ignore all problems related to flux noise. + +The cross resonance effect was first proposed :cite:p:`CR_First` in and later +independently discovered in :cite:p:`CR_Righetti, CR_Second`. + +The CR effect can be showed by starting with the Hamiltonian of a two-qubit system with +a drive term on the first qubit :cite:p:`Manenti:2023zzn` + +.. math:: + + H = b_1^\dagger b_1 \omega_1 + \frac{\alpha_1}{2} b_1^\dagger b_1^\dagger b_1 b_1 + + b_2^\dagger b_2 \omega_2 + \frac{\alpha_2}{2} b_2^\dagger b_2^\dagger b_2 b_2 + + g (b_1 b_2^\dagger + b_1^\dagger b_2) + \Omega(t) (b_1 + b_1^\dagger) + +If we are in a dispersive regime (i.e. :math:`|\omega_1 - \omega_2| \gg g`), through a +Schrieffer-Wolff transformation we can obtain the effective Hamiltonian: + +.. math:: + + H_\text{eff} = - \frac{\tilde{\omega_1}}{2} \sigma_1^z - \frac{\tilde{\omega_2}}{2} \sigma_2^z + + \frac{\zeta}{4} \sigma_1^z \sigma_2^z + + \Omega(t) \Big[ \sigma_1^x + \nu \sigma_2^x + \mu \sigma_1^z \sigma_2^x\Big] + +where :math:`\zeta` is the ZZ coupling, :math:`\nu` is quantum crosstalk factor and :math:`\mu` is the +cross-resonance factor. From the equation above we can see that by driving the first qubit +at the frequency of the second qubit . + +By tuning the amplitude and the duration of this drive pulse it is possible to calibrate a +:math:`RZX` rotation to rotate exactly by :math:`- \pi/2`. This is done because starting +from a :math:`ZX_{frac{\pi}{2}}` we can obtain a CNOT gate using single qubit rotations. + +.. math:: + + \text{CNOT} = \text{R}_\text{ZX}(-\pi/2) \text{R}_\text{IX}(\pi/2) \text{R}_\text{ZI}(\pi/2) + +In Qibocal we provide protocols to calibrate CR pulses. + +Sweeping the duration of the CR pulse +------------------------------------- + +In a first experiment we can sweep the duration of the CR pulse and measure both the +target and control qubit. The measurement is performed while preparing the control +qubit in state :math:`\ket{0}` and :math:`\ket{1}`. + +Parameters +^^^^^^^^^^ + +.. autoclass:: qibocal.protocols.two_qubit_interaction.cross_resonance.length.CrossResonanceLengthParameters + :noindex: + +Example +^^^^^^^ + +A possible runcard to launch the experiment could be the following: + +.. code-block:: yaml + + - id: CR length + operation: cross_resonance_length + parameters: + targets: [[0,1]] + pulse_duration_start: 10 + pulse_duration_end: 200 + pulse_duration_step: 10 + flux_pulse_amplitude: 0.1 + nshots: 2000 + relaxation_time: 50000 + + +The expected output is the following: + +.. image:: length.png + +Post-processing +^^^^^^^^^^^^^^^ + +The probability of the target qubit is fitted in both cases to a dumped cosine functions. +It is possible to extract the effective coupling as + +.. math:: + + \text{J}_\text{eff}/ 2 \pi = \frac{f^{\pi}_\text{Rabi} - f_\text{Rabi}}{2} + + +where :math:`f^{\pi}_\text{Rabi}` and :math:`f_\text{Rabi}` are the frequencies of the +fitted Rabi oscillations on the target qubit. + +Sweeping amplitude of the CR pulse +---------------------------------- + +Similarly it is possible to sweep the amplitude of the CR pulse and measure both the +target and control qubit. + + +Parameters +^^^^^^^^^^ + +.. autoclass:: qibocal.protocols.two_qubit_interaction.cross_resonance.length.CrossResonanceLengthParameters + :noindex: + +Example +^^^^^^^ + +A possible runcard to launch the experiment could be the following: + +.. code-block:: yaml + + - id: CR amplitude + operation: cross_resonance_amplitude + parameters: + targets: [[0,1]] + max_amp: 0.05 + min_amp: 0.01 + step_amp: 0.005 + pulse_duration: 100 + nshots: 2000 + relaxation_time: 50000 + + +The expected output is the following: + +.. image:: amplitude.png + +Post-processing +^^^^^^^^^^^^^^^ + +The probability of the target qubit is fitted in both cases to a cosine function. + +Hamiltonian Tomography measurement +---------------------------------- + +Although from the two previous experiments it is possible to perform an initial +calibration of the CR gate, by performing a state tomography on the target qubit it is +possible to reconstruct the effective Hamiltonian of the system :cite:p:`CRDrag`: + +.. math:: + + H_\text{eff} = \frac{\nu_\text{ZX}}{2} \text{ZX} + \frac{\nu_\text{ZY}}{2} \text{ZY} + + \frac{\nu_\text{ZZ}}{2} \text{ZZ} + \frac{\nu_\text{IX}}{2} \text{IX} + + \frac{\nu_\text{IY}}{2} \text{IY} + \frac{\nu_\text{IZ}}{2} \text{IZ} + +In particular, by sweeping the duration of the CR pulse and measuring the expectation +values of the target qubit :math:`\langle X \rangle`, :math:`\langle Y \rangle` and :math:`\langle Z \rangle` +when the control qubit is prepared in :math:`\ket{0}` and :math:`\ket{1}` we can compute all terms in +the effective Hamiltonian following the procedure in :cite:p:`CRDrag`. + +Parameters +^^^^^^^^^^ + + +.. autoclass:: qibocal.protocols.two_qubit_interaction.cross_resonance.hamiltonian_tomography.length.HamiltonianTomographyCRLengthParameters + :noindex: + +Example +^^^^^^^ + +A possible runcard to launch the experiment could be the following: + +.. code-block:: yaml + + - id: Hamiltonian tomography CR + operation: cross_resonance_amplitude + parameters: + targets: [[0,1]] + nshots: 2000 + pulse_amplitude: 0.1 + pulse_duration_end: 400 + pulse_duration_start: 10 + pulse_duration_step: 20 + + +The expected output is the following: + +.. image:: tomography_length.png + + + +Requirements +^^^^^^^^^^^^ + +To run these experiments single qubit gates for both target and control qubit needs +to be calibrated. diff --git a/doc/source/protocols/cross_resonance/length.png b/doc/source/protocols/cross_resonance/length.png new file mode 100644 index 0000000000..a04dfc7595 Binary files /dev/null and b/doc/source/protocols/cross_resonance/length.png differ diff --git a/doc/source/protocols/cross_resonance/tomography_length.png b/doc/source/protocols/cross_resonance/tomography_length.png new file mode 100644 index 0000000000..229d70cbe3 Binary files /dev/null and b/doc/source/protocols/cross_resonance/tomography_length.png differ diff --git a/doc/source/protocols/index.rst b/doc/source/protocols/index.rst index 34409aa87e..a989a4a6ce 100644 --- a/doc/source/protocols/index.rst +++ b/doc/source/protocols/index.rst @@ -44,4 +44,5 @@ In this section we introduce the basics of all protocols supported by ``qibocal` coherence/index chsh twpa_calibration/twpa + cross_resonance/cross_resonance references diff --git a/doc/source/refs.bib b/doc/source/refs.bib index d8dfdd2e63..cc3e08ff1f 100644 --- a/doc/source/refs.bib +++ b/doc/source/refs.bib @@ -284,6 +284,7 @@ @misc{pasquale2024qibocalopensourceframeworkcalibration primaryClass={quant-ph}, url={https://arxiv.org/abs/2410.00101}, } + @misc{shin2024inoperandomicrowavescatteringparametercalibrated, title={In-operando microwave scattering-parameter calibrated measurement of a Josephson travelling wave parametric amplifier}, author={S. H. Shin and M. Stanley and W. N. Wong and T. Sweetnam and A. Elarabi and T. Lindström and N. M. Ridler and S. E. de Graaf}, @@ -293,3 +294,65 @@ @misc{shin2024inoperandomicrowavescatteringparametercalibrated primaryClass={quant-ph}, url={https://arxiv.org/abs/2406.03063}, } + +@article{CR_First, + title = {Microwave-induced coupling of superconducting qubits}, + author = {Paraoanu, G. S.}, + journal = {Phys. Rev. B}, + volume = {74}, + issue = {14}, + pages = {140504}, + numpages = {4}, + year = {2006}, + month = {Oct}, + publisher = {American Physical Society}, + doi = {10.1103/PhysRevB.74.140504}, + url = {https://link.aps.org/doi/10.1103/PhysRevB.74.140504} +} + +@article{CR_Righetti, + title = {Fully microwave-tunable universal gates in superconducting qubits with linear couplings and fixed transition frequencies}, + author = {Rigetti, Chad and Devoret, Michel}, + journal = {Phys. Rev. B}, + volume = {81}, + issue = {13}, + pages = {134507}, + numpages = {7}, + year = {2010}, + month = {Apr}, + publisher = {American Physical Society}, + doi = {10.1103/PhysRevB.81.134507}, + url = {https://link.aps.org/doi/10.1103/PhysRevB.81.134507} +} + +@article{CR_Second, +author = {Groot, P. and Lisenfeld, Jürgen and Schouten, R. and Ashhab, S. and Lupascu, Adrian and Harmans, Kees and Mooij, Hans}, +year = {2010}, +month = {08}, +pages = {}, +title = {Selective darkening of degenerate transitions demonstrated with two +superconducting quantum bits}, +volume = {6}, +journal = {Nature Physics}, +doi = {10.1038/nphys1733} +} + +@book{Manenti:2023zzn, + author = "Manenti, Riccardo and Motta, Mario", + title = "{Quantum Information Science}", + isbn = "978-0-19-878748-8", + publisher = "Oxford University Press", + month = "8", + year = "2023" +} + +@article{CRDrag, +author = {Li, Boxi and Calarco, Tommaso and Motzoi, Felix}, +year = {2024}, +month = {07}, +pages = {}, +title = {Experimental error suppression in Cross-Resonance gates via multi-derivative pulse shaping}, +volume = {10}, +journal = {npj Quantum Information}, +doi = {10.1038/s41534-024-00863-4} +} diff --git a/platforms/qutrits/calibration.json b/platforms/qutrits/calibration.json new file mode 100644 index 0000000000..def074c85e --- /dev/null +++ b/platforms/qutrits/calibration.json @@ -0,0 +1,75 @@ +{ + "single_qubits": { + "0": { + "resonator": { + "bare_frequency": 0.0, + "dressed_frequency": 0.0, + "depletion_time": 0, + "bare_frequency_amplitude": null + }, + "qubit": { + "frequency_01": 5114000000.0, + "frequency_12": 4784000000.0, + "maximum_frequency": 5114000000.0, + "asymmetry": 0.0, + "sweetspot": 0.0, + "flux_coefficients": null + }, + "readout": { + "fidelity": 0.0, + "coupling": null, + "effective_temperature": null, + "ground_state": [ + 0.0, + 1.0 + ], + "excited_state": [ + 1.0, + 0.0 + ], + "qudits_frequency": {} + }, + "t1": null, + "t2": null, + "t2_spin_echo": null, + "rb_fidelity": null + }, + "1": { + "resonator": { + "bare_frequency": 0.0, + "dressed_frequency": 0.0, + "depletion_time": 0, + "bare_frequency_amplitude": null + }, + "qubit": { + "frequency_01": 4914000000.0, + "frequency_12": 4584000000.0, + "maximum_frequency": 4914000000.0, + "asymmetry": 0.0, + "sweetspot": 0.0, + "flux_coefficients": null + }, + "readout": { + "fidelity": 0.0, + "coupling": null, + "effective_temperature": null, + "ground_state": [ + 0.0, + 1.0 + ], + "excited_state": [ + 1.0, + 0.0 + ], + "qudits_frequency": {} + }, + "t1": null, + "t2": null, + "t2_spin_echo": null, + "rb_fidelity": null + } + }, + "two_qubits": {}, + "readout_mitigation_matrix": null, + "flux_crosstalk_matrix": null +} diff --git a/platforms/qutrits/parameters.json b/platforms/qutrits/parameters.json new file mode 100644 index 0000000000..fdd91e5840 --- /dev/null +++ b/platforms/qutrits/parameters.json @@ -0,0 +1,357 @@ +{ + "settings": { + "nshots": 1024, + "relaxation_time": 0 + }, + "configs": { + "emulator/bounds": { + "kind": "bounds", + "waveforms": 1000000.0, + "readout": 50, + "instructions": 200 + }, + "hamiltonian": { + "kind": "hamiltonian", + "transmon_levels": 3, + "single_qubit": { + "0": { + "frequency": 5114000000.0, + "anharmonicity": -330000000.0, + "t1": { + "0-1": 38000.0, + "1-2": 19000.0 + }, + "t2": { + "0-1": 50000.0 + } + }, + "1": { + "frequency": 4914000000.0, + "anharmonicity": -330000000.0, + "t1": { + "0-1": 42000.0, + "1-2": 21000.0 + }, + "t2": { + "0-1": 61000.0 + } + } + }, + "two_qubit": { + "0-1": { + "coupling": 3800000.0 + } + } + }, + "0/drive": { + "kind": "drive-emulator", + "frequency": 5114000000.0, + "rabi_frequency": 50000000.0, + "scale_factor": 10.0 + }, + "0/drive1": { + "kind": "drive-emulator", + "frequency": 4914000000.0, + "rabi_frequency": 50000000.0, + "scale_factor": 10.0 + }, + "0/drive12": { + "kind": "drive-emulator", + "frequency": 4784000000.0, + "rabi_frequency": 50000000.0, + "scale_factor": 10.0 + }, + "1/drive": { + "kind": "drive-emulator", + "frequency": 4914000000.0, + "rabi_frequency": 50000000.0, + "scale_factor": 10.0 + }, + "1/drive12": { + "kind": "drive-emulator", + "frequency": 4584000000.0, + "rabi_frequency": 50000000.0, + "scale_factor": 10.0 + }, + "0/probe": { + "kind": "iq", + "frequency": 5500000000.0 + }, + "0/acquisition": { + "kind": "acquisition", + "delay": 0.0, + "smearing": 0.0, + "threshold": 0.0, + "iq_angle": 0.0, + "kernel": null + }, + "1/probe": { + "kind": "iq", + "frequency": 5500000000.0 + }, + "1/acquisition": { + "kind": "acquisition", + "delay": 0.0, + "smearing": 0.0, + "threshold": 0.0, + "iq_angle": 0.0, + "kernel": null + } + }, + "native_gates": { + "single_qubit": { + "0": { + "RX": [ + [ + "0/drive", + { + "kind": "pulse", + "duration": 20.0, + "amplitude": 0.200449, + "envelope": { + "kind": "drag", + "rel_sigma": 0.1, + "beta": 0.25 + }, + "relative_phase": 0.0 + } + ] + ], + "RX90": [ + [ + "0/drive", + { + "kind": "pulse", + "duration": 20.0, + "amplitude": 0.1002245, + "envelope": { + "kind": "drag", + "rel_sigma": 0.1, + "beta": 0.25 + }, + "relative_phase": 0.0 + } + ] + ], + "RX12": [ + [ + "0/drive12", + { + "kind": "pulse", + "duration": 20.0, + "amplitude": 0.1426586, + "envelope": { + "kind": "gaussian", + "rel_sigma": 0.1 + }, + "relative_phase": 0.0 + } + ] + ], + "MZ": [ + [ + "0/acquisition", + { + "kind": "readout", + "acquisition": { + "kind": "acquisition", + "duration": 10.0 + }, + "probe": { + "kind": "pulse", + "duration": 10.0, + "amplitude": 0.1, + "envelope": { + "kind": "rectangular" + }, + "relative_phase": 0.0 + } + } + ] + ], + "CP": null + }, + "1": { + "RX": [ + [ + "1/drive", + { + "kind": "pulse", + "duration": 20.0, + "amplitude": 0.200792, + "envelope": { + "kind": "drag", + "rel_sigma": 0.1, + "beta": 0.25 + }, + "relative_phase": 0.0 + } + ] + ], + "RX90": [ + [ + "1/drive", + { + "kind": "pulse", + "duration": 20.0, + "amplitude": 0.100396, + "envelope": { + "kind": "drag", + "rel_sigma": 0.1, + "beta": 0.25 + }, + "relative_phase": 0.0 + } + ] + ], + "RX12": [ + [ + "1/drive12", + { + "kind": "pulse", + "duration": 20.0, + "amplitude": 0.138545, + "envelope": { + "kind": "gaussian", + "rel_sigma": 0.1 + }, + "relative_phase": 0.0 + } + ] + ], + "MZ": [ + [ + "1/acquisition", + { + "kind": "readout", + "acquisition": { + "kind": "acquisition", + "duration": 10.0 + }, + "probe": { + "kind": "pulse", + "duration": 10.0, + "amplitude": 0.1, + "envelope": { + "kind": "rectangular" + }, + "relative_phase": 0.0 + } + } + ] + ], + "CP": null + } + }, + "coupler": {}, + "two_qubit": { + "0-1": { + "CZ": null, + "CNOT": [ + [ + "0/drive1", + { + "kind": "pulse", + "duration": 50.0, + "amplitude": 0.1, + "envelope": { + "kind": "rectangular" + }, + "relative_phase": 0.0 + } + ], + [ + "0/drive", + { + "kind": "delay", + "duration": 50.0 + } + ], + [ + "0/drive", + { + "kind": "pulse", + "duration": 20.0, + "amplitude": 0.200449, + "envelope": { + "kind": "drag", + "rel_sigma": 0.1, + "beta": 0.25 + }, + "relative_phase": 0.0 + } + ], + [ + "0/drive1", + { + "kind": "delay", + "duration": 20.0 + } + ], + [ + "01/drive", + { + "kind": "pulse", + "duration": 50.0, + "amplitude": -0.1, + "envelope": { + "kind": "rectangular" + }, + "relative_phase": 0.0 + } + ], + [ + "0/drive", + { + "kind": "delay", + "duration": 50.0 + } + ], + [ + "0/drive", + { + "kind": "pulse", + "duration": 20.0, + "amplitude": 0.200449, + "envelope": { + "kind": "drag", + "rel_sigma": 0.1, + "beta": 0.25 + }, + "relative_phase": 0.0 + } + ], + [ + "1/drive", + { + "kind": "delay", + "duration": 140.0 + } + ], + [ + "1/drive", + { + "kind": "pulse", + "duration": 20.0, + "amplitude": 0.100396, + "envelope": { + "kind": "drag", + "rel_sigma": 0.1, + "beta": 0.25 + }, + "relative_phase": 0.0 + } + ], + [ + "0/drive", + { + "kind": "virtualz", + "phase": -1.5707963267948966 + } + ] + ], + "iSWAP": null + } + } + } +} diff --git a/platforms/qutrits/platform.py b/platforms/qutrits/platform.py new file mode 100644 index 0000000000..00dc63d7b2 --- /dev/null +++ b/platforms/qutrits/platform.py @@ -0,0 +1,45 @@ +import pathlib + +from qibolab import ConfigKinds +from qibolab._core.components import IqChannel +from qibolab._core.instruments.emulator.emulator import EmulatorController +from qibolab._core.instruments.emulator.hamiltonians import ( + DriveEmulatorConfig, + HamiltonianConfig, +) +from qibolab._core.platform import Platform +from qibolab._core.qubits import Qubit + +FOLDER = pathlib.Path(__file__).parent + +ConfigKinds.extend([HamiltonianConfig, DriveEmulatorConfig]) + + +def create() -> Platform: + """Create a dummy platform using the dummy instrument.""" + qubits = {} + channels = {} + + qubits[0] = qubit = Qubit.default( + 0, drive_extra={(1, 2): "0/drive12", 1: "0/drive1"} + ) + channels |= { + qubit.drive: IqChannel(mixer=None, lo=None), + qubits[0].drive_extra[1, 2]: IqChannel(mixer=None, lo=None), + qubits[0].drive_extra[1]: IqChannel(mixer=None, lo=None), + } + qubits[1] = qubit = Qubit.default(1, drive_extra={(1, 2): "1/drive12"}) + channels |= { + qubit.drive: IqChannel(mixer=None, lo=None), + qubits[1].drive_extra[1, 2]: IqChannel(mixer=None, lo=None), + } + # register the instruments + instruments = { + "dummy": EmulatorController(address="0.0.0.0", channels=channels), + } + + return Platform.load( + path=FOLDER, + instruments=instruments, + qubits=qubits, + ) diff --git a/src/qibocal/auto/transpile.py b/src/qibocal/auto/transpile.py index fac8f20067..c426f9784f 100644 --- a/src/qibocal/auto/transpile.py +++ b/src/qibocal/auto/transpile.py @@ -297,7 +297,11 @@ def _natives(platform: Platform) -> dict[str, NativeContainer]: # add two qubit natives only if there are pairs pair = next(iter(platform.pairs)) two_qubit_natives_container = platform.natives.two_qubit[pair] - two_qubit_natives = list(two_qubit_natives_container.model_fields) + two_qubit_natives = [ + i + for i in list(two_qubit_natives_container.model_fields) + if getattr(two_qubit_natives_container, i) is not None + ] else: two_qubit_natives = [] # Solve Qibo-Qibolab mismatch diff --git a/src/qibocal/protocols/rabi/amplitude.py b/src/qibocal/protocols/rabi/amplitude.py index 275c6d52f4..e30d7e00de 100644 --- a/src/qibocal/protocols/rabi/amplitude.py +++ b/src/qibocal/protocols/rabi/amplitude.py @@ -118,7 +118,7 @@ def _fit(data: RabiAmplitudeData) -> RabiAmplitudeResults: signal=False, ) pi_pulse_amplitudes[qubit] = [pi_pulse_parameter, perr[2] / 2] - fitted_parameters[qubit] = popt.tolist() + fitted_parameters[qubit] = popt durations = {key: [value, 0] for key, value in data.durations.items()} chi2[qubit] = [ chi2_reduced( diff --git a/src/qibocal/protocols/rabi/utils.py b/src/qibocal/protocols/rabi/utils.py index a49cb2b33a..047e1b7e21 100644 --- a/src/qibocal/protocols/rabi/utils.py +++ b/src/qibocal/protocols/rabi/utils.py @@ -5,17 +5,15 @@ from scipy.optimize import curve_fit from qibocal.auto.operation import Parameters, QubitId -from qibocal.update import replace - -from ..utils import ( +from qibocal.protocols.utils import ( COLORBAND, COLORBAND_LINE, angle_wrap, - fallback_period, guess_period, table_dict, table_html, ) +from qibocal.update import replace QUANTILE_CONSTANT = 1.5 """Scaling factor to recover signal amplitude from quantiles. @@ -53,7 +51,7 @@ def rabi_length_function(x, offset, amplitude, period, phase, t2_inv): def rabi_initial_guess(x, y, experiment: str, signal: bool): - period = fallback_period(guess_period(x, y)) + period = guess_period(x, y) median_sig = np.median(y) q80 = np.quantile(y, 0.8) q20 = np.quantile(y, 0.2) @@ -338,7 +336,7 @@ def sequence_length( def fit_length_function( x, y, guess, sigma=None, signal=True, x_limits=(None, None), y_limits=(None, None) -): +) -> tuple[list[float], list[float], float]: popt, perr = curve_fit( rabi_length_function, x, @@ -374,12 +372,12 @@ def fit_length_function( ] pi_pulse_parameter = popt[2] / 2 * period_correction_factor(phase=popt[3]) - return popt, perr, pi_pulse_parameter + return popt, perr.tolist(), pi_pulse_parameter def fit_amplitude_function( x, y, guess, sigma=None, signal=True, x_limits=(None, None), y_limits=(None, None) -): +) -> tuple[list[float], list[float], float]: popt, perr = curve_fit( rabi_amplitude_function, x, @@ -404,5 +402,9 @@ def fit_amplitude_function( - 2 * np.pi * x_limits[0] / (x_limits[1] - x_limits[0]) / popt[2] ), ] + else: + popt = popt.tolist() + pi_pulse_parameter = popt[2] / 2 * period_correction_factor(phase=popt[3]) - return popt, perr, pi_pulse_parameter + + return popt, perr.tolist(), pi_pulse_parameter diff --git a/src/qibocal/protocols/ramsey/processing.py b/src/qibocal/protocols/ramsey/processing.py index e7a187e4fa..546d051fab 100644 --- a/src/qibocal/protocols/ramsey/processing.py +++ b/src/qibocal/protocols/ramsey/processing.py @@ -15,8 +15,7 @@ from qibocal.protocols.utils import ( GHZ_TO_HZ, angle_wrap, - fallback_period, - guess_period, + quinn_fernandes_algorithm, table_dict, table_html, ) @@ -77,8 +76,7 @@ def fitting(x: list, y: list) -> tuple[list[float], list[float]]: y = (y - y_min) / delta_y x = (x - x_min) / delta_x - period = fallback_period(guess_period(x, y)) - omega = 2 * np.pi / period + omega = quinn_fernandes_algorithm(y, x, speedup_flag=True) median_sig = np.median(y) q80 = np.quantile(y, 0.8) q20 = np.quantile(y, 0.2) diff --git a/src/qibocal/protocols/two_qubit_interaction/__init__.py b/src/qibocal/protocols/two_qubit_interaction/__init__.py index 3997d14787..4918f768f2 100644 --- a/src/qibocal/protocols/two_qubit_interaction/__init__.py +++ b/src/qibocal/protocols/two_qubit_interaction/__init__.py @@ -1,5 +1,11 @@ from .chevron import chevron, chevron_signal from .chsh import chsh +from .cross_resonance import ( + cancellation_amplitude_tuning, + cancellation_phase_tuning, + cr_amplitude, + cr_length, +) from .optimize import optimize_two_qubit_gate from .snz_optimize import snz_optimize from .snz_optimize_t_idle import snz_optimize_t_idle @@ -15,4 +21,8 @@ "snz_optimize_t_idle", "snz_optimize_t_idle_vs_t_tot", "chsh", + "cr_length", + "cr_amplitude", + "cancellation_amplitude_tuning", + "cancellation_phase_tuning", ] diff --git a/src/qibocal/protocols/two_qubit_interaction/cross_resonance/__init__.py b/src/qibocal/protocols/two_qubit_interaction/cross_resonance/__init__.py new file mode 100644 index 0000000000..3b3b30ed72 --- /dev/null +++ b/src/qibocal/protocols/two_qubit_interaction/cross_resonance/__init__.py @@ -0,0 +1,11 @@ +from .cancellation_amplitude import cancellation_amplitude_tuning +from .cancellation_phase import cancellation_phase_tuning +from .cross_resonance_amplitude import cr_amplitude +from .cross_resonance_length import cr_length + +__all__ = [ + "cr_length", + "cr_amplitude", + "cancellation_phase_tuning", + "cancellation_amplitude_tuning", +] diff --git a/src/qibocal/protocols/two_qubit_interaction/cross_resonance/cancellation_amplitude.py b/src/qibocal/protocols/two_qubit_interaction/cross_resonance/cancellation_amplitude.py new file mode 100644 index 0000000000..eddc581df7 --- /dev/null +++ b/src/qibocal/protocols/two_qubit_interaction/cross_resonance/cancellation_amplitude.py @@ -0,0 +1,437 @@ +"""Hamiltonian tomography protocol for CR gate calibration. + +This protocol computes the expectation values for X, Y and Z for the target qubit +after the application of a cross resonance sequence. The CR pulses are played on the control drive +channel with frequency set to the frequency of the target drive channel. +""" + +from dataclasses import dataclass, field + +import numpy as np +from qibolab import ( + AcquisitionType, + AveragingMode, + ParallelSweepers, + Parameter, + Sweeper, +) +from scipy.constants import kilo + +from qibocal.auto.operation import ( + QubitId, + QubitPairId, + Routine, +) +from qibocal.calibration import CalibrationPlatform +from qibocal.protocols.utils import table_dict, table_html + +from .cr_parent_classes import ( + Basis, + HamiltonianTerm, + HamiltonianTomographyData, + HamiltonianTomographyParameters, + HamiltonianTomographyResults, + SetControl, + check_qubit_overlap, +) +from .cross_resonance_processing import ( + cancellation_amplitude_fit, + reconstruct_full_hamiltonian_terms, +) +from .plotting import cancellation_calibration_plot +from .utils import ( + cross_resonance_experiment, + cross_resonance_pulses, + ro_delay_range, + update_cnot_from_fit, +) + +__all__ = ["cancellation_amplitude_tuning"] + +HamiltonianTomographyCANCAmplType = np.dtype( + [ + ("prob_target", np.float64), + ("error_target", np.float64), + ("prob_control", np.float64), + ("error_control", np.float64), + ("amp", np.float64), + ("x", np.float64), + ] +) +"""Custom dtype for Cancellation amplitude.""" + + +@dataclass(kw_only=True) +class HamiltonianTomographyCANCAmplParameters(HamiltonianTomographyParameters): + """HamiltonianTomographyCANCAmplitude runcard inputs.""" + + target_ampl_range: tuple[float, float, float] + """Amplitude range of cancellation pulse.""" + verbose_plot: bool = False + """If `True` in the report all the single Hamiltonian tomographies are plotted.""" + + +@dataclass(kw_only=True) +class HamiltonianTomographyCANCAmplResults(HamiltonianTomographyResults): + """HamiltonianTomographyCANCAmpl outputs.""" + + hamiltonian_terms: dict[ + tuple[QubitId, QubitId], list[tuple[float, dict[HamiltonianTerm, float]]] + ] = field(default_factory=dict) + """Terms in effective Hamiltonian.""" + + fitted_parameters: dict[tuple[QubitId, QubitId], dict[HamiltonianTerm, list]] = ( + field(default_factory=dict) + ) + """Fitted parameters for Hamiltonian Terms values for different amplitudes.""" + + cancellation_pulse_amplitudes: dict[QubitPairId, dict[str, float]] = field( + default_factory=dict + ) + """Fitted parameters for cancellation pulse amplitudes.""" + + hamiltonian_tom_params: dict[ + float, dict[tuple[QubitId, QubitId, SetControl], list[float]] + ] = field(default_factory=dict) + """Fitted parameters for of Hamiltonian Tomography experiment for each qubit and per amplitude value. + Used for plotting with the `verbose_plot` option""" + + cr_lengths: dict[float, dict[tuple[QubitId, QubitId, SetControl], list[float]]] = ( + field(default_factory=dict) + ) + """Cross resonance pulse duration found in Hamiltonian Tomography experiment for each qubit and per amplitude value. + Used for plotting with the `verbose_plot` option""" + + verbose_plot: bool = False + """If `True` in the report all the single Hamiltonian tomographies are plotted.""" + + def select_pair_and_ampl_ham_params(self, amplitude: float, pair: QubitPairId): + """Data + Select and refactor Hamiltonian tomography parameters for a given amplitude and qubit pair. + """ + selected_ham_tom_params = next( + ( + val_ham_params + for (amp, val_ham_params) in self.hamiltonian_terms[pair] + if amp == amplitude + ), + None, + ) + + return reconstruct_full_hamiltonian_terms(selected_ham_tom_params, pair) + + +@dataclass(kw_only=True) +class HamiltonianTomographyCANCAmplData(HamiltonianTomographyData): + """Data structure for CANC Amplitude.""" + + verbose_plot: bool = False + """If `True` in the report all the single Hamiltonian tomographies are plotted.""" + + @property + def amplitudes(self) -> list[float]: + first_key = next(iter(self.data.keys())) + return np.unique(self.data[first_key].amp).tolist() + + def select_amplitude(self, amplitude: float) -> HamiltonianTomographyData: + new_data = HamiltonianTomographyData(echo=self.echo) + new_data.data = {k: d[d.amp == amplitude] for k, d in self.data.items()} + return new_data + + def register_qubit(self, dtype, data_keys, data_dict): + """Store output for single qubit.""" + duration_list = data_dict["x"] + amp_list = data_dict["amp"] + size = len(duration_list) * len(amp_list) + ar = np.empty(size, dtype=dtype) + amplitudes, durations = np.meshgrid(amp_list, duration_list) + ar["x"] = durations.ravel() + ar["amp"] = amplitudes.ravel() + ar["prob_target"] = data_dict["prob_target"].ravel() + ar["error_target"] = data_dict["error_target"].ravel() + ar["prob_control"] = data_dict["prob_control"].ravel() + ar["error_control"] = data_dict["error_control"].ravel() + + self.data[data_keys] = np.rec.array(ar) + + +def _acquisition( + params: HamiltonianTomographyCANCAmplParameters, + platform: CalibrationPlatform, + targets: list[QubitPairId], +) -> HamiltonianTomographyCANCAmplData: + """Data acquisition for Hamiltonian tomography CR protocol. + + We measure the expectation values X,Y and Z on the target qubit after + applying the CR sequence specified by the input parameters. We repeat the + measurement twice for each target qubit, once with the control qubit in state 0 + and once with the control qubit in state 1. + + We store the probability of the control qubit and the expectation value of the target qubit. + + """ + + data = HamiltonianTomographyCANCAmplData( + echo=params.echo, + verbose_plot=params.verbose_plot, + ) + + # check validity of input + check_qubit_overlap(targets) + + updates = [] + control_ampls: dict[QubitPairId, float] = {} + control_phases: dict[QubitPairId, float] = {} + target_phases: dict[QubitPairId, float] = {} + for pair in targets: + control, target = pair + + updates.append( + { + platform.qubits[control].drive_extra[target]: { + "frequency": platform.config( + platform.qubits[target].drive + ).frequency + } + } + ) + + cr_pulse, canc_pulse = cross_resonance_pulses(platform, control, target) + control_ampls |= {pair: cr_pulse.amplitude} + control_phases |= {pair: cr_pulse.relative_phase} + target_phases |= {pair: canc_pulse.relative_phase} + + for basis in Basis: + for setup in SetControl: + sequence, cr_pulses, cr_target_pulses, cr_delays, ro_delays = ( + cross_resonance_experiment( + platform=platform, + pair_list=targets, + duration=0.0, + ctrl_ampl=control_ampls, + ctrl_phase=control_phases, + targ_ampl=0.0, + targ_phase=target_phases, + basis=basis, + setup=setup, + echo=params.echo, + interpolated_sweeper=params.interpolated_sweeper, + ) + ) + + if params.interpolated_sweeper: + duration_parallel_sweeper = ParallelSweepers( + [ + Sweeper( + parameter=Parameter.duration_interpolated, + range=params.duration_range, + pulses=cr_pulses[pair] + cr_target_pulses[pair], + ) + for pair in targets + ] + ) + else: + duration_parallel_sweeper = ParallelSweepers( + [ + Sweeper( + parameter=Parameter.duration, + range=params.duration_range, + pulses=cr_pulses[pair] + + cr_target_pulses[pair] + + cr_delays[pair], + ) + for pair in targets + ] + ) + ro_sweeper = ParallelSweepers( + [ + Sweeper( + parameter=Parameter.duration, + range=ro_delay_range( + cr_pulse_duration_range=params.duration_range, + echo=params.echo, + cntl_setup=setup, + control=c, + platform=platform, + ), + pulses=ro_delays[(c, t)], + ) + for c, t in targets + ] + ) + duration_parallel_sweeper += ro_sweeper + + amp_sweepers = ParallelSweepers( + [ + Sweeper( + parameter=Parameter.amplitude, + range=params.target_ampl_range, + pulses=[cr_target_pulses[pair][0] for pair in targets], + ) + ] + ) + if params.echo: + # sweeping over the out-of-phase signal (refocusing) of the cancellation pulse + echo_ampl_range = tuple(-x for x in params.target_ampl_range) + amp_sweepers += ParallelSweepers( + [ + Sweeper( + parameter=Parameter.amplitude, + range=echo_ampl_range, + pulses=[cr_target_pulses[pair][1] for pair in targets], + ) + ] + ) + + results = platform.execute( + [sequence], + [ + duration_parallel_sweeper, + amp_sweepers, + ], + nshots=params.nshots, + relaxation_time=params.relaxation_time, + acquisition_type=AcquisitionType.DISCRIMINATION, + averaging_mode=AveragingMode.CYCLIC, + updates=updates, + ) + + for ctrl, targ in targets: + target_acq_handle = list( + sequence.channel(platform.qubits[targ].acquisition) + )[-1].id + control_acq_handle = list( + sequence.channel(platform.qubits[ctrl].acquisition) + )[-1].id + + prob_target = results[target_acq_handle].ravel() + prob_control = results[control_acq_handle].ravel() + + data.register_qubit( + HamiltonianTomographyCANCAmplType, + (ctrl, targ, basis, setup), + dict( + x=np.arange(*params.duration_range), + amp=np.arange(*params.target_ampl_range), + prob_target=1 - 2 * prob_target, + error_target=2 + * np.sqrt(prob_target * (1 - prob_target) / params.nshots), + prob_control=1 - 2 * prob_control, + error_control=2 + * np.sqrt(prob_control * (1 - prob_control) / params.nshots), + ), + ) + + return data + + +def _fit( + data: HamiltonianTomographyCANCAmplData, +) -> HamiltonianTomographyCANCAmplResults: + """Post-processing function for HamiltonianTomographyCANCAmpl. + + We fit the expectation values using the Eq. S10 from the paper https://arxiv.org/pdf/2303.01427. + Afterwards, we extract the Hamiltonian terms from the fitted parameters. + + """ + hamiltonian_terms, fitted_parameters, cal_amplitudes, ham_tom_params, cr_lengths = ( + cancellation_amplitude_fit( + data=data, + ) + ) + + return HamiltonianTomographyCANCAmplResults( + echo=data.echo, + hamiltonian_terms=hamiltonian_terms, + fitted_parameters=fitted_parameters, + cancellation_pulse_amplitudes=cal_amplitudes, + hamiltonian_tom_params=ham_tom_params, + cr_lengths=cr_lengths, + verbose_plot=data.verbose_plot, + ) + + +def _plot( + data: HamiltonianTomographyCANCAmplData, + target: QubitPairId, + fit: HamiltonianTomographyCANCAmplResults, +): + """Plotting function for HamiltonianTomographyCANCAmpl.""" + figs, fitting_report = cancellation_calibration_plot(data, target, fit) + + plot_ham_tom = True if fit is None else fit.verbose_plot + + if plot_ham_tom: + from .plotting import ( + tomography_cr_plot, + ) + + for a in data.amplitudes: + ampl_data = data.select_amplitude(a) + if fit is not None: + selected_ham_terms = fit.select_pair_and_ampl_ham_params(a, target) + ham_tom_fit = HamiltonianTomographyResults( + echo=fit.echo, + hamiltonian_terms=selected_ham_terms, + fitted_parameters=fit.hamiltonian_tom_params[a], + cr_lengths=fit.cr_lengths[a], + ) + + fitting_report += "\n" + table_html( + table_dict( + 8 * [target], + [f"{term.name} [MHz]" for term in HamiltonianTerm] + + ["CR duration (ns)", "Control amplitude [a.u.]"], + [ + ham_tom_fit.hamiltonian_terms[target[0], target[1], term] + * kilo + for term in HamiltonianTerm + ] + + [ + fit.cr_lengths[a][target] + if target in fit.cr_lengths[a] + else None + ] + + [a], + ) + ) + else: + ham_tom_fit = None + fitting_report = "" + + f, _ = tomography_cr_plot(ampl_data, target, ham_tom_fit) + figs += f + + return figs, fitting_report + + +def _update( + results: HamiltonianTomographyCANCAmplResults, + platform: CalibrationPlatform, + target: QubitPairId, +): + + # now no check is needed since the acquisition was executed correctly, + # which means we have all parameters defined. + cr_pulse, canc_pulse = cross_resonance_pulses(platform, target[0], target[1]) + gate_duration = cr_pulse.duration + control_amplitude = cr_pulse.amplitude + control_phase = cr_pulse.relative_phase + target_phase = canc_pulse.relative_phase + + # check if the resulting fit was succsessfull + if target in results.cancellation_pulse_amplitudes: + update_cnot_from_fit( + platform=platform, + pair=target, + cr_duration=gate_duration, + cr_ampl=control_amplitude, + control_phase=control_phase, + canc_ampl=results.cancellation_pulse_amplitudes[target]["ampl_iy"], + canc_phase=target_phase, + echo_flag=results.echo, + ) + + +cancellation_amplitude_tuning = Routine(_acquisition, _fit, _plot, _update) +"""HamiltonianTomographyCANCAmpl Routine object.""" diff --git a/src/qibocal/protocols/two_qubit_interaction/cross_resonance/cancellation_phase.py b/src/qibocal/protocols/two_qubit_interaction/cross_resonance/cancellation_phase.py new file mode 100644 index 0000000000..020327d8fe --- /dev/null +++ b/src/qibocal/protocols/two_qubit_interaction/cross_resonance/cancellation_phase.py @@ -0,0 +1,433 @@ +"""Hamiltonian tomography protocol for CR gate calibration. + +This protocol computes the expectation values for X, Y and Z for the target qubit +after the application of a cross resonance sequence. The CR pulses are played on the control drive +channel with frequency set to the frequency of the target drive channel. +""" + +from dataclasses import dataclass, field + +import numpy as np +from qibolab import ( + AcquisitionType, + AveragingMode, + ParallelSweepers, + Parameter, + Sweeper, +) +from scipy.constants import kilo + +from qibocal.auto.operation import ( + QubitId, + QubitPairId, + Routine, +) +from qibocal.calibration import CalibrationPlatform +from qibocal.protocols.utils import table_dict, table_html + +from .cr_parent_classes import ( + Basis, + HamiltonianTerm, + HamiltonianTomographyData, + HamiltonianTomographyParameters, + HamiltonianTomographyResults, + SetControl, + check_qubit_overlap, +) +from .cross_resonance_processing import ( + cancellation_phase_fit, + reconstruct_full_hamiltonian_terms, +) +from .plotting import cancellation_calibration_plot +from .utils import ( + cross_resonance_experiment, + cross_resonance_pulses, + ro_delay_range, + update_cnot_from_fit, +) + +__all__ = ["cancellation_phase_tuning"] + +HamiltonianTomographyCANCPhaseType = np.dtype( + [ + ("prob_target", np.float64), + ("error_target", np.float64), + ("prob_control", np.float64), + ("error_control", np.float64), + ("phase", np.float64), + ("x", np.int64), + ] +) +"""Custom dtype for Cancellation phase.""" + + +@dataclass(kw_only=True) +class HamiltonianTomographyCANCPhaseParameters(HamiltonianTomographyParameters): + """HamiltonianTomographyCANCPhase runcard inputs.""" + + phase_range: tuple[float, float, float] + """CR pulse phase range for control qubit.""" + verbose_plot: bool = False + """If `True` in the report all the single Hamiltonian tomographies are plotted.""" + + +@dataclass(kw_only=True) +class HamiltonianTomographyCANCPhaseResults(HamiltonianTomographyResults): + """HamiltonianTomographyCANCPhase outputs.""" + + hamiltonian_terms: dict[ + QubitPairId, list[tuple[float, dict[HamiltonianTerm, float]]] + ] = field(default_factory=dict) + """Terms in effective Hamiltonian.""" + fitted_parameters: dict[tuple[QubitId, QubitId], dict[HamiltonianTerm, list]] = ( + field(default_factory=dict) + ) + """Fitted parameters from Hamiltonian Terms values for different phases.""" + + cancellation_pulse_phases: dict[QubitPairId, dict[str, float]] = field( + default_factory=dict + ) + """Fitted parameters for cancellation pulse phases.""" + hamiltonian_tom_params: dict[ + float, dict[tuple[QubitId, QubitId, SetControl], list[float]] + ] = field(default_factory=dict) + """Fitted parameters for of Hamiltonian Tomography experiment for each qubit and per amplitude value. + Used for plotting with the `verbose_plot` option""" + + cr_lengths: dict[float, dict[tuple[QubitId, QubitId, SetControl], list[float]]] = ( + field(default_factory=dict) + ) + """Cross resonance pulse duration found in Hamiltonian Tomography experiment for each qubit and per amplitude value. + Used for plotting with the `verbose_plot` option""" + + verbose_plot: bool = False + """If `True` in the report all the single Hamiltonian tomographies are plotted.""" + + def select_pair_and_phase_ham_params(self, phase: float, pair: QubitPairId): + """Data + Select and refactor Hamiltonian tomography parameters for a given amplitude and qubit pair. + """ + selected_ham_tom_params = next( + ( + val_ham_params + for (phi, val_ham_params) in self.hamiltonian_terms[pair] + if phi == phase + ), + None, + ) + + return reconstruct_full_hamiltonian_terms(selected_ham_tom_params, pair) + + +@dataclass +class HamiltonianTomographyCANCPhaseData(HamiltonianTomographyData): + """Data structure for CANC Phase.""" + + verbose_plot: bool = False + """If `True` in the report all the single Hamiltonian tomographies are plotted.""" + + @property + def phases(self) -> list[float]: + first_key = next(iter(self.data.keys())) + return np.unique(self.data[first_key].phase).tolist() + + def select_phase(self, phase: float) -> HamiltonianTomographyData: + new_data = HamiltonianTomographyData(echo=self.echo) + new_data.data = {k: d[d.phase == phase] for k, d in self.data.items()} + return new_data + + def register_qubit(self, dtype, data_keys, data_dict): + """Store output for single qubit.""" + duration_list = data_dict["x"] + phase_list = data_dict["phase"] + size = len(duration_list) * len(phase_list) + ar = np.empty(size, dtype=dtype) + phases, durations = np.meshgrid(phase_list, duration_list) + ar["x"] = durations.ravel() + ar["phase"] = phases.ravel() + ar["prob_target"] = data_dict["prob_target"].ravel() + ar["error_target"] = data_dict["error_target"].ravel() + ar["prob_control"] = data_dict["prob_control"].ravel() + ar["error_control"] = data_dict["error_control"].ravel() + + self.data[data_keys] = np.rec.array(ar) + + +def _acquisition( + params: HamiltonianTomographyCANCPhaseParameters, + platform: CalibrationPlatform, + targets: list[QubitPairId], +) -> HamiltonianTomographyCANCPhaseData: + """Data acquisition for Hamiltonian tomography CR ctrl_phaseprotocol. + + We measure the expectation values X,Y and Z on the target qubit after + applying the CR sequence specified by the input parameters. We repeat the + measurement twice for each target qubit, once with the control qubit in state 0 + and once with the control qubit in state 1. + + We store the probability of the control qubit and the expectation value of the target qubit. + + """ + + data = HamiltonianTomographyCANCPhaseData( + echo=params.echo, + verbose_plot=params.verbose_plot, + ) + + # check validity of input + check_qubit_overlap(targets) + + updates = [] + control_ampls: dict[QubitPairId, float] = {} + target_ampls: dict[QubitPairId, float] = {} + for pair in targets: + control, target = pair + + updates.append( + { + platform.qubits[control].drive_extra[target]: { + "frequency": platform.config( + platform.qubits[target].drive + ).frequency + } + } + ) + + cr_pulse, canc_pulse = cross_resonance_pulses(platform, control, target) + if cr_pulse is None: + raise ValueError( + "CR pulse not found for control {control} and target {target}. " + "Please check the CR pulse configuration or first run previous protocols." + ) + control_ampls |= {pair: cr_pulse.amplitude} + target_ampls |= {pair: canc_pulse.amplitude if canc_pulse is not None else 0.0} + + for basis in Basis: + for setup in SetControl: + sequence, cr_pulses, cr_target_pulses, cr_delays, ro_delays = ( + cross_resonance_experiment( + platform=platform, + pair_list=targets, + duration=0.0, + ctrl_ampl=control_ampls, + ctrl_phase=0.0, + targ_ampl=target_ampls, + targ_phase=0.0, + basis=basis, + setup=setup, + echo=params.echo, + interpolated_sweeper=params.interpolated_sweeper, + ) + ) + + if params.interpolated_sweeper: + duration_parallel_sweeper = ParallelSweepers( + [ + Sweeper( + parameter=Parameter.duration_interpolated, + range=params.duration_range, + pulses=cr_pulses[pair] + cr_target_pulses[pair], + ) + for pair in targets + ] + ) + else: + duration_parallel_sweeper = ParallelSweepers( + [ + Sweeper( + parameter=Parameter.duration, + range=params.duration_range, + pulses=cr_pulses[pair] + + cr_target_pulses[pair] + + cr_delays[pair], + ) + for pair in targets + ] + ) + ro_sweeper = ParallelSweepers( + [ + Sweeper( + parameter=Parameter.duration, + range=ro_delay_range( + cr_pulse_duration_range=params.duration_range, + echo=params.echo, + cntl_setup=setup, + control=c, + platform=platform, + ), + pulses=ro_delays[(c, t)], + ) + for c, t in targets + ] + ) + duration_parallel_sweeper += ro_sweeper + + phase_sweepers = ParallelSweepers( + [ + Sweeper( + parameter=Parameter.relative_phase, + range=params.phase_range, + pulses=cr_pulses[pair], + ) + for pair in targets + ] + ) + + results = platform.execute( + [sequence], + [ + duration_parallel_sweeper, + phase_sweepers, + ], + nshots=params.nshots, + relaxation_time=params.relaxation_time, + acquisition_type=AcquisitionType.DISCRIMINATION, + averaging_mode=AveragingMode.CYCLIC, + updates=updates, + ) + + for ctrl, targ in targets: + target_acq_handle = list( + sequence.channel(platform.qubits[targ].acquisition) + )[-1].id + control_acq_handle = list( + sequence.channel(platform.qubits[ctrl].acquisition) + )[-1].id + + prob_target = results[target_acq_handle].ravel() + prob_control = results[control_acq_handle].ravel() + + data.register_qubit( + HamiltonianTomographyCANCPhaseType, + (ctrl, targ, basis, setup), + dict( + x=np.arange(*params.duration_range), + phase=np.arange(*params.phase_range), + prob_target=1 - 2 * prob_target, + error_target=2 + * np.sqrt(prob_target * (1 - prob_target) / params.nshots), + prob_control=1 - 2 * prob_control, + error_control=2 + * np.sqrt(prob_control * (1 - prob_control) / params.nshots), + ), + ) + + return data + + +def _fit( + data: HamiltonianTomographyCANCPhaseData, +) -> HamiltonianTomographyCANCPhaseResults: + """Post-processing function for HamiltonianTomographyCANCPhase. + + We fit the expectation values using the Eq. S10 from the paper https://arxiv.org/pdf/2303.01427. + Afterwards, we extract the Hamiltonian terms from the fitted parameters. + + """ + hamiltonian_terms, fitted_parameters, pulses_phases, ham_tom_params, cr_lengths = ( + cancellation_phase_fit( + data=data, + ) + ) + + return HamiltonianTomographyCANCPhaseResults( + echo=data.echo, + hamiltonian_terms=hamiltonian_terms, + fitted_parameters=fitted_parameters, + cancellation_pulse_phases=pulses_phases, + hamiltonian_tom_params=ham_tom_params, + cr_lengths=cr_lengths, + verbose_plot=data.verbose_plot, + ) + + +def _plot( + data: HamiltonianTomographyCANCPhaseData, + target: QubitPairId, + fit: HamiltonianTomographyCANCPhaseResults | None = None, +): + """Plotting function for HamiltonianTomographyCANCPhase.""" + figs, fitting_report = cancellation_calibration_plot(data, target, fit) + + plot_ham_tom = True if fit is None else fit.verbose_plot + + if plot_ham_tom: + from .plotting import ( + tomography_cr_plot, + ) + + for phi in data.phases: + ampl_data = data.select_phase(phi) + if fit is not None: + selected_ham_terms = fit.select_pair_and_phase_ham_params(phi, target) + ham_tom_fit = HamiltonianTomographyResults( + echo=fit.echo, + hamiltonian_terms=selected_ham_terms, + fitted_parameters=fit.hamiltonian_tom_params[phi], + cr_lengths=fit.cr_lengths[phi], + ) + + fitting_report += "\n" + table_html( + table_dict( + 8 * [target], + [f"{term.name} [MHz]" for term in HamiltonianTerm] + + ["CR duration (ns)", "Cancellation phase [rad.]"], + [ + ham_tom_fit.hamiltonian_terms[target[0], target[1], term] + * kilo + for term in HamiltonianTerm + ] + + [ + fit.cr_lengths[phi][target] + if target in fit.cr_lengths[phi] + else None + ] + + [phi], + ) + ) + else: + ham_tom_fit = None + fitting_report = "" + + f, _ = tomography_cr_plot(ampl_data, target, ham_tom_fit) + figs += f + + return figs, fitting_report + + +def _update( + results: HamiltonianTomographyCANCPhaseResults, + platform: CalibrationPlatform, + target: QubitPairId, +): + # here is updated the full CNOT Pulse Sequence, which is composed by a CR sequence followe by a X_pi/2 and Z_(-pi/2) rotations on + # target and control qubit respectively + + target = target[::-1] if target not in results.cancellation_pulse_phases else target + + # now no check is needed since the acquisition was executed correctly, + # which means we have all parameters defined. + cr_pulse, canc_pulse = cross_resonance_pulses(platform, target[0], target[1]) + gate_duration = cr_pulse.duration + control_amplitude = cr_pulse.amplitude + target_amplitude = 0.0 if canc_pulse is None else canc_pulse.amplitude + + # check if the resulting fit was succsessfull + if target in results.cancellation_pulse_phases: + update_cnot_from_fit( + platform=platform, + pair=target, + cr_duration=gate_duration, + cr_ampl=control_amplitude, + control_phase=results.cancellation_pulse_phases[target]["control"], + canc_ampl=target_amplitude, + canc_phase=results.cancellation_pulse_phases[target]["target"], + echo_flag=results.echo, + ) + + +cancellation_phase_tuning = Routine(_acquisition, _fit, _plot, _update) +"""HamiltonianTomographyCANCPhase Routine object.""" + +"""See http://login.qrccluster.com:9000/bf6RezP1SpCnI861v6UNpA== for an example run on the emulator. +""" diff --git a/src/qibocal/protocols/two_qubit_interaction/cross_resonance/cr_parent_classes.py b/src/qibocal/protocols/two_qubit_interaction/cross_resonance/cr_parent_classes.py new file mode 100644 index 0000000000..afaaf73289 --- /dev/null +++ b/src/qibocal/protocols/two_qubit_interaction/cross_resonance/cr_parent_classes.py @@ -0,0 +1,111 @@ +from dataclasses import dataclass, field +from enum import Enum +from itertools import combinations +from typing import Any + +import numpy as np + +from qibocal.auto.operation import ( + Data, + Parameters, + QubitId, + QubitPairId, + Results, +) + +HamiltonianTomographyType = np.dtype( + [ + ("prob_target", np.float64), + ("error_target", np.float64), + ("prob_control", np.float64), + ("error_control", np.float64), + ("x", np.float64), + ] +) +"""Custom dtype for Hamiltonian Tomography.""" + + +class OverlappingQubitPairError(Exception): + """Raised when the target qubit pairs are not independent.""" + + pass + + +def check_qubit_overlap(targets: list[QubitPairId]) -> None: + """This function checks if the input qubit pairs are independent, i.e. they don't share any qubit.""" + if any(set(t1) & set(t2) for t1, t2 in combinations(targets, 2)): + raise OverlappingQubitPairError( + "Target pairs must be independent, but are overlapping." + ) + + +class SetControl(str, Enum): + """Helper to create sequence with control set to X or I.""" + + Id = "Id" + X = "X" + + +class Basis(str, Enum): + """Measurement basis.""" + + X = "X" + Y = "Y" + Z = "Z" + + +class HamiltonianTerm(str, Enum): + """Hamiltonian terms for CR effective Hamiltonian.""" + + IX = "IX" + IY = "IY" + IZ = "IZ" + ZX = "ZX" + ZY = "ZY" + ZZ = "ZZ" + + +@dataclass +class HamiltonianTomographyParameters(Parameters): + """Parent class for parameters in all time-sweeping hamiltonian tomography experiments.""" + + duration_range: tuple[float, float, float] + """Range of amplitudes for CR pulse (start, end, step).""" + echo: bool = False + """Apply echo sequence or not. + + The ECR is described in https://arxiv.org/pdf/1210.7011 + """ + interpolated_sweeper: bool = False + """Use real-time interpolation if supported by instruments.""" + + +@dataclass +class HamiltonianTomographyResults(Results): + """Results for Hamiltonian Tomography CR Length experiment.""" + + echo: bool + cr_lengths: dict[tuple[QubitId, QubitId], float] = field(default_factory=dict) + """Estimated durations of CR gate.""" + hamiltonian_terms: dict = field(default_factory=dict) + """Terms in effective Hamiltonian.""" + fitted_parameters: dict = field(default_factory=dict) + """Fitted parameters from X,Y,Z expectation values.""" + + def __contains__(self, pair: QubitPairId) -> bool: + return all(key[:2] == pair for key in list(self.fitted_parameters)) + + +@dataclass +class HamiltonianTomographyData(Data): + """Data for Hamiltonian Tomography CR Amplitude experiment.""" + + echo: bool + data: dict[tuple[QubitId, QubitId, Basis, SetControl], Any] = field( + default_factory=dict + ) + """Raw data acquired.""" + + @property + def pairs(self): + return {(i[0], i[1]) for i in self.data} diff --git a/src/qibocal/protocols/two_qubit_interaction/cross_resonance/cross_resonance_amplitude.py b/src/qibocal/protocols/two_qubit_interaction/cross_resonance/cross_resonance_amplitude.py new file mode 100644 index 0000000000..762ebe2f21 --- /dev/null +++ b/src/qibocal/protocols/two_qubit_interaction/cross_resonance/cross_resonance_amplitude.py @@ -0,0 +1,301 @@ +from dataclasses import dataclass, field + +import numpy as np +from qibolab import ( + AcquisitionType, + AveragingMode, + ParallelSweepers, + Parameter, + Sweeper, +) +from scipy.constants import kilo + +from qibocal.auto.operation import ( + Parameters, + QubitId, + QubitPairId, + Routine, +) +from qibocal.calibration import CalibrationPlatform +from qibocal.protocols.utils import table_dict, table_html + +from .cr_parent_classes import ( + Basis, + HamiltonianTerm, + HamiltonianTomographyData, + HamiltonianTomographyResults, + HamiltonianTomographyType, + SetControl, + check_qubit_overlap, +) +from .cross_resonance_processing import ( + extract_hamiltonian_terms, + tomography_cr_fit, +) +from .plotting import tomography_cr_plot +from .utils import ( + cross_resonance_experiment, + update_cnot_from_fit, +) + +__all__ = ["cr_amplitude"] + + +@dataclass(kw_only=True) +class HamiltonianTomographyCRAmplParameters(Parameters): + """HamiltonianTomographyCRAmplitude runcard inputs.""" + + pulse_duration: float + """Duration of CR pulse [ns].""" + amplitude_range: tuple[float, float, float] + """Range of amplitudes for CR pulse (start, end, step).""" + control_phase: float = 0.0 + """Phase of the CR pulse on the control qubit""" + target_amplitude: float | None = None + """Amplitude of the Cancellation pulse on the target qubit""" + target_phase: float = 0.0 + """Phase of the Cancellation pulse on the target qubit""" + echo: bool = False + """Apply echo sequence or not. + + The ECR is described in https://arxiv.org/pdf/1210.7011 + """ + + +@dataclass(kw_only=True) +class HamiltonianTomographyCRAmplResults(HamiltonianTomographyResults): + """HamiltonianTomographyCRAmpl outputs.""" + + cr_duration: float + control_phase: float + target_amplitude: float | None + target_phase: float + cr_amplitudes: dict[tuple[QubitId, QubitId], float] = field(default_factory=dict) + """Estimated amplitudes of CR gate.""" + hamiltonian_terms: dict[tuple[QubitId, QubitId, HamiltonianTerm], float] = field( + default_factory=dict + ) + """Terms in effective Hamiltonian.""" + fitted_parameters: dict[tuple[QubitId, QubitId, SetControl], list] = field( + default_factory=dict + ) + """Fitted parameters for Hamiltonian Terms values for different amplitudes.""" + + +@dataclass(kw_only=True) +class HamiltonianTomographyCRAmplData(HamiltonianTomographyData): + """Data structure for CR Amplitude.""" + + cr_duration: float + control_phase: float + target_amplitude: float | None + target_phase: float + + +def _acquisition( + params: HamiltonianTomographyCRAmplParameters, + platform: CalibrationPlatform, + targets: list[QubitPairId], +) -> HamiltonianTomographyCRAmplData: + """Data acquisition for Hamiltonian tomography CR protocol. + + We store the probability of the control qubit and the expectation value of the target qubit. + """ + + data = HamiltonianTomographyCRAmplData( + echo=params.echo, + cr_duration=params.pulse_duration, + control_phase=params.control_phase, + target_amplitude=params.target_amplitude, + target_phase=params.target_phase, + ) + + # check validity of input + check_qubit_overlap(targets) + + # update the CR channel with the target qubit frequency. + updates = [ + { + platform.qubits[c].drive_extra[t]: { + "frequency": platform.config(platform.qubits[t].drive).frequency + } + for c, t in targets + } + ] + + # We create one sequence for each combination of basis and control state, + # and we sweep the amplitude of the CR pulse for each of them; unrolling is not performed + for basis in Basis: + for setup in SetControl: + sequence, cr_pulses, _, _, _ = cross_resonance_experiment( + platform=platform, + pair_list=targets, + duration=params.pulse_duration, + ctrl_ampl=0.0, # this is the swept param + ctrl_phase=params.control_phase, + targ_ampl=params.target_amplitude, + targ_phase=params.target_phase, + basis=basis, + setup=setup, + echo=params.echo, + ) + + ampl_parsweepers = ParallelSweepers( + [ + Sweeper( + parameter=Parameter.amplitude, + range=params.amplitude_range, + pulses=[cr_pulses[pair][0] for pair in targets], + ) + ] + ) + if params.echo: + # sweeping over the out-of-phase signal (refocusing) + echo_ampl_range = tuple(-x for x in params.amplitude_range) + ampl_parsweepers += ParallelSweepers( + [ + Sweeper( + parameter=Parameter.amplitude, + range=echo_ampl_range, + pulses=[cr_pulses[pair][1] for pair in targets], + ) + ] + ) + + results = platform.execute( + [sequence], + [ampl_parsweepers], + nshots=params.nshots, + relaxation_time=params.relaxation_time, + acquisition_type=AcquisitionType.DISCRIMINATION, + averaging_mode=AveragingMode.CYCLIC, + updates=updates, + ) + for ctrl, trg in targets: + target_acq_handle = list( + sequence.channel(platform.qubits[trg].acquisition) + )[-1].id + control_acq_handle = list( + sequence.channel(platform.qubits[ctrl].acquisition) + )[-1].id + + prob_target = results[target_acq_handle].ravel() + prob_control = results[control_acq_handle].ravel() + + data.register_qubit( + HamiltonianTomographyType, + (ctrl, trg, basis, setup), + dict( + x=np.arange(*params.amplitude_range), + prob_target=1 - 2 * prob_target, + error_target=( + 2 * np.sqrt(prob_target * (1 - prob_target) / params.nshots) + ).tolist(), + prob_control=1 - 2 * prob_control, + error_control=( + 2 + * np.sqrt(prob_control * (1 - prob_control) / params.nshots) + ).tolist(), + ), + ) + return data + + +def _fit( + data: HamiltonianTomographyCRAmplData, +) -> HamiltonianTomographyCRAmplResults: + """Post-processing function for HamiltonianTomographyCRAmpl. + + We fit the expectation values using the Eq. S10 from the paper https://arxiv.org/pdf/2303.01427. + Afterwards, we extract the Hamiltonian terms from the fitted parameters. + + """ + fitted_parameters, cr_gate_ampls = tomography_cr_fit(data=data) + hamiltonian_terms = {} + for pair in data.pairs: + hamiltonian_terms |= extract_hamiltonian_terms( + pair=pair, fitted_parameters=fitted_parameters + ) + + return HamiltonianTomographyCRAmplResults( + echo=data.echo, + hamiltonian_terms=hamiltonian_terms, + fitted_parameters=fitted_parameters, + cr_amplitudes=cr_gate_ampls, + cr_duration=data.cr_duration, + control_phase=data.control_phase, + target_amplitude=data.target_amplitude, + target_phase=data.target_phase, + ) + + +def _plot( + data: HamiltonianTomographyCRAmplData, + target: QubitPairId, + fit: HamiltonianTomographyCRAmplResults | None = None, +): + """Plotting function for HamiltonianTomographyCRAmpl.""" + figs, fitting_report = tomography_cr_plot(data, target, fit) + + if fit is not None: + fitting_report = table_html( + table_dict( + 11 * [target], + ( + [f"{term.name} [MHz]" for term in HamiltonianTerm] + + [ + "CR duration (ns)", + "Control amplitude (a.u.)", + "Control phase (rad)", + "Target amplitude (a.u.)", + "Target phase (rad)", + ] + ), + ( + [ + fit.hamiltonian_terms[target[0], target[1], term] * kilo + for term in HamiltonianTerm + ] + + [fit.cr_duration] + + [ + fit.cr_amplitudes[target] + if target in fit.cr_amplitudes + else None + ] + + [fit.control_phase, fit.target_amplitude, fit.target_phase] + ), + ) + ) + else: + fitting_report = "" + return figs, fitting_report + + +def _update( + results: HamiltonianTomographyCRAmplResults, + platform: CalibrationPlatform, + target: QubitPairId, +): + # check if the resulting fit was succsessfull + if target in results.cr_amplitudes: + update_cnot_from_fit( + platform=platform, + pair=target, + cr_duration=results.cr_duration, + cr_ampl=results.cr_amplitudes[target], + control_phase=results.control_phase, + canc_ampl=results.target_amplitude, + canc_phase=results.target_phase, + echo_flag=results.echo, + ) + + +cr_amplitude = Routine(_acquisition, _fit, _plot, _update, two_qubit_gates=True) +"""Hamiltonian tomography protocol for CR gate calibration. + +This protocol computes the expectation values for X, Y and Z for the target qubit +after the application of a cross resonance sequence. The CR pulses are played on the control drive +channel with frequency set to the frequency of the target drive channel. +We repeat the measurement twice for each target qubit, once with the control qubit in state 0 +and once with the control qubit in state 1. +""" diff --git a/src/qibocal/protocols/two_qubit_interaction/cross_resonance/cross_resonance_length.py b/src/qibocal/protocols/two_qubit_interaction/cross_resonance/cross_resonance_length.py new file mode 100644 index 0000000000..80d5c7383c --- /dev/null +++ b/src/qibocal/protocols/two_qubit_interaction/cross_resonance/cross_resonance_length.py @@ -0,0 +1,327 @@ +"""Hamiltonian tomography protocol for CR gate calibration. + +This protocol computes the expectation values for X, Y and Z for the target qubit after the application of a cross resonance sequence. +The CR pulses are played on the control drive channel with frequency set to the frequency of the target drive channel. +""" + +from dataclasses import dataclass, field + +import numpy as np +from qibolab import ( + AcquisitionType, + AveragingMode, + ParallelSweepers, + Parameter, + Sweeper, +) +from scipy.constants import kilo + +from qibocal.auto.operation import ( + QubitPairId, + Routine, +) +from qibocal.calibration import CalibrationPlatform +from qibocal.protocols.utils import table_dict, table_html + +from .cr_parent_classes import ( + Basis, + HamiltonianTerm, + HamiltonianTomographyData, + HamiltonianTomographyParameters, + HamiltonianTomographyResults, + HamiltonianTomographyType, + SetControl, + check_qubit_overlap, +) +from .cross_resonance_processing import ( + extract_hamiltonian_terms, + tomography_cr_fit, +) +from .plotting import tomography_cr_plot +from .utils import ( + QubitId, + cross_resonance_experiment, + ro_delay_range, + update_cnot_from_fit, +) + +__all__ = ["cr_length"] + + +@dataclass(kw_only=True) +class HamiltonianTomographyCRLengthParameters(HamiltonianTomographyParameters): + """HamiltonianTomographyCRLength runcard inputs.""" + + pulse_amplitude: float + """CR pulse amplitude""" + phase: float = 0.0 + """Phase of CR pulse.""" + target_amplitude: float | None = None + """Amplitude of cancellation pulse.""" + target_phase: float = 0 + """Phase of target pulse.""" + + +@dataclass(kw_only=True) +class HamiltonianTomographyCRLengthResults(HamiltonianTomographyResults): + """HamiltonianTomographyCRLength outputs.""" + + control_amplitude: float + control_phase: float = 0 + target_amplitude: float | None = None + target_phase: float = 0 + hamiltonian_terms: dict[tuple[QubitId, QubitId, HamiltonianTerm], float] = field( + default_factory=dict + ) + """Terms in effective Hamiltonian.""" + fitted_parameters: dict[tuple[QubitId, QubitId, SetControl], list] = field( + default_factory=dict + ) + """Fitted parameters from X,Y,Z expectation values.""" + + +@dataclass(kw_only=True) +class HamiltonianTomographyCRLengthData(HamiltonianTomographyData): + """Data structure for CR length.""" + + control_amplitude: float + control_phase: float + target_amplitude: float | None + target_phase: float + + +def _acquisition( + params: HamiltonianTomographyCRLengthParameters, + platform: CalibrationPlatform, + targets: list[QubitPairId], +) -> HamiltonianTomographyCRLengthData: + """Data acquisition for Hamiltonian tomography CR protocol.""" + + data = HamiltonianTomographyCRLengthData( + echo=params.echo, + control_amplitude=params.pulse_amplitude, + control_phase=params.phase, + target_amplitude=params.target_amplitude, + target_phase=params.target_phase, + ) + + # check validity of input + check_qubit_overlap(targets) + + # update the CR channel with the target qubit frequency. + updates = [ + { + platform.qubits[c].drive_extra[t]: { + "frequency": platform.config(platform.qubits[t].drive).frequency + } + for c, t in targets + } + ] + + # We create one sequence for each combination of basis and control state, + # and we sweep the duration of the CR pulse for each of them; unrolling is not performed + for basis in Basis: + for setup in SetControl: + sequence, cr_pulses, cr_target_pulses, cr_delays, ro_delays = ( + cross_resonance_experiment( + platform=platform, + pair_list=targets, + duration=0.0, # this is the swept param + ctrl_ampl=params.pulse_amplitude, + ctrl_phase=params.phase, + targ_ampl=params.target_amplitude, + targ_phase=params.target_phase, + basis=basis, + setup=setup, + echo=params.echo, + interpolated_sweeper=params.interpolated_sweeper, + ) + ) + + if params.interpolated_sweeper: + duration_parsweepers = ParallelSweepers( + [ + Sweeper( + parameter=Parameter.duration_interpolated, + range=params.duration_range, + pulses=cr_pulses[pair] + cr_target_pulses[pair], + ) + for pair in targets + ] + ) + + else: + duration_parsweepers = ParallelSweepers( + [ + Sweeper( + parameter=Parameter.duration, + range=params.duration_range, + pulses=cr_pulses[pair] + + cr_target_pulses[pair] + + cr_delays[pair], + ) + for pair in targets + ] + ) + ro_parsweepers = ParallelSweepers( + [ + Sweeper( + parameter=Parameter.duration, + range=ro_delay_range( + cr_pulse_duration_range=params.duration_range, + echo=params.echo, + cntl_setup=setup, + control=c, + platform=platform, + ), + pulses=ro_delays[(c, t)], + ) + for c, t in targets + ] + ) + duration_parsweepers += ro_parsweepers + + results = platform.execute( + [sequence], + [duration_parsweepers], + nshots=params.nshots, + relaxation_time=params.relaxation_time, + acquisition_type=AcquisitionType.DISCRIMINATION, + averaging_mode=AveragingMode.CYCLIC, + updates=updates, + ) + + for control, target in targets: + target_acq_handle = list( + sequence.channel(platform.qubits[target].acquisition) + )[-1].id + control_acq_handle = list( + sequence.channel(platform.qubits[control].acquisition) + )[-1].id + + prob_target = results[target_acq_handle].ravel() + prob_control = results[control_acq_handle].ravel() + + data.register_qubit( + HamiltonianTomographyType, + (control, target, basis, setup), + dict( + x=np.arange(*params.duration_range), + prob_target=1 - 2 * prob_target, + error_target=( + 2 * np.sqrt(prob_target * (1 - prob_target) / params.nshots) + ).tolist(), + prob_control=1 - 2 * prob_control, + error_control=( + 2 + * np.sqrt(prob_control * (1 - prob_control) / params.nshots) + ).tolist(), + ), + ) + + return data + + +def _fit( + data: HamiltonianTomographyCRLengthData, +) -> HamiltonianTomographyCRLengthResults: + """Post-processing function for HamiltonianTomographyCRLength. + + We fit the expectation values using the Eq. S10 from the paper https://arxiv.org/pdf/2303.01427. + Afterwards, we extract the Hamiltonian terms from the fitted parameters. + + """ + fitted_parameters, cr_gate_lengths = tomography_cr_fit( + data=data, + ) + hamiltonian_terms = {} + for pair in data.pairs: + hamiltonian_terms |= extract_hamiltonian_terms( + pair=pair, fitted_parameters=fitted_parameters + ) + + return HamiltonianTomographyCRLengthResults( + echo=data.echo, + hamiltonian_terms=hamiltonian_terms, + fitted_parameters=fitted_parameters, + cr_lengths=cr_gate_lengths, + control_amplitude=data.control_amplitude, + control_phase=data.control_phase, + target_amplitude=data.target_amplitude, + target_phase=data.target_phase, + ) + + +def _plot( + data: HamiltonianTomographyCRLengthData, + target: QubitPairId, + fit: HamiltonianTomographyCRLengthResults | None = None, +): + """Plotting function for HamiltonianTomographyCRLength.""" + figs, fitting_report = tomography_cr_plot(data, target, fit) + + if fit is not None: + fitting_report = table_html( + table_dict( + 11 * [target], + ( + [f"{term.name} [MHz]" for term in HamiltonianTerm] + + [ + "CR duration (ns)", + "Control amplitude (a.u.)", + "Control phase (rad)", + "Target amplitude (a.u.)", + "Target phase (rad)", + ] + ), + ( + [ + fit.hamiltonian_terms[target[0], target[1], term] * kilo + for term in HamiltonianTerm + ] + + [fit.cr_lengths[target] if target in fit.cr_lengths else None] + + [ + fit.control_amplitude, + fit.control_phase, + fit.target_amplitude, + fit.target_phase, + ] + ), + ) + ) + else: + fitting_report = "" + return figs, fitting_report + + +def _update( + results: HamiltonianTomographyCRLengthResults, + platform: CalibrationPlatform, + target: QubitPairId, +): + + # check if the resulting fit was succsessfull + if target in results.cr_lengths: + update_cnot_from_fit( + platform=platform, + pair=target, + cr_duration=results.cr_lengths[target], + cr_ampl=results.control_amplitude, + control_phase=results.control_phase, + canc_ampl=results.target_amplitude, + canc_phase=results.target_phase, + echo_flag=results.echo, + ) + + +cr_length = Routine(_acquisition, _fit, _plot, _update, two_qubit_gates=True) +"""HamiltonianTomographyCRLength Routine object. + +We measure the expectation values X,Y and Z on the target qubit after +applying the CR sequence specified by the input parameters. We repeat the +measurement twice for each target qubit, once with the control qubit in state 0 +and once with the control qubit in state 1. + +Check http://login.qrccluster.com:9000/u7Xw0C4_Rti6s2JdJ_HI4g== for a good experiment result found on emulator with 1% of classical crosstalk between +the two qubits. +""" diff --git a/src/qibocal/protocols/two_qubit_interaction/cross_resonance/cross_resonance_processing.py b/src/qibocal/protocols/two_qubit_interaction/cross_resonance/cross_resonance_processing.py new file mode 100644 index 0000000000..ffcd1114fa --- /dev/null +++ b/src/qibocal/protocols/two_qubit_interaction/cross_resonance/cross_resonance_processing.py @@ -0,0 +1,723 @@ +from collections.abc import Callable +from typing import Union + +import numpy as np +from numpy.typing import NDArray +from scipy.optimize import ( + Bounds, + NonlinearConstraint, + curve_fit, + differential_evolution, + minimize, +) + +from qibocal.auto.operation import ( + QubitId, + QubitPairId, +) +from qibocal.config import log +from qibocal.protocols.utils import angle_wrap, quinn_fernandes_algorithm + +from . import fitting +from .cr_parent_classes import ( + Basis, + HamiltonianTerm, + HamiltonianTomographyData, + SetControl, +) + +QUANTILE_CONSTANT = 1.5 + + +def dynamic_evolution_optimizer( + signals_id: NDArray, + x: NDArray, + init_omega_guess: float, + use_constraints: bool = False, +) -> NDArray: + """Optimizer for sinusoidal fitting; it exploits evolution algorithm to find the best fit parameters. + This algorithm is a gradient-free optimization, hence might be less affected by local minima in the cost + function landscape, but it can be computationally expensive, especially for large datasets or complex models. + """ + + assert x.ndim == 1, f"Expected 1D array, got array with shape {x.shape}" + + def func_to_minimize(z): + return np.sum( + (fitting.simultaneous_expectation(x, *z) - signals_id.ravel()) ** 2 + ) + + def constraint(x): + return np.sqrt(np.sum(x**2)) + + bounds = Bounds( + [-init_omega_guess, -init_omega_guess, -init_omega_guess, 0], + [init_omega_guess, init_omega_guess, init_omega_guess, np.inf], + ) + res = differential_evolution( + func_to_minimize, + bounds, + maxiter=int(1e6), + constraints=( + NonlinearConstraint(constraint, 0, init_omega_guess * 2) + if use_constraints + else () + ), + ) + popt = res.x + return popt + + +def scipy_curve_fit_optimizer( + concatenated_signal: NDArray, + vector_x: NDArray, + init_omega_guess: float, +) -> NDArray: + """Optimizer for sinusoidal fitting; it exploits gradient-based algorithms to find the best fit parameters. + This algorithm is a gradient-base optimization, hence might be affected by local minima in the cost + function landscape when dealing with sinusoidal functions due to their periodicity. + """ + + sampling_rate = 1 / abs(vector_x[1] - vector_x[0]) + period = int(2 * np.pi / init_omega_guess * sampling_rate) + + offsets = np.median(concatenated_signal[:, :period], axis=-1) * init_omega_guess**2 + + omega_guesses = np.zeros(offsets.shape) + omega_guesses[-1] = np.sqrt(np.abs(offsets[-1])) + omega_guesses[1] = offsets[1] / omega_guesses[-1] + omega_guesses[0] = offsets[0] / omega_guesses[-1] + omega_guesses = np.concatenate([omega_guesses, [0]]) + + # even though vector_x does not have the same shape as concatenated_signal, + # curve_fit will still work (despite it's documented the opposite) because it's just + # necesary fitting.simultaneous_expectation's output to have the same shape as concatenated_signal. + popt, _ = curve_fit( + fitting.simultaneous_expectation, + vector_x, + concatenated_signal.ravel(), + maxfev=int(1e6), + p0=omega_guesses, + bounds=([-np.inf, -np.inf, -np.inf, 0], [np.inf, np.inf, np.inf, np.inf]), + ) + + return popt + + +def numerical_root_finder( + root_func: Callable[..., float], + x_range: NDArray | list[float], + tol: float, + **kwargs, +): + """Finds the root of a generic function :data:`root_func` by minimizing its absolute value. + + The solver first performs a coarse grid search across the provided range to + identify a starting candidate, then refines the result using the + Nelder-Mead simplex algorithm; then, if the search is successful, the + algorithm returns the value where the function is closest to zero. + This function minimizes `abs(root_func(x))`. If the function does not + cross zero, it will return the local minimum of the absolute value. + """ + + def func_to_solve(x): + return np.abs(root_func(x, **kwargs)) + + x_grid = np.linspace(x_range[0], x_range[-1], 100 * len(x_range)) + + y_vals = np.array([func_to_solve(xi) for xi in x_grid]).flatten() + + cr_sig = x_grid[y_vals - np.min(y_vals) <= tol][0] + + # Use a numerical minimizer starting from our guess + res = minimize( + func_to_solve, + x0=cr_sig, + method="Nelder-Mead", + ) + + if res.success: + cr_sig = min(res.x) + + return cr_sig if cr_sig <= np.max(x_grid) else np.max(x_grid) + + +def compute_total_expectation_value( + data: HamiltonianTomographyData, + pair: QubitPairId, +) -> tuple[NDArray, NDArray]: + """Given a Qubit Pair :data:`pair`=(control, target), it computes the expectation + values for each Pauli Basis for the target and control qubits respectively. + """ + + tot_exp_vals_target = [] + tot_exp_vals_control = [] + for basis in Basis: + tot_exp_vals_target.append( + data.data[pair[0], pair[1], basis, SetControl.Id].prob_target + + data.data[pair[0], pair[1], basis, SetControl.X].prob_target + ) + tot_exp_vals_control.append( + data.data[pair[0], pair[1], basis, SetControl.Id].prob_control + + data.data[pair[0], pair[1], basis, SetControl.X].prob_control + ) + + return np.vstack(tot_exp_vals_target), np.vstack(tot_exp_vals_control) + + +def bloch_func( + x: list[float] | NDArray, pair: QubitPairId, fitted_parameters: dict +) -> NDArray: + """Given the fitted parameters for the target's Pauli expectation values + for either control qubit in |0> or in |1>, computes the estimated Bloch vector. + """ + + id_blochfit = fitting.simultaneous_expectation( + x, *fitted_parameters[pair[0], pair[1], SetControl.Id] + ).reshape((3, -1)) + x_blochfit = fitting.simultaneous_expectation( + x, *fitted_parameters[pair[0], pair[1], SetControl.X] + ).reshape((3, -1)) + return np.sqrt(np.sum((id_blochfit + x_blochfit) ** 2, axis=0)) / 2 + + +def compute_bloch_vector( + data: Union[ + "HamiltonianTomographyCRLengthData", # noqa: F821 + "HamiltonianTomographyCRAmplData", # noqa: F821 + ], + pair: QubitPairId, + fitted_parameters: dict | None = None, +) -> tuple[NDArray, NDArray | None, NDArray]: + """For a given qubit pair :data:`pair`, it computes the Bloch vector R for each data point and + also estimates the Bloch vector using the fitted parameters of the Hamiltonian Tomography. + See `arXiv:1603.04821 `__ for further information. + """ + + bloch_exp_targ, bloch_exp_ctrl = compute_total_expectation_value(data, pair) + bloch_exp_targ = np.sqrt(np.sum((bloch_exp_targ) ** 2, axis=0)) / 2 + bloch_exp_ctrl = np.sqrt(np.sum((bloch_exp_ctrl) ** 2, axis=0)) / 2 + + bloch_fit_targ = None + if fitted_parameters is not None: + times = data.data[pair[0], pair[1], Basis.Z, SetControl.Id].x + times_range = np.linspace(min(times), max(times), 2 * len(times)) + + bloch_fit_targ = bloch_func(times_range, pair, fitted_parameters) + + return bloch_exp_targ, bloch_fit_targ, bloch_exp_ctrl + + +def estimate_cr_param( + x_range: NDArray, + data: HamiltonianTomographyData, + pair: QubitPairId, + fitted_parameters: dict, + tol: float = 1e-6, +) -> float: + """Function for estimating important parameters for the cross resonance, depending on the + specific experiment run; if :data:`data` is type :class:`HamiltonianTomographyData` it finds the + thuned pulse duration, while if :data:`data` is type :class:`HamiltonianTomographyCRAmplData` it + finds the tuned pulse amplitude. + + The cross resonance parameter is computed by finding the value that solves the Bloch vector R. + To do so, the function checks whether all the fits in the Hamiltonian Tomography experiment succeeded; + if so the Bloch vector is computed by using the fitted parameter, otherwise it simply computes + R only for the acquired datapoints. + """ + + if all([(pair[0], pair[1], s) in fitted_parameters for s in SetControl]): + bloch_data, _, _ = compute_bloch_vector(data, pair, fitted_parameters) + x_range = data.data[pair[0], pair[1], Basis.Z, SetControl.Id].x + param = numerical_root_finder( + root_func=bloch_func, + x_range=x_range, + tol=tol, + pair=pair, + fitted_parameters=fitted_parameters, + ) + else: + bloch_data, _, _ = compute_bloch_vector(data, pair) + idx = np.argmin(bloch_data) + param = x_range[idx] + + return float(param) + + +def tune_cancellation_sequence( + x: list[float] | NDArray, + function_to_tune: Callable[..., float], + interactions_to_analyze: list[HamiltonianTerm], + ham_term: dict, + fit_params: dict, + tuned_keys: list[str], + tol: float, +) -> dict[str, float]: + """Function for estimating parameters for the cancellation pulse sequence, depending on the + specific experiment (either pulses phase tuning or cancellation amplitude tuning). + + The specific parameter is computed by finding the value that solves an input function :data:`function_to_tune`. + To do so, the function checks whether all the fits in the Hamiltonian Tomography experiment succeeded; + if so roots are solved by using the fitted parameter, otherwise it simply computes :data:`function_to_tune` only + for the acquired datapoints and the closest value to 0 is selected. + + These parameters are computed for every term of the Cross Resonance Hamiltonian: IX, ZX, IY, ZY, IZ, ZZ and saved + into a dictionary. + """ + + assert len(tuned_keys) == len(interactions_to_analyze), ( + """tuned_keys and interactions_to_analyze must be equally long.""" + ) + + # converting list into numpy array + x = np.array(x) + + tuned_parameters = {} + for ham_int, k in zip(interactions_to_analyze, tuned_keys): + if ham_int in fit_params: + tuned_parameters[k] = float( + numerical_root_finder( + root_func=function_to_tune, + x_range=x, + tol=tol, + **fit_params[ham_int], + ) + ) + else: + selected_ham_term = np.abs(np.array(ham_term[ham_int])) + min_idx = np.argmin(selected_ham_term) + tuned_parameters[k] = float(x[min_idx]) + + return tuned_parameters + + +def estimate_cancellation_amplitudes( + amplitudes: list[float] | NDArray, + ham_term: dict, + ampl_params: dict, + tol: float = 1e-8, +) -> dict[str, float]: + """Extrapolates the cancellation pulse amplitude for the cross resonance pulse sequence. + + This function estimates the optimal amplitudes for cancelling unwanted Hamiltonian terms + (IX and IY) in the cross resonance pulse sequence. It uses linear fitting to determine + the amplitude values that minimize these interaction terms. + """ + + interaction_terms = [HamiltonianTerm.IX, HamiltonianTerm.IY] + amps_names = ["ampl_ix", "ampl_iy"] + + tuned_amplitudes = tune_cancellation_sequence( + x=amplitudes, + function_to_tune=fitting.linear_func, + interactions_to_analyze=interaction_terms, + ham_term=ham_term, + fit_params=ampl_params, + tuned_keys=amps_names, + tol=tol, + ) + + return tuned_amplitudes + + +def estimate_cr_phases( + phases: list[float] | NDArray, + ham_term: dict, + phase_params: dict, + tol: float = 1e-6, +) -> tuple[float, float]: + """Extrapolates the cancellation pulse phases for the cross resonance pulse sequence. + + This function estimates the optimal phases for cancelling unwanted Hamiltonian terms + (ZY and IY) in the cross resonance pulse sequence. It uses sinusoidal fitting to determine + the phase values that minimize these interaction terms. + """ + + interaction_terms = [HamiltonianTerm.ZY, HamiltonianTerm.IY] + phases_names = ["phi0", "phi1"] + + tuned_phases = tune_cancellation_sequence( + x=phases, + function_to_tune=fitting.sin_func, + interactions_to_analyze=interaction_terms, + ham_term=ham_term, + fit_params=phase_params, + tuned_keys=phases_names, + tol=tol, + ) + + zx_value = fitting.sin_func( + tuned_phases["phi0"], + **phase_params[HamiltonianTerm.ZX], + ) + + zx_semiperiod_value = fitting.sin_func( + tuned_phases["phi0"] + np.pi, + **phase_params[HamiltonianTerm.ZX], + ) + + if zx_value > zx_semiperiod_value: + # in https://journals.aps.org/pra/pdf/10.1103/PhysRevA.93.060302 + # it is said we need to choose the CR phase that minimizes ZY components + # and maximizes ZX interaction. + tuned_phases["phi0"] += np.pi + + ix_value = fitting.sin_func( + tuned_phases["phi1"], + **phase_params[HamiltonianTerm.IX], + ) + + ix_semiperiod_value = fitting.sin_func( + tuned_phases["phi1"] + np.pi, + **phase_params[HamiltonianTerm.IX], + ) + + if ix_value > ix_semiperiod_value: + # same as above, but this time is more euristic. + tuned_phases["phi1"] += np.pi + + return angle_wrap(tuned_phases["phi0"]), angle_wrap( + tuned_phases["phi0"] - tuned_phases["phi1"] + ) + + +def tomography_cr_fit( + data: HamiltonianTomographyData, + fit_with_evolution: bool = False, +) -> tuple[ + dict[tuple[QubitId, QubitId, SetControl], list[float]], + dict[tuple[QubitId, QubitId], float], +]: + """Fit Hamiltonian tomography data for cross-resonance gates. + + This function performs sinusoidal fitting on tomography data collected for + cross-resonance interactions between qubit pairs. It fits the measurement + probabilities in X, Y, and Z bases across different control settings. + """ + + fitted_parameters = {} + cr_gate_x = {} + + for pair in data.pairs: + vector_x = data.data[pair[0], pair[1], Basis.X, SetControl.Id].x + for setup in SetControl: + concatenated_signal = np.concatenate( + [ + data.data[pair[0], pair[1], Basis.X, setup].prob_target, + data.data[pair[0], pair[1], Basis.Y, setup].prob_target, + data.data[pair[0], pair[1], Basis.Z, setup].prob_target, + ] + ).reshape(len(Basis), -1) + + total_omega_guess = quinn_fernandes_algorithm( + concatenated_signal, vector_x, speedup_flag=True + ) + + if fit_with_evolution: + popt = dynamic_evolution_optimizer( + concatenated_signal, + vector_x, + total_omega_guess, + ) + else: + popt = scipy_curve_fit_optimizer( + concatenated_signal, + vector_x, + total_omega_guess, + ) + + fitted_parameters[pair[0], pair[1], setup] = popt.tolist() + + cr_gate_x[pair[0], pair[1]] = estimate_cr_param( + vector_x, data, pair, fitted_parameters + ) + + return fitted_parameters, cr_gate_x + + +def extract_hamiltonian_terms( + pair: QubitPairId, + fitted_parameters: dict[tuple[QubitId, QubitId, SetControl], list[float]], +) -> dict[QubitPairId, dict[HamiltonianTerm, float]]: + """Extract Hamiltonian terms from fitted parameters. + + We follow the procedure presented in the paper https://arxiv.org/pdf/2303.01427. + """ + hamiltonian_terms = {} + hamiltonian_terms[pair[0], pair[1], HamiltonianTerm.ZX] = 0.5 * ( + fitted_parameters[pair[0], pair[1], SetControl.Id][0] + - fitted_parameters[pair[0], pair[1], SetControl.X][0] + ) + hamiltonian_terms[pair[0], pair[1], HamiltonianTerm.IX] = 0.5 * ( + fitted_parameters[pair[0], pair[1], SetControl.Id][0] + + fitted_parameters[pair[0], pair[1], SetControl.X][0] + ) + hamiltonian_terms[pair[0], pair[1], HamiltonianTerm.ZY] = 0.5 * ( + fitted_parameters[pair[0], pair[1], SetControl.Id][1] + - fitted_parameters[pair[0], pair[1], SetControl.X][1] + ) + hamiltonian_terms[pair[0], pair[1], HamiltonianTerm.IY] = 0.5 * ( + fitted_parameters[pair[0], pair[1], SetControl.Id][1] + + fitted_parameters[pair[0], pair[1], SetControl.X][1] + ) + hamiltonian_terms[pair[0], pair[1], HamiltonianTerm.ZZ] = 0.5 * ( + fitted_parameters[pair[0], pair[1], SetControl.Id][2] + - fitted_parameters[pair[0], pair[1], SetControl.X][2] + ) + hamiltonian_terms[pair[0], pair[1], HamiltonianTerm.IZ] = 0.5 * ( + fitted_parameters[pair[0], pair[1], SetControl.Id][2] + + fitted_parameters[pair[0], pair[1], SetControl.X][2] + ) + return hamiltonian_terms + + +def refactor_hamiltonian_terms( + ham_terms: dict[tuple[QubitId, QubitId, HamiltonianTerm], float], + pair: QubitPairId, +) -> dict[HamiltonianTerm, float]: + """Refactor Hamiltonian terms by removing qubit pair information from keys. + + Converts dictionary keys from (qubit_id_0, qubit_id_1, HamiltonianTerm) tuples + to just HamiltonianTerm, simplifying the data structure. + """ + + for term in HamiltonianTerm: + ham_terms[term] = ham_terms.pop((pair[0], pair[1], term)) + + return ham_terms + + +def reconstruct_full_hamiltonian_terms( + ham_terms: dict[HamiltonianTerm, float], + pair: QubitPairId, +) -> dict[tuple[QubitId, QubitId, HamiltonianTerm], float]: + """ + Reconstruct full Hamiltonian terms dictionary by adding qubit pair information to keys. + + Converts dictionary keys from HamiltonianTerm to (qubit_id_0, qubit_id_1, HamiltonianTerm) tuples, + restoring the original data structure used for Hamiltonian terms. + """ + for term in HamiltonianTerm: + ham_terms[(pair[0], pair[1], term)] = ham_terms.pop(term) + + return ham_terms + + +def amp_tom_fit( + x: list[float] | NDArray, + y: list[float] | NDArray, + q_pair: QubitPairId, + term: HamiltonianTerm, + result_dict: dict[HamiltonianTerm, dict[str, float]], +) -> dict[HamiltonianTerm, dict[str, float]]: + """Fit linear function to amplitude vs Hamiltonian term data. + + Performs a linear fit on the provided data to extract amplitude-dependent + parameters for a specific Hamiltonian term. + """ + try: + pguess = [0, 0] + popt, _ = curve_fit( + fitting.linear_func, + x, + y, + p0=pguess, + maxfev=int(1e6), + absolute_sigma=True, + bounds=([-np.inf, -np.inf], [np.inf, np.inf]), + ) + result_dict[term] = {"a": popt[0], "b": popt[1]} + + except Exception as e: + log.warning(f"{term} term vs amplitudes fit failed for {q_pair} due to {e}.") + + return result_dict + + +def cancellation_amplitude_fit( + data: HamiltonianTomographyData, +) -> tuple[ + dict[QubitPairId, list[tuple[float, dict[HamiltonianTerm, float]]]], + dict[QubitPairId, dict[HamiltonianTerm, dict[str, float]]], + dict[QubitPairId, dict[str, float]], + dict[float, dict[tuple[QubitId, QubitId, SetControl], list[float]]], + dict[float, dict[tuple[QubitId, QubitId], float]], +]: + """Perform amplitude-dependent tomography fitting for calibrating + cross resonance cancellation pulse. + + Fits the dependence of Hamiltonian term parameters on the CR pulse amplitude. + Extracts Hamiltonian terms at different amplitudes and fits their variation + with amplitude to obtain linear parameters and cancellation amplitudes. + """ + + amp_hamiltonian_params: dict[ + QubitPairId, list[tuple[float, dict[HamiltonianTerm, float]]] + ] = {} + ham_tomography_dict: dict[ + float, dict[tuple[QubitId, QubitId, SetControl], list[float]] + ] = {} + gate_duration_dict: dict[float, dict[tuple[QubitId, QubitId], float]] = {} + for amp in data.amplitudes: + amp_data = data.select_amplitude(amp) + length_params, cr_durations = tomography_cr_fit(amp_data) + gate_duration_dict[amp] = cr_durations + ham_tomography_dict[amp] = length_params + for pair in amp_data.pairs: + terms = extract_hamiltonian_terms(pair, length_params) + terms = refactor_hamiltonian_terms(terms, pair) + res_tuple = (amp, terms) + if pair not in amp_hamiltonian_params: + amp_hamiltonian_params[pair] = [res_tuple] + else: + amp_hamiltonian_params[pair].append(res_tuple) + + amp_lin_fit_params = {} + calibrated_amplitudes = {} + for pair_key, pair_value in amp_hamiltonian_params.items(): + num_terms = {t: [] for t in HamiltonianTerm} + amplitudes = [] + for vals in pair_value: + amplitudes.append(vals[0]) + for t in HamiltonianTerm: + num_terms[t].append(vals[1][t]) + + fit_params_pair = {} + for t in HamiltonianTerm: + fit_params_pair = amp_tom_fit( + x=amplitudes, + y=num_terms[t], + q_pair=pair_key, + term=t, + result_dict=fit_params_pair, + ) + + amp_lin_fit_params[pair_key] = fit_params_pair + target_amplitudes = estimate_cancellation_amplitudes( + amplitudes=amplitudes, ham_term=num_terms, ampl_params=fit_params_pair + ) + calibrated_amplitudes[pair_key] = target_amplitudes + + return ( + amp_hamiltonian_params, + amp_lin_fit_params, + calibrated_amplitudes, + ham_tomography_dict, + gate_duration_dict, + ) + + +def phase_tom_fit( + x: list[float] | NDArray, + y: list[float] | NDArray, + q_pair: QubitPairId, + term: HamiltonianTerm, + result_dict: dict[HamiltonianTerm, dict[str, float]], +) -> dict[HamiltonianTerm, dict[str, float]]: + """Fit sinusoidal function to phase vs Hamiltonian term data. + + Performs a sinusoidal fit on the provided data to extract phase-dependent + parameters for a specific Hamiltonian term. + """ + + median_sig = np.median(y) + q80 = np.quantile(y, 0.8) + q20 = np.quantile(y, 0.2) + amplitude_guess = abs(q80 - q20) / QUANTILE_CONSTANT + phase_guess = 0 + pguess = [amplitude_guess, median_sig, phase_guess] + try: + popt, _ = curve_fit( + lambda x, a, b, phi: fitting.sin_func(x, a, b, 1, phi), + x, + y, + p0=pguess, + maxfev=int(1e6), + absolute_sigma=True, + bounds=( + [-np.inf, -np.inf, -np.inf], + [np.inf, np.inf, np.inf], + ), + ) + result_dict[term] = { + "a": popt[0], + "b": popt[1], + "omega": 1, + "phi": popt[2], + } + except Exception as e: + log.warning(f"{term} term vs amplitudes fit failed for {q_pair} due to {e}.") + + return result_dict + + +def cancellation_phase_fit( + data: HamiltonianTomographyData, +) -> tuple[ + dict[QubitPairId, list[tuple[float, dict[HamiltonianTerm, float]]]], + dict[QubitPairId, dict[HamiltonianTerm, dict[str, float]]], + dict[QubitPairId, dict[str, float]], + dict[float, dict[tuple[QubitId, QubitId, SetControl], list[float]]], + dict[float, dict[tuple[QubitId, QubitId], float]], +]: + """Fit phase-dependent Hamiltonian parameters using cross-resonance tomography. + + Extracts and fits Hamiltonian terms for different phases across all qubit pairs, + performing sinusoidal fits on phase-dependent data and estimating cancelling phases + for control and target qubits. + """ + + phase_hamiltonian_params: dict[ + QubitPairId, list[tuple[float, dict[HamiltonianTerm, float]]] + ] = {} + ham_tomography_dict: dict[ + float, dict[tuple[QubitId, QubitId, SetControl], list[float]] + ] = {} + gate_duration_dict: dict[float, dict[tuple[QubitId, QubitId], float]] = {} + for p in data.phases: + phase_data = data.select_phase(p) + length_params, cr_durations = tomography_cr_fit(phase_data) + gate_duration_dict[p] = cr_durations + ham_tomography_dict[p] = length_params + for pair in phase_data.pairs: + terms = extract_hamiltonian_terms(pair, length_params) + terms = refactor_hamiltonian_terms(terms, pair) + res_tuple = (p, terms) + if pair not in phase_hamiltonian_params: + phase_hamiltonian_params[pair] = [res_tuple] + else: + phase_hamiltonian_params[pair].append(res_tuple) + + phase_sin_fit_params = {} + cancellating_phases = {} + for pair_key, pair_value in phase_hamiltonian_params.items(): + num_terms = {t: [] for t in HamiltonianTerm} + phases = [] + for vals in pair_value: + phases.append(vals[0]) + for t in HamiltonianTerm: + num_terms[t].append(vals[1][t]) + + fit_params_pair = {} + for t in HamiltonianTerm: + fit_params_pair = phase_tom_fit( + x=phases, + y=num_terms[t], + q_pair=pair_key, + term=t, + result_dict=fit_params_pair, + ) + + phase_sin_fit_params[pair_key] = fit_params_pair + ctrl_phase, trgt_phase = estimate_cr_phases( + phases=phases, ham_term=num_terms, phase_params=fit_params_pair + ) + cancellating_phases[pair_key] = {"control": ctrl_phase, "target": trgt_phase} + + return ( + phase_hamiltonian_params, + phase_sin_fit_params, + cancellating_phases, + ham_tomography_dict, + gate_duration_dict, + ) diff --git a/src/qibocal/protocols/two_qubit_interaction/cross_resonance/fitting.py b/src/qibocal/protocols/two_qubit_interaction/cross_resonance/fitting.py new file mode 100644 index 0000000000..1f589d3bd4 --- /dev/null +++ b/src/qibocal/protocols/two_qubit_interaction/cross_resonance/fitting.py @@ -0,0 +1,88 @@ +"""Fitting function for CR tomography.""" + +import numpy as np + + +def pauli_z_expectation( + t: np.typing.NDArray, wx: float, wy: float, wz: float, gamma: float +) -> np.typing.NDArray: + """Pauli Z expectation value for CR tomography. + + See https://arxiv.org/pdf/2303.01427 Eq. S10. + """ + return vectorized_simultaneous_expectation(t, wx, wy, wz, gamma).reshape(3, -1)[ + 2, : + ] + + +def pauli_x_expectation( + t: np.typing.NDArray, wx: float, wy: float, wz: float, gamma: float +) -> np.typing.NDArray: + """Pauli X expectation value for CR tomography. + + See https://arxiv.org/pdf/2303.01427 Eq. S10. + """ + return vectorized_simultaneous_expectation(t, wx, wy, wz, gamma)[0, :] + + +def pauli_y_expectation( + t: np.typing.NDArray, wx: float, wy: float, wz: float, gamma: float +) -> np.typing.NDArray: + """Pauli Y expectation value for CR tomography. + + See https://arxiv.org/pdf/2303.01427 Eq. S10. + """ + return vectorized_simultaneous_expectation(t, wx, wy, wz, gamma)[1, :] + + +def vectorized_simultaneous_expectation( + t: np.typing.NDArray, + wx: float, + wy: float, + wz: float, + gamma: float, +) -> np.typing.NDArray: + + # Compute the norm of the w vector with components wx, wy, wz + w = np.sqrt(wx**2 + wy**2 + wz**2) + + # Map cosine, sine, and constant time components to X, Y, Z expectation values. + system_matrix = np.array( + [ + [-wx * wz, w * wy, wx * wz], + [-wy * wz, -w * wx, wy * wz], + [wx**2 + wy**2, 0, wz**2], + ] + ) + + # Create time array with cosine, sine, and constant components + time_array = np.array([np.cos(t * w), np.sin(t * w), np.ones_like(t)]) + + # Apply system matrix transformation, multiply by exponential decay, normalize + return (system_matrix @ time_array) * np.exp(-gamma * t) / w**2 + + +def simultaneous_expectation( + t: np.typing.NDArray, + wx: float, + wy: float, + wz: float, + gamma: float, +) -> np.typing.NDArray: + """Simulateneous computation for X, Y and Z expectation values. + + We also constrain w to be the norm of the w vector with component + wx, wy, wz. + See https://arxiv.org/pdf/2303.01427 Eq. S10. + """ + + # flatten the vectorized solution, necessary for curve_fit + return vectorized_simultaneous_expectation(t, wx, wy, wz, gamma).ravel() + + +def linear_func(x, a, b): + return a * x + b + + +def sin_func(x, a, b, omega, phi): + return a * np.sin(x * omega + phi) + b diff --git a/src/qibocal/protocols/two_qubit_interaction/cross_resonance/plotting.py b/src/qibocal/protocols/two_qubit_interaction/cross_resonance/plotting.py new file mode 100644 index 0000000000..cb98430909 --- /dev/null +++ b/src/qibocal/protocols/two_qubit_interaction/cross_resonance/plotting.py @@ -0,0 +1,311 @@ +import logging + +import numpy as np +import plotly.graph_objects as go +from plotly.subplots import make_subplots + +from qibocal.auto.operation import QubitPairId +from qibocal.protocols.utils import ( + angle_wrap, + table_dict, + table_html, +) + +from . import fitting +from .cr_parent_classes import ( + Basis, + HamiltonianTerm, + HamiltonianTomographyData, + HamiltonianTomographyResults, + SetControl, +) +from .cross_resonance_processing import compute_bloch_vector + + +def tomography_cr_plot( + data: HamiltonianTomographyData, + target: QubitPairId, + fit: HamiltonianTomographyResults | None = None, +) -> tuple[list[go.Figure], str]: + """Create cross-resonance Hamiltonian tomography plots given the acquired data and + the fit results for to a pair of qubits.""" + + fig = make_subplots( + rows=4, + cols=2, + horizontal_spacing=0.1, + vertical_spacing=0.05, + shared_xaxes=False, + shared_yaxes=False, + column_titles=["Target Qubit Evolution", "Control Qubit Evolution"], + ) + + if type(data).__name__ == "HamiltonianTomographyCRAmplData": + annotation = "CR gate amplitude [a.u.]" + else: + annotation = "CR gate duration [ns]" + + target = target if target in data.pairs else (target[1], target[0]) + for i, basis in enumerate(Basis): + for setup in SetControl: + pair_data = data.data[target[0], target[1], basis, setup] + fig.add_trace( + go.Scatter( + x=pair_data.x, + y=pair_data.prob_target, + name=f"Target when Control at {0 if setup is SetControl.Id else 1}", + showlegend=True if basis is Basis.Z else False, + legendgroup=f"Target when Control at {0 if setup is SetControl.Id else 1}", + mode="markers", + marker=dict(color="blue" if setup is SetControl.Id else "red"), + error_y=dict( + type="data", + array=pair_data.error_target, + visible=True, + ), + ), + row=i + 1, + col=1, + ) + fig.add_trace( + go.Scatter( + x=pair_data.x, + y=pair_data.prob_control, + name=f"Control in {0 if setup is SetControl.Id else 1}", + showlegend=True if basis is Basis.Z else False, + legendgroup=f"Control in {0 if setup is SetControl.Id else 1}", + mode="markers", + marker=dict(color="blue" if setup is SetControl.Id else "red"), + error_y=dict( + type="data", + array=pair_data.error_control, + visible=True, + ), + ), + row=i + 1, + col=2, + ) + if fit is not None and (*target, setup) in fit.fitted_parameters: + x = np.linspace(pair_data.x.min(), pair_data.x.max(), 100) + fig.add_trace( + go.Scatter( + x=x, + y=getattr(fitting, f"pauli_{basis.name.lower()}_expectation")( + x, + wx=fit.fitted_parameters[target[0], target[1], setup][0], + wy=fit.fitted_parameters[target[0], target[1], setup][1], + wz=fit.fitted_parameters[target[0], target[1], setup][2], + gamma=fit.fitted_parameters[target[0], target[1], setup][3], + ), + name=f"Simultaneous Fit of target when control at {0 if setup is SetControl.Id else 1}", + showlegend=True if basis is Basis.Z else False, + legendgroup=f"Simultaneous Fit target when control at {0 if setup is SetControl.Id else 1}", + mode="lines", + line=dict( + color="green" if setup is SetControl.Id else "orange", + ), + ), + row=i + 1, + col=1, + ) + + if type(fit).__name__ == "HamiltonianTomographyCRAmplResults": + fit_dict = fit.cr_amplitudes + else: + fit_dict = fit.cr_lengths + + if target in fit_dict: + fig.add_vline( + x=fit_dict[target], + line_width=2, + line_dash="dash", + line_color="red", + annotation_text=annotation, + row=i + 1, + col=1, + ) + + if fit is not None and fit.fitted_parameters: + bloch_vect_targ, bloch_fit_targ, bloch_vect_ctrl = compute_bloch_vector( + data, target, fit.fitted_parameters + ) + fig.add_traces( + [ + go.Scatter( + x=pair_data.x, + y=bloch_vect_targ, + name="Bloch vector |R(t)|", + legendgroup="Bloch vector |R(t)|", + showlegend=True, + mode="markers", + ), + go.Scatter( + x=pair_data.x, + y=bloch_vect_ctrl, + name="Bloch vector |R(t)|", + legendgroup="Bloch vector |R(t)|", + showlegend=True, + mode="markers", + ), + ], + rows=[4, 4], + cols=[1, 2], + ) + + if bloch_fit_targ is not None: + x = np.linspace(pair_data.x.min(), pair_data.x.max(), len(bloch_fit_targ)) + fig.add_trace( + go.Scatter( + x=x, + y=bloch_fit_targ, + name="Fitted Bloch vector |R(t)|", + showlegend=True, + legendgroup="Fitted Bloch vector |R(t)|", + mode="lines", + ), + row=4, + col=1, + ) + + if type(fit).__name__ == "HamiltonianTomographyCRAmplResults": + fit_dict = fit.cr_amplitudes + else: + fit_dict = fit.cr_lengths + + if target in fit_dict: + fig.add_vline( + x=fit_dict[target], + line_width=2, + line_dash="dash", + line_color="red", + annotation_text=annotation, + row=4, + col=1, + ) + + fig.update_layout(height=1000) + fig.update_xaxes(title_text=annotation, row=4, col=1) + fig.update_xaxes(title_text=annotation, row=4, col=2) + fig.update_yaxes(title_text="", range=[-1.2, 1.2], row=1) + fig.update_yaxes(title_text="", range=[-1.2, 1.2], row=2) + fig.update_yaxes(title_text="", range=[-1.2, 1.2], row=3) + fig.update_yaxes(title_text="|R(t)|", range=[-0.2, 1.2], row=4) + + return [fig], "" + + +def cancellation_calibration_plot( + data: HamiltonianTomographyData, + target: QubitPairId, + fit: HamiltonianTomographyResults | None = None, +) -> tuple[list[go.Figure], str]: + """Plot calibration results for cross-resonance Hamiltonian tomography when + tuning cancellation pulses. + + Generates plots for either amplitude or phase calibration data of cancellation pulses + in cross-resonance interactions. Fits effective Hamiltonian terms and visualizes the + results with fitted curves overlaid on experimental data. + """ + + figures = [] + fitting_report = "" + + if fit is None: + logging.warning("Fit failed, plotting only data.") + else: + if type(data).__name__ == "HamiltonianTomographyCANCAmplData": + fit_func = fitting.linear_func + x_title = "amplitude [a.u.]" + tunable_params = fit.cancellation_pulse_amplitudes[target] + plotting_terms = { + HamiltonianTerm.IX: "ampl_ix", + HamiltonianTerm.IY: "ampl_iy", + } + + if type(data).__name__ == "HamiltonianTomographyCANCPhaseData": + fit_func = fitting.sin_func + x_title = "phase [rad.]" + tunable_params = {} + tunable_params["phi0"] = fit.cancellation_pulse_phases[target]["control"] + tunable_params["phi1"] = angle_wrap( + fit.cancellation_pulse_phases[target]["control"] + - fit.cancellation_pulse_phases[target]["target"] + ) + plotting_terms = { + HamiltonianTerm.ZY: "phi0", + HamiltonianTerm.IY: "phi1", + } + + for t in HamiltonianTerm: + eff_ham_term = [] + exp_sweeper = [] + for f in fit.hamiltonian_terms[target]: + eff_ham_term.append(f[1][t]) + exp_sweeper.append(f[0]) + fig = go.Figure( + [ + go.Scatter( + x=exp_sweeper, + y=eff_ham_term, + opacity=1, + name=f"{t.name}", + showlegend=True, + legendgroup="Probability", + line=go.scatter.Line(dash="dot"), + ) + ] + ) + + if target in fit.fitted_parameters and t in fit.fitted_parameters[target]: + amp_range = np.linspace( + min(exp_sweeper), + max(exp_sweeper), + 2 * len(exp_sweeper) if len(exp_sweeper) >= 100 else 200, + ) + params = fit.fitted_parameters[target][t] + fig.add_trace( + go.Scatter( + x=amp_range, + y=fit_func(amp_range, **params), + name=f"{t.name} Fit", + mode="lines", + ), + ) + + if t in plotting_terms: + params_name = plotting_terms[t] + fig.add_vline( + x=tunable_params[params_name], + name=f"{params_name}", + line_dash="dash", + ) + + fig.update_layout( + showlegend=True, + xaxis_title=(f"{x_title}"), + yaxis_title="Interaction strength [MHz]", + ) + + fitting_report = table_html( + table_dict( + 8 * [target], + ( + [ + f"{term.name}: Fitted_parameters" + for term in HamiltonianTerm + ] + + [k for k in tunable_params.keys()] + ), + ( + [ + fit.fitted_parameters[target][term] + for term in HamiltonianTerm + ] + + [v for v in tunable_params.values()] + ), + ) + ) + + figures.append(fig) + + return figures, fitting_report diff --git a/src/qibocal/protocols/two_qubit_interaction/cross_resonance/utils.py b/src/qibocal/protocols/two_qubit_interaction/cross_resonance/utils.py new file mode 100644 index 0000000000..55ca375603 --- /dev/null +++ b/src/qibocal/protocols/two_qubit_interaction/cross_resonance/utils.py @@ -0,0 +1,430 @@ +import numpy as np +from qibolab import Delay, Platform, Pulse, PulseSequence, Rectangular, VirtualZ + +from qibocal.auto.operation import ( + QubitId, + QubitPairId, +) +from qibocal.calibration import CalibrationPlatform +from qibocal.protocols.utils import angle_wrap +from qibocal.update import cnot_sequence, replace + +from .cr_parent_classes import Basis, SetControl + + +def cross_resonance_pulses( + platform: CalibrationPlatform, control: QubitId, target: QubitId +) -> tuple[Pulse | None, Pulse | None]: + """Retrieve cross-resonance (CR) pulse parameters from platform calibration. + This function extracts the CR pulse and its corresponding cancellation pulse + from the platform's native gates configuration for a given control-target qubit pair. + """ + + # extracting the channels involved in the CR gate + target_channel, _ = platform.parameters.native_gates.single_qubit[target].RX()[0] + cr_channel = platform.qubits[control].drive_extra[target] + + # CNOT sequence present in the platform parameters.json + cnot_cal_seq = platform.parameters.native_gates.two_qubit[(control, target)].CNOT + + cr_params = None + canc_params = None + if cnot_cal_seq is None and len(cnot_cal_seq) != 0: + for p in cnot_cal_seq[2:]: + if cr_params is None and p[0] == cr_channel and isinstance(p[1], Pulse): + cr_params = p[1] + if ( + canc_params is None + and p[1] == target_channel + and isinstance(p[1], Pulse) + ): + canc_params = p[1] + + if cr_params is not None and canc_params is not None: + # here we identified both the CR and cancellation pulse, + # we can stop looking at the remaining pulses + break + + return cr_params, canc_params + + +def ro_delay_range( + cr_pulse_duration_range: tuple[float, float, float], + echo: bool, + cntl_setup: SetControl, + control: QubitId, + platform: Platform, +) -> tuple[float, float, float]: + """Delay range for RO pulses. + + add the number of the pi-pulses if we are in echo mode + or if we want to set the control to 1. + num_cr_pulses is the number of cross resonance pulses: + == 1 if there is no echo sequence + == 2 if there is echo sequence + num_pi_pulses is the number of pi-pulses: + == 0 if there is no echo and control at 0 + == 1 if there is no echo and control at 1 + == 2 if there is echo and control at 0 + == 3 if there is echo and control at 1 + """ + + # Calculate number of CR pulses: 1 if no echo, 2 if echo is enabled + num_cr_pulses = 1 + int(echo) + # Calculate number of pi-pulses: 2 per echo, plus 1 if control is set to X (SetControl.X) + num_pi_pulses = 2 * int(echo) + int(cntl_setup == SetControl.X) + + # Get the duration of the single-qubit RX (pi/2) pulse for the control qubit + pi_pulse_duration = platform.natives.single_qubit[control].RX()[0][1].duration + + # the total duration of the cross resonance sequence is given by: + # num_cr_pulses * cr_duration + num_pi_pulses * pi_pulse_duration + + # Calculate minimum delay: start of CR range + all pi-pulses + tot_delay_start = ( + num_cr_pulses * cr_pulse_duration_range[0] + num_pi_pulses * pi_pulse_duration + ) + # Calculate maximum delay: end of CR range + all pi-pulses + tot_delay_end = ( + num_cr_pulses * cr_pulse_duration_range[1] + num_pi_pulses * pi_pulse_duration + ) + # Calculate step size: CR pulse range step multiplied by number of CR pulses + tot_delay_step = num_cr_pulses * cr_pulse_duration_range[2] + + return (tot_delay_start, tot_delay_end, tot_delay_step) + + +def cross_res_sequence( + platform: Platform, + control: QubitId, + target: QubitId, + duration: float, + control_amplitude: float, + control_phase: float, + target_amplitude: float | None, + target_phase: float, + echo: bool, + interpolated_sweeper: bool = False, +) -> tuple[PulseSequence, list[Pulse], list[Pulse], list[Delay]]: + """Creates sequence for CR experiment on ``control`` and ``target`` qubits. + + With ``setup`` it is possible to set the control qubit to 1 or keep it at 0. + If ``echo`` is set to ``True`` a ECR gate will be played. + With ``basis`` it is possible to set the measurement basis. If it is not provided + the default is Z.""" + + cr_control_pulses = [] + cr_target_pulses = [] + # delays introduced by the cross resonance sequence + cr_delays = [] + + sequence = PulseSequence() + + natives_control = platform.natives.single_qubit[control] + natives_target = platform.natives.single_qubit[target] + + control_drive_channel, control_drive_pulse = natives_control.RX()[0] + target_drive_channel, _ = natives_target.RX()[0] + cr_channel = platform.qubits[control].drive_extra[target] + + cr_drive_pulse, target_drive_pulse = cross_resonance_pulses( + platform, control, target + ) + + cr_drive_pulse = Pulse( + duration=duration, + amplitude=control_amplitude, + relative_phase=angle_wrap(control_phase), + envelope=Rectangular() if cr_drive_pulse is None else cr_drive_pulse.envelope, + ) + cr_control_pulses.append(cr_drive_pulse) + + if target_amplitude is None: + target_drive_pulse = Delay(duration=duration) + else: + target_drive_pulse = Pulse( + duration=duration, + amplitude=target_amplitude, + relative_phase=angle_wrap(target_phase), + envelope=Rectangular() + if target_drive_pulse is None + else target_drive_pulse.envelope, + ) + cr_target_pulses.append(target_drive_pulse) + + # delays introduced by cross resonance pulses + cr_delays.append(Delay(duration=cr_drive_pulse.duration)) + # first cross resonance pulse + sequence.append((cr_channel, cr_drive_pulse)) + sequence.append((target_drive_channel, target_drive_pulse)) + if interpolated_sweeper: + _ = sequence.align([control_drive_channel, cr_channel, target_drive_channel]) + else: + sequence.append((control_drive_channel, cr_delays[-1])) + + if echo: + # phase-flipped cross resonance pulses + cr_drive_pulse_flipped = replace( + cr_drive_pulse.new(), + amplitude=-cr_drive_pulse.amplitude, + ) + cr_control_pulses.append(cr_drive_pulse_flipped) + + if target_amplitude is None: + target_pulse_flipped = target_drive_pulse.new() + else: + target_pulse_flipped = replace( + target_drive_pulse.new(), + amplitude=-target_drive_pulse.amplitude, + ) + cr_target_pulses.append(target_pulse_flipped) + + # delay introduced by echo sequence + echo_delay = Delay(duration=control_drive_pulse.duration) + + # first echo pulse + sequence.append((control_drive_channel, control_drive_pulse)) + sequence.append((cr_channel, echo_delay)) + sequence.append((target_drive_channel, echo_delay)) + + # delays introduced by cross resonance pulses + cr_delays.append(Delay(duration=cr_drive_pulse.duration)) + # second cross resonance pulse with flipped phase + sequence.append((cr_channel, cr_drive_pulse_flipped)) + sequence.append((target_drive_channel, target_pulse_flipped)) + if interpolated_sweeper: + _ = sequence.align( + [control_drive_channel, cr_channel, target_drive_channel] + ) + else: + sequence.append((control_drive_channel, cr_delays[-1])) + + # second echo pulse + sequence.append((control_drive_channel, control_drive_pulse)) + sequence.append((cr_channel, echo_delay)) + sequence.append((target_drive_channel, echo_delay)) + + return sequence, cr_control_pulses, cr_target_pulses, cr_delays + + +def appending_ro_sequence( + platform: Platform, + control: QubitId, + target: QubitId, + exp_sequence: PulseSequence, + basis: Basis, + interpolated_sweeper: bool, +) -> tuple[PulseSequence, list[Delay]]: + """Append a readout pulse sequence for two qubits with a specified delay. + + This function constructs a pulse sequence that applies readout operations to both + control and target qubits with an initial delay on each readout channel. + """ + natives_control = platform.natives.single_qubit[control] + natives_target = platform.natives.single_qubit[target] + + # platform acquisition channels and pulses for control and target qubits + ro_channel_control, ro_pulse_control = natives_control.MZ()[0] + ro_channel_target, ro_pulse_target = natives_target.MZ()[0] + + # switching measurement basis for target and control lines and align all the others + target_drive_channel, _ = natives_target.RX()[0] + control_drive_channel, _ = natives_control.RX()[0] + # delay + ro_delays = [Delay(duration=exp_sequence.duration) for _ in range(2)] + + # basis on which measuring the state + if basis == Basis.X: + tom_angle = np.pi / 2 + else: # only relevant for Y + tom_angle = 0 + + target_rotation = natives_target.R(theta=np.pi / 2, phi=tom_angle)[0][1] + control_rotation = natives_control.R(theta=np.pi / 2, phi=tom_angle)[0][1] + + if interpolated_sweeper: + # if using interpolated_sweepers I need all the channels to be aligned + cr_channel = platform.qubits[control].drive_extra[target] + + # align all the channels of the sequence + _ = exp_sequence.align( + [ + cr_channel, + target_drive_channel, + control_drive_channel, + ro_channel_control, + ro_channel_target, + ] + ) + else: + # delay to align the readout pulses of control and target qubits + # done in every case, even when measuring Z-basis (no additional pulse on target) + # to have the same timing for all the measurements + flip_duration = max(target_rotation.duration, control_rotation.duration) + flip_delay = Delay( + duration=flip_duration, + ) + exp_sequence.append((ro_channel_control, flip_delay)) + exp_sequence.append((ro_channel_target, flip_delay.new())) + + # wait the whole cr-sequence duration before starting the readout pulses + exp_sequence.append((ro_channel_target, ro_delays[0])) + exp_sequence.append((ro_channel_control, ro_delays[1])) + + # applying the tomography rotations on X and Y basis + if basis != Basis.Z: + exp_sequence += PulseSequence( + [ + (target_drive_channel, target_rotation), + (control_drive_channel, control_rotation), + ] + ) + + # adding readout pulses + exp_sequence.append((ro_channel_target, ro_pulse_target)) + exp_sequence.append((ro_channel_control, ro_pulse_control)) + + return exp_sequence, ro_delays + + +def cross_resonance_experiment( + platform: Platform, + pair_list: list[QubitPairId], + duration: float | dict[QubitPairId, float], + ctrl_ampl: float | dict[QubitPairId, float], + ctrl_phase: float | dict[QubitPairId, float], + targ_ampl: float | dict[QubitPairId, float | None] | None, + targ_phase: float | dict[QubitPairId, float], + basis: Basis, + setup: SetControl, + echo: bool = False, + interpolated_sweeper: bool = False, +) -> tuple[ + PulseSequence, + dict[QubitPairId, Pulse], + dict[QubitPairId, Pulse], + dict[QubitPairId, Delay], + dict[QubitPairId, Delay], +]: + """Build the pulse sequence for a cross-resonance experiment. + + The sequence is created for one or more control-target qubit pairs. The + control qubit can be prepared in either |0> or |1> depending on ``setup``. + If ``echo`` is True, an echoed cross-resonance sequence is generated. + The target measurement basis is selected with ``basis``. + """ + + parallel_cr_sequences = PulseSequence() + parallel_cr_pulses: dict[QubitPairId, Pulse] = {} + parallel_cr_target_pulses: dict[QubitPairId, Pulse] = {} + parallel_cr_delays: dict[QubitPairId, Delay] = {} + parallel_ro_delays: dict[QubitPairId, Delay] = {} + for pair in pair_list: + control, target = pair + + # adding pi-pulse if we want to set control to 1 + cntl_setup_sequence = PulseSequence() + if setup == SetControl.X: + control_drive_channel, control_drive_pulse = platform.natives.single_qubit[ + control + ].RX()[0] + target_drive_channel, _ = platform.natives.single_qubit[target].RX()[0] + cr_channel = platform.qubits[control].drive_extra[target] + + control_delay = Delay(duration=control_drive_pulse.duration) + + cntl_setup_sequence.append((control_drive_channel, control_drive_pulse)) + cntl_setup_sequence.append((target_drive_channel, control_delay)) + cntl_setup_sequence.append((cr_channel, control_delay)) + + cr_sequence, cr_pulses, cr_target_pulses, cr_delays = cross_res_sequence( + platform=platform, + control=control, + target=target, + duration=duration[pair] if isinstance(duration, dict) else duration, + control_amplitude=ctrl_ampl[pair] + if isinstance(ctrl_ampl, dict) + else ctrl_ampl, + control_phase=ctrl_phase[pair] + if isinstance(ctrl_phase, dict) + else ctrl_phase, + target_amplitude=targ_ampl[pair] + if isinstance(targ_ampl, dict) + else targ_ampl, + target_phase=targ_phase[pair] + if isinstance(targ_phase, dict) + else targ_phase, + echo=echo, + interpolated_sweeper=interpolated_sweeper, + ) + cr_sequence = cntl_setup_sequence | cr_sequence + + total_sequence, ro_delays = appending_ro_sequence( + platform=platform, + control=control, + target=target, + exp_sequence=cr_sequence, + basis=basis, + interpolated_sweeper=interpolated_sweeper, + ) + + parallel_cr_sequences += total_sequence + parallel_cr_pulses |= {pair: cr_pulses} + parallel_cr_target_pulses |= {pair: cr_target_pulses} + parallel_cr_delays |= {pair: cr_delays} + parallel_ro_delays |= {pair: ro_delays} + + return ( + parallel_cr_sequences, + parallel_cr_pulses, + parallel_cr_target_pulses, + parallel_cr_delays, + parallel_ro_delays, + ) + + +def update_cnot_from_fit( + platform: CalibrationPlatform, + pair: QubitPairId, + cr_duration: float, + cr_ampl: float, + control_phase: float, + canc_ampl: float | None, + canc_phase: float, + echo_flag: bool, +) -> None: + """Update CNOT gate calibration from cross-resonance fit parameters. + + Constructs and updates the CNOT gate using cross-resonance pulses with + fitted parameters, including single-qubit rotations and virtual Z phases. + """ + ctrl, targ = pair + + cr_seq, _, _, _ = cross_res_sequence( + platform=platform, + control=ctrl, + target=targ, + duration=cr_duration, + control_amplitude=cr_ampl, + control_phase=control_phase, + target_amplitude=canc_ampl, + target_phase=canc_phase, + echo=echo_flag, + ) + + target_single_qubit_operation = ( + platform.qubits[targ].drive, + platform.natives.single_qubit[targ].R(theta=np.pi / 2, phi=0)[0][1], + ) + + control_single_qubit_operation = ( + platform.qubits[ctrl].drive, + VirtualZ(phase=np.pi / 2), + ) + + new_cr_seq = ( + PulseSequence([control_single_qubit_operation, target_single_qubit_operation]) + | cr_seq + ) + + cnot_sequence(new_cr_seq, platform, pair) diff --git a/src/qibocal/protocols/utils.py b/src/qibocal/protocols/utils.py index 3ca859ea3b..1b251e1a9f 100644 --- a/src/qibocal/protocols/utils.py +++ b/src/qibocal/protocols/utils.py @@ -2,6 +2,7 @@ from collections.abc import Sequence from colorsys import hls_to_rgb from enum import Enum +from typing import Any import numpy as np import pandas as pd @@ -1013,24 +1014,6 @@ def reshaping_raw_signal(x, y, z): return x_, y_, z_ -def guess_period(x, y): - """Return fft period estimation given a sinusoidal plot.""" - - fft = np.fft.rfft(y) - fft_freqs = np.fft.rfftfreq(len(y), d=(x[1] - x[0])) - mags = abs(fft) - mags[0] = 0 - local_maxima, _ = find_peaks(mags) - if len(local_maxima) > 0: - return 1 / fft_freqs[np.argmax(mags)] - return None - - -def fallback_period(period): - """Function to estimate period if guess_period fails.""" - return period if period is not None else 4 - - def angle_wrap(angle: float): """Wrap an angle from [-np.inf,np.inf] into the [0,2*np.pi] domain""" return angle % (2 * np.pi) @@ -1064,3 +1047,143 @@ def baseline_als(data: NDArray, lamda: float, p: float, niter: int = 10) -> NDAr z = sparse.linalg.spsolve(a, b) weights = p * (data > z) + (1 - p) * (data < z) return z + + +def guess_frequency(x: NDArray, y: NDArray, axis: int = -1) -> NDArray: + """Return FFT frequency estimation given a sinusoidal plot. + + Computes the dominant frequency component from the FFT of the signal. + axis represents the dimension along which to compute the FFT. + This parameter specifies which dimension of y contains the sinusoidal + data to analyze. For multi-dimensional arrays, the FFT is computed + along this axis, and the dominant frequency is determined for each + slice along other dimensions. + """ + assert x.ndim == 1, f"Expected 1D array, got array with shape {x.shape}" + + y = np.moveaxis(y, axis, -1) + fft = np.fft.rfft(y, axis=-1) + fft_freqs = np.fft.rfftfreq(y.shape[-1], d=(x[1] - x[0])) + mags = np.abs(fft) + mags[:, 0] = 0 + + selected_freqs = fft_freqs[np.argmax(mags, axis=-1)] + + return np.moveaxis(selected_freqs, -1, axis) + + +def fallback_frequency(frequency: NDArray) -> NDArray: + """Return a numeric frequency array with NaNs replaced by a fallback value. + + This helper is used to ensure a valid numeric frequency is available for + downstream processing when the frequency estimation returns missing values. + """ + assert frequency.ndim <= 1, ( + f"Expected 1D array or scalar, got array with shape {frequency.shape}" + ) + + return np.where(np.isnan(frequency), 0.25, frequency) + + +def quinn_fernandes_algorithm( + signal: Any, + x: Any, + axis: int = -1, + speedup_flag: bool = False, + iterations: int = 100, + tol: int = 1e-8, +) -> NDArray: + """This is a custom implementation of the Quinn-Fernandes algorithm. + We compute the signal sampling rate from :param:x, hence this function assumes x to be + ordered. + + The Quinn–Fernandes method is a high-accuracy frequency estimator based on + phase interpolation of the discrete Fourier transform (DFT). It refines the + peak frequency obtained from the FFT by analyzing the phase evolution of the + complex spectrum, achieving super-resolution beyond the FFT bin spacing. + If :const:`speedup_flag` is set to `True`, the algorithm will change the updating rule, + can lead to faster convergence, especially when the initial guess is close to the true frequency. + Link for the original paper: https://www.jstor.org/stable/2337018?seq=3 + """ + + fs = 1 / abs(x[0] - x[1]) + + if not isinstance(x, np.ndarray): + x = np.array(x) + + if not isinstance(signal, np.ndarray): + signal = np.array(signal) + + if signal.ndim == 1: + signal = signal[np.newaxis, :] + + omegas = 2 * np.pi * fallback_frequency(guess_frequency(x, signal, axis=axis)) + alpha = 2 * np.cos(omegas) + + signal = signal - np.mean(signal, axis=axis, keepdims=True) + + sig_shape = list(signal.shape) + sig_shape[axis] += 2 + buffer_beta = [] + for _ in range(iterations): + xi = np.zeros(sig_shape) + for t in range(2, xi.shape[axis]): + xi[..., t] = signal[..., t - 2] + alpha * xi[..., t - 1] - xi[..., t - 2] + + beta = np.sum((xi[..., 2:] + xi[..., :-2]) * xi[..., 1:-1], axis=axis) / np.sum( + xi[..., :-1] ** 2, axis=axis + ) + beta[np.isnan(beta)] = 0 + if len(buffer_beta) >= 5: + buffer_beta.pop(0) + buffer_beta.append(beta) + if np.all(np.abs(np.mean(buffer_beta, axis=0) - alpha) < tol): + alpha = beta + break + + if speedup_flag: + alpha = beta + else: + alpha = 2 * beta - alpha + + alpha = np.clip(alpha, -2, 2) + omega_est = np.arccos(alpha / 2) + med_omega = np.median(omega_est) + + return med_omega * fs + + +def guess_period( + x: NDArray, y: NDArray, axis: int = -1, speedup_flag: bool = True +) -> NDArray: + """Estimate the period(s) of a (set of) sinusoidal signal(s). + + This is a thin wrapper around :func:`guess_frequency` that returns the + reciprocal of the estimated dominant frequency. For multi-dimensional + signals, the period is computed along ``axis`` of ``y`` and returned with + the same shape semantics as :func:`guess_frequency`. + """ + + return ( + 2 + * np.pi + / quinn_fernandes_algorithm( + signal=y, + x=x, + axis=axis, + speedup_flag=speedup_flag, + ) + ) + + +def fallback_period(period: NDArray) -> NDArray: + """Return a numeric period array with NaNs replaced by a fallback value. + + This helper is used to ensure a valid numeric period is available for + downstream processing when the period estimation returns missing values. + """ + assert period.ndim <= 1, ( + f"Expected 1D array or scalar, got array with shape {period.shape}" + ) + + return np.where(np.isnan(period), 4, period) diff --git a/src/qibocal/result.py b/src/qibocal/result.py index f1be8afc29..7737868135 100644 --- a/src/qibocal/result.py +++ b/src/qibocal/result.py @@ -76,6 +76,6 @@ def phase(iq: npt.NDArray): return np.unwrap(np.arctan2(iq[..., 0], iq[..., 1])) -def probability(values: npt.NDArray, state: int = 0): +def probability(values: npt.NDArray, state: int = 0) -> npt.NDArray: """Return the statistical frequency of the specified state.""" - return np.sum(values == state, axis=0) / len(values) + return np.mean(values == state, axis=0) diff --git a/src/qibocal/update.py b/src/qibocal/update.py index 8903edf47f..d1f9a3de41 100644 --- a/src/qibocal/update.py +++ b/src/qibocal/update.py @@ -4,7 +4,7 @@ from typing import Literal import numpy as np -from pydantic import BaseModel +from pydantic import BaseModel, TypeAdapter from qibolab import Platform, PulseSequence, VirtualZ from qibocal.auto.operation import QubitId, QubitPairId @@ -18,6 +18,14 @@ ] +_QPAIR_DUMPER = TypeAdapter(QubitPairId) + + +def _dump_pair(pair: QubitPairId) -> str: + """Convert a QubitPairId to its string representation.""" + return _QPAIR_DUMPER.dump_python(pair) + + def replace(model: BaseModel, **update): """Replace interface for pydantic models.""" return model.model_copy(update=update) @@ -153,40 +161,43 @@ def virtual_phases( for qubit, phase in phases.items(): new_native.append((platform.qubits[qubit].drive, VirtualZ(phase=phase))) - platform.update( - {f"native_gates.two_qubit.{f'{pair[0]}-{pair[1]}'}.{native}": new_native} - ) + platform.update({f"native_gates.two_qubit.{_dump_pair(pair)}.{native}": new_native}) def CZ_duration(duration: float, platform: Platform, pair: QubitPairId): """Update CZ duration for specific pair.""" platform.update( - {f"native_gates.two_qubit.{f'{pair[0]}-{pair[1]}'}.CZ.0.1.duration": duration} + {f"native_gates.two_qubit.{_dump_pair(pair)}.CZ.0.1.duration": duration} ) def CZ_amplitude(amp: float, platform: Platform, pair: QubitPairId): """Update CZ amplitude for specific pair.""" platform.update( - {f"native_gates.two_qubit.{f'{pair[0]}-{pair[1]}'}.CZ.0.1.amplitude": amp} + {f"native_gates.two_qubit.{_dump_pair(pair)}.CZ.0.1.amplitude": amp} ) def iSWAP_duration(duration: float, platform: Platform, pair: QubitPairId): """Update iSWAP_duration duration for specific pair.""" platform.update( - {f"native_gates.two_qubit.{f'{pair[0]}-{pair[1]}'}.CZ.0.1.duration": duration} + {f"native_gates.two_qubit.{_dump_pair(pair)}.CZ.0.1.duration": duration} ) def iSWAP_amplitude(amp: float, platform: Platform, pair: QubitPairId): """Update iSWAP_duration amplitude for specific pair.""" platform.update( - {f"native_gates.two_qubit.{f'{pair[0]}-{pair[1]}'}.CZ.0.1.amplitude": amp} + {f"native_gates.two_qubit.{_dump_pair(pair)}.CZ.0.1.amplitude": amp} ) -def t1(t1: float, platform: Platform, qubit: QubitId): +def cnot_sequence(cr_sequence: PulseSequence, platform: Platform, pair: QubitPairId): + """Update CNOT gate for specific pair.""" + platform.update({f"native_gates.two_qubit.{_dump_pair(pair)}.CNOT": cr_sequence}) + + +def t1(t1: int, platform: Platform, qubit: QubitId): """Update t1 value in platform for specific qubit.""" platform.calibration.single_qubits[qubit].t1 = tuple(t1) diff --git a/tests/platforms/mock/parameters.json b/tests/platforms/mock/parameters.json index 0f087ccded..cf0e074555 100644 --- a/tests/platforms/mock/parameters.json +++ b/tests/platforms/mock/parameters.json @@ -8,6 +8,10 @@ "kind": "iq", "frequency": 4000000000.0 }, + "01/drive": { + "kind": "iq", + "frequency": 4200000000.0 + }, "1/drive": { "kind": "iq", "frequency": 4200000000.0 diff --git a/tests/platforms/mock/platform.py b/tests/platforms/mock/platform.py index f336d9aa87..9bda139fe7 100644 --- a/tests/platforms/mock/platform.py +++ b/tests/platforms/mock/platform.py @@ -16,19 +16,28 @@ def create_mock_hardware() -> Hardware: channels = {} # attach the channels pump_name = "twpa_pump" - for q in range(2): - drive12 = f"{q}/drive12" - qubits[q] = qubit = Qubit.default(q, drive_extra={(1, 2): drive12}) - channels |= { - qubit.probe: IqChannel(mixer=None, lo="01/probe_lo"), - qubit.acquisition: AcquisitionChannel( - twpa_pump=pump_name, probe=qubit.probe - ), - qubit.drive: IqChannel(mixer=None, lo=f"{q}/drive_lo"), - drive12: IqChannel(mixer=None, lo=f"{q}/drive_lo"), - qubit.flux: DcChannel(), - } - + qubits[0] = Qubit.default(0, drive_extra={(1, 2): "0/drive12", 1: "01/drive"}) + + channels |= { + qubits[0].probe: IqChannel(mixer=None, lo="01/probe_lo"), + qubits[0].acquisition: AcquisitionChannel( + twpa_pump=pump_name, probe=qubits[0].probe + ), + qubits[0].drive: IqChannel(mixer=None, lo="0/drive_lo"), + qubits[0].drive_extra[1, 2]: IqChannel(mixer=None, lo="0/drive_lo"), + qubits[0].drive_extra[1]: IqChannel(mixer=None, lo="0/drive_lo"), + qubits[0].flux: DcChannel(), + } + qubits[1] = Qubit.default(1, drive_extra={(1, 2): "1/drive12"}) + channels |= { + qubits[1].probe: IqChannel(mixer=None, lo="01/probe_lo"), + qubits[1].acquisition: AcquisitionChannel( + twpa_pump=pump_name, probe=qubits[1].probe + ), + qubits[1].drive: IqChannel(mixer=None, lo="1/drive_lo"), + qubits[1].drive_extra[1, 2]: IqChannel(mixer=None, lo="1/drive_lo"), + qubits[1].flux: DcChannel(), + } couplers = {} couplers["01"] = coupler = Qubit(flux="coupler_01/flux") channels |= {coupler.flux: DcChannel()} diff --git a/tests/runcards/protocols.yml b/tests/runcards/protocols.yml index 06d2e9f74a..1d24143843 100644 --- a/tests/runcards/protocols.yml +++ b/tests/runcards/protocols.yml @@ -678,9 +678,9 @@ actions: amplitude_step: 0.01 duration_min: 10 duration_max: 50 - duration_step: 10 - nshots: 10 + duration_step: 1 native: CZ + nshots: 1000 parking: True - id: chevron cz signal @@ -935,3 +935,29 @@ actions: twpa_pow_width: 10 twpa_pow_step: 3 nshots: 20000 + + - id: CR tomography amplitude + operation: cr_amplitude + targets: [[0, 1]] + parameters: + amplitude_range: [0.01, 0.1, 0.01] + echo: true + nshots: 2000 + pulse_duration: 10 + + - id: CR tomography length + operation: cr_length + targets: [[0, 1]] + parameters: + duration_range: [10, 20, 1] + pulse_amplitude: 0.1 + nshots: 2000 + + - id: CR tomography interpolated + operation: cr_length + targets: [[0, 1]] + parameters: + duration_range: [10, 20, 1] + nshots: 2000 + pulse_amplitude: 0.1 + interpolated_sweeper: true diff --git a/tests/test_fft.py b/tests/test_fft.py index 2eadc5b32c..09eeb50296 100644 --- a/tests/test_fft.py +++ b/tests/test_fft.py @@ -4,14 +4,23 @@ from qibocal.protocols.utils import guess_period -@pytest.mark.parametrize("period", [0.2, 0.5, 1, 2]) +@pytest.mark.parametrize("period", [0.05, 0.08, 0.2, 0.5, 1, 1.2]) def test_guess_period(period): """Testing guess period function with different periods.""" t = np.linspace(0, 1, 100, endpoint=False) signal = np.sin(2 * np.pi / period * t) extracted_period = guess_period(t, signal) - if period < 1: - np.testing.assert_allclose(period, extracted_period) - else: - extracted_period == 1 + assert ( + pytest.approx( + extracted_period, + rel=3e-2 + if period < 0.3 + else 20e-2 + if period < 0.8 + else 0.5 + if period < 1.5 + else None, # raise error, window too small + ) + == period + ) diff --git a/tests/test_fit_functions.py b/tests/test_fit_functions.py index 4522c03dad..781ddd465f 100644 --- a/tests/test_fit_functions.py +++ b/tests/test_fit_functions.py @@ -8,9 +8,11 @@ fit_amplitude_function as rabi_fit_amplitude_function, ) from qibocal.protocols.rabi.utils import fit_length_function as rabi_fit_length_function +from qibocal.protocols.rabi.utils import ( + rabi_initial_guess, +) from qibocal.protocols.ramsey.processing import fitting as ramsey_fitting from qibocal.protocols.ramsey.processing import process_fit as ramsey_process_fit -from qibocal.protocols.utils import fallback_period, guess_period def test_ramsey_fit(): @@ -80,28 +82,18 @@ def test_rabi_fit(): signal = (raw_signal - sig_min) / (sig_max - sig_min) x_lims = (x_min, x_max) signal_lims = (sig_min, sig_max) - + signal_flag = True else: signal = raw_signal x = raw_x x_lims = (None, None) signal_lims = (None, None) - - period = fallback_period(guess_period(x, signal)) - median_sig = np.median(signal) - q80 = np.quantile(signal, 0.8) - q20 = np.quantile(signal, 0.2) - amplitude_guess = abs(q80 - q20) + signal_flag = False if "amp" in str_sub: - signal_flag = "signal" in str_sub - pguess = [ - median_sig, - amplitude_guess, - period, - np.pi / 2 if signal_flag else np.pi, - ] - _, _, pi_pulse_parameter = rabi_fit_amplitude_function( + pguess = rabi_initial_guess(x, signal, "amplitude", signal_flag) + + fit_params, _, pi_pulse_parameter = rabi_fit_amplitude_function( x, signal, pguess, @@ -121,15 +113,9 @@ def test_rabi_fit(): assert math.isclose(true_amplitude, new_amplitude, rel_tol=2.5e-2) if "length" in str_sub: - signal_flag = "signal" in str_sub - pguess = [ - median_sig, - amplitude_guess, - period, - np.pi / 2 if signal_flag else np.pi, - 0, - ] - _, _, pi_pulse_parameter = rabi_fit_length_function( + pguess = rabi_initial_guess(x, signal, "length", signal_flag) + + fit_params, _, pi_pulse_parameter = rabi_fit_length_function( x, signal, pguess,