Skip to content

Commit

Permalink
Merge pull request #567 from qiboteam/rabi_prob
Browse files Browse the repository at this point in the history
Rabi with probabilities
  • Loading branch information
andrea-pasquale authored Oct 26, 2023
2 parents 099476b + 9f28412 commit 6fa349d
Show file tree
Hide file tree
Showing 11 changed files with 569 additions and 111 deletions.
4 changes: 4 additions & 0 deletions src/qibocal/protocols/characterization/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -26,8 +26,10 @@
from .qubit_spectroscopy_ef import qubit_spectroscopy_ef
from .qutrit_classification import qutrit_classification
from .rabi.amplitude import rabi_amplitude
from .rabi.amplitude_msr import rabi_amplitude_msr
from .rabi.ef import rabi_amplitude_ef
from .rabi.length import rabi_length
from .rabi.length_msr import rabi_length_msr
from .rabi.length_sequences import rabi_length_sequences
from .ramsey import ramsey
from .ramsey_msr import ramsey_msr
Expand Down Expand Up @@ -60,6 +62,8 @@ class Operation(Enum):
rabi_amplitude = rabi_amplitude
rabi_length = rabi_length
rabi_length_sequences = rabi_length_sequences
rabi_amplitude_msr = rabi_amplitude_msr
rabi_length_msr = rabi_length_msr
ramsey = ramsey
ramsey_msr = ramsey_msr
ramsey_sequences = ramsey_sequences
Expand Down
71 changes: 34 additions & 37 deletions src/qibocal/protocols/characterization/rabi/amplitude.py
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,7 @@
from qibocal.auto.operation import Data, Parameters, Qubits, Results, Routine
from qibocal.config import log

from ..utils import chi2_reduced
from . import utils


Expand All @@ -36,16 +37,17 @@ class RabiAmplitudeParameters(Parameters):
class RabiAmplitudeResults(Results):
"""RabiAmplitude outputs."""

amplitude: dict[QubitId, float] = field(metadata=dict(update="drive_amplitude"))
amplitude: dict[QubitId, tuple[float, Optional[float]]]
"""Drive amplitude for each qubit."""
length: dict[QubitId, float] = field(metadata=dict(update="drive_length"))
length: dict[QubitId, tuple[float, Optional[float]]]
"""Drive pulse duration. Same for all qubits."""
fitted_parameters: dict[QubitId, dict[str, float]]
"""Raw fitted parameters."""
chi2: dict[QubitId, tuple[float, Optional[float]]] = field(default_factory=dict)


RabiAmpType = np.dtype(
[("amp", np.float64), ("msr", np.float64), ("phase", np.float64)]
[("amp", np.float64), ("prob", np.float64), ("error", np.float64)]
)
"""Custom dtype for rabi amplitude."""

Expand Down Expand Up @@ -100,9 +102,6 @@ def _acquisition(
type=SweeperType.FACTOR,
)

# create a DataUnits object to store the results,
# DataUnits stores by default MSR, phase, i, q
# additionally include qubit drive pulse amplitude
data = RabiAmplitudeData(durations=durations)

# sweep the parameter
Expand All @@ -111,21 +110,20 @@ def _acquisition(
ExecutionParameters(
nshots=params.nshots,
relaxation_time=params.relaxation_time,
acquisition_type=AcquisitionType.INTEGRATION,
averaging_mode=AveragingMode.CYCLIC,
acquisition_type=AcquisitionType.DISCRIMINATION,
averaging_mode=AveragingMode.SINGLESHOT,
),
sweeper,
)
for qubit in qubits:
# average msr, phase, i and q over the number of shots defined in the runcard
result = results[ro_pulses[qubit].serial]
prob = results[qubit].probability(state=1)
data.register_qubit(
RabiAmpType,
(qubit),
dict(
amp=qd_pulses[qubit].amplitude * qd_pulse_amplitude_range,
msr=result.magnitude,
phase=result.phase,
prob=prob.tolist(),
error=np.sqrt(prob * (1 - prob) / params.nshots).tolist(),
),
)
return data
Expand All @@ -137,19 +135,13 @@ def _fit(data: RabiAmplitudeData) -> RabiAmplitudeResults:

pi_pulse_amplitudes = {}
fitted_parameters = {}
chi2 = {}

for qubit in qubits:
qubit_data = data[qubit]

rabi_parameter = qubit_data.amp
voltages = qubit_data.msr

y_min = np.min(voltages)
y_max = np.max(voltages)
x_min = np.min(rabi_parameter)
x_max = np.max(rabi_parameter)
x = (rabi_parameter - x_min) / (x_max - x_min)
y = (voltages - y_min) / (y_max - y_min)
x = qubit_data.amp
y = qubit_data.prob

# Guessing period using fourier transform
ft = np.fft.rfft(y)
Expand All @@ -158,9 +150,9 @@ def _fit(data: RabiAmplitudeData) -> RabiAmplitudeResults:
index = local_maxima[0] if len(local_maxima) > 0 else None
# 0.5 hardcoded guess for less than one oscillation
f = x[index] / (x[1] - x[0]) if index is not None else 0.5
pguess = [0.5, 1, f, np.pi / 2]
pguess = [0.5, 0.5, 1 / f, np.pi / 2]
try:
popt, _ = curve_fit(
popt, perr = curve_fit(
utils.rabi_amplitude_fit,
x,
y,
Expand All @@ -170,29 +162,34 @@ def _fit(data: RabiAmplitudeData) -> RabiAmplitudeResults:
[0, 0, 0, -np.pi],
[1, 1, np.inf, np.pi],
),
sigma=qubit_data.error,
)
translated_popt = [
y_min + (y_max - y_min) * popt[0],
(y_max - y_min) * popt[1],
popt[2] / (x_max - x_min),
popt[3] - 2 * np.pi * x_min / (x_max - x_min) * popt[2],
]
pi_pulse_parameter = np.abs((1.0 / translated_popt[2]) / 2)
perr = np.sqrt(np.diag(perr))
pi_pulse_parameter = np.abs(popt[2] / 2)

except:
log.warning("rabi_fit: the fitting was not succesful")
pi_pulse_parameter = 0
fitted_parameters = [0] * 4

pi_pulse_amplitudes[qubit] = pi_pulse_parameter
fitted_parameters[qubit] = translated_popt

return RabiAmplitudeResults(pi_pulse_amplitudes, data.durations, fitted_parameters)
popt = [0] * 4
perr = [1] * 4

pi_pulse_amplitudes[qubit] = (pi_pulse_parameter, perr[2] / 2)
fitted_parameters[qubit] = popt.tolist()
durations = {key: (value, 0) for key, value in data.durations.items()}
chi2[qubit] = (
chi2_reduced(
y,
utils.rabi_amplitude_fit(x, *popt),
qubit_data.error,
),
np.sqrt(2 / len(y)),
)
return RabiAmplitudeResults(pi_pulse_amplitudes, durations, fitted_parameters, chi2)


def _plot(data: RabiAmplitudeData, qubit, fit: RabiAmplitudeResults = None):
"""Plotting function for RabiAmplitude."""
return utils.plot(data, qubit, fit)
return utils.plot_probabilities(data, qubit, fit)


def _update(results: RabiAmplitudeResults, platform: Platform, qubit: QubitId):
Expand Down
183 changes: 183 additions & 0 deletions src/qibocal/protocols/characterization/rabi/amplitude_msr.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,183 @@
from dataclasses import dataclass

import numpy as np
from qibolab import AcquisitionType, AveragingMode, ExecutionParameters
from qibolab.platform import Platform
from qibolab.pulses import PulseSequence
from qibolab.qubits import QubitId
from qibolab.sweeper import Parameter, Sweeper, SweeperType
from scipy.optimize import curve_fit
from scipy.signal import find_peaks

from qibocal import update
from qibocal.auto.operation import Qubits, Routine
from qibocal.config import log
from qibocal.protocols.characterization.rabi.amplitude import (
RabiAmplitudeData,
RabiAmplitudeParameters,
RabiAmplitudeResults,
)

from . import utils


@dataclass
class RabiAmplitudeVoltParameters(RabiAmplitudeParameters):
"""RabiAmplitude runcard inputs."""


@dataclass
class RabiAmplitudeVoltResults(RabiAmplitudeResults):
"""RabiAmplitude outputs."""


RabiAmpVoltType = np.dtype(
[("amp", np.float64), ("msr", np.float64), ("phase", np.float64)]
)
"""Custom dtype for rabi amplitude."""


@dataclass
class RabiAmplitudeVoltData(RabiAmplitudeData):
"""RabiAmplitude data acquisition."""


def _acquisition(
params: RabiAmplitudeVoltParameters, platform: Platform, qubits: Qubits
) -> RabiAmplitudeVoltData:
r"""
Data acquisition for Rabi experiment sweeping amplitude.
In the Rabi experiment we apply a pulse at the frequency of the qubit and scan the drive pulse amplitude
to find the drive pulse amplitude that creates a rotation of a desired angle.
"""

# create a sequence of pulses for the experiment
sequence = PulseSequence()
qd_pulses = {}
ro_pulses = {}
durations = {}
for qubit in qubits:
qd_pulses[qubit] = platform.create_RX_pulse(qubit, start=0)
if params.pulse_length is not None:
qd_pulses[qubit].duration = params.pulse_length

durations[qubit] = qd_pulses[qubit].duration
ro_pulses[qubit] = platform.create_qubit_readout_pulse(
qubit, start=qd_pulses[qubit].finish
)
sequence.add(qd_pulses[qubit])
sequence.add(ro_pulses[qubit])

# define the parameter to sweep and its range:
# qubit drive pulse amplitude
qd_pulse_amplitude_range = np.arange(
params.min_amp_factor,
params.max_amp_factor,
params.step_amp_factor,
)
sweeper = Sweeper(
Parameter.amplitude,
qd_pulse_amplitude_range,
[qd_pulses[qubit] for qubit in qubits],
type=SweeperType.FACTOR,
)

data = RabiAmplitudeVoltData(durations=durations)

# sweep the parameter
results = platform.sweep(
sequence,
ExecutionParameters(
nshots=params.nshots,
relaxation_time=params.relaxation_time,
acquisition_type=AcquisitionType.INTEGRATION,
averaging_mode=AveragingMode.CYCLIC,
),
sweeper,
)
for qubit in qubits:
result = results[ro_pulses[qubit].serial]
data.register_qubit(
RabiAmpVoltType,
(qubit),
dict(
amp=qd_pulses[qubit].amplitude * qd_pulse_amplitude_range,
msr=result.magnitude,
phase=result.phase,
),
)
return data


def _fit(data: RabiAmplitudeVoltData) -> RabiAmplitudeVoltResults:
"""Post-processing for RabiAmplitude."""
qubits = data.qubits

pi_pulse_amplitudes = {}
fitted_parameters = {}

for qubit in qubits:
qubit_data = data[qubit]

rabi_parameter = qubit_data.amp
voltages = qubit_data.msr

y_min = np.min(voltages)
y_max = np.max(voltages)
x_min = np.min(rabi_parameter)
x_max = np.max(rabi_parameter)
x = (rabi_parameter - x_min) / (x_max - x_min)
y = (voltages - y_min) / (y_max - y_min)

# Guessing period using fourier transform
ft = np.fft.rfft(y)
mags = abs(ft)
local_maxima = find_peaks(mags, threshold=10)[0]
index = local_maxima[0] if len(local_maxima) > 0 else None
# 0.5 hardcoded guess for less than one oscillation
f = x[index] / (x[1] - x[0]) if index is not None else 0.5
pguess = [0.5, 1, 1 / f, np.pi / 2]
try:
popt, _ = curve_fit(
utils.rabi_amplitude_fit,
x,
y,
p0=pguess,
maxfev=100000,
bounds=(
[0, 0, 0, -np.pi],
[1, 1, np.inf, np.pi],
),
)
translated_popt = [ # Change it according to fit function changes
y_min + (y_max - y_min) * popt[0],
(y_max - y_min) * popt[1],
popt[2] * (x_max - x_min),
popt[3] - 2 * np.pi * x_min / (x_max - x_min) / popt[2],
]
pi_pulse_parameter = np.abs((translated_popt[2]) / 2)

except:
log.warning("rabi_fit: the fitting was not succesful")
pi_pulse_parameter = 0
fitted_parameters = [0] * 4

pi_pulse_amplitudes[qubit] = pi_pulse_parameter
fitted_parameters[qubit] = translated_popt

return RabiAmplitudeVoltResults(
pi_pulse_amplitudes, data.durations, fitted_parameters
)


def _plot(data: RabiAmplitudeVoltData, qubit, fit: RabiAmplitudeVoltResults = None):
"""Plotting function for RabiAmplitude."""
return utils.plot(data, qubit, fit)


def _update(results: RabiAmplitudeVoltResults, platform: Platform, qubit: QubitId):
update.drive_amplitude(results.amplitude[qubit], platform, qubit)


rabi_amplitude_msr = Routine(_acquisition, _fit, _plot, _update)
"""RabiAmplitude Routine object."""
Loading

0 comments on commit 6fa349d

Please sign in to comment.