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
[docs]
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
[docs]
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
[docs]
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)
[docs]
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)
[docs]
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
[docs]
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 <https://arxiv.org/abs/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
[docs]
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)
[docs]
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
[docs]
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
[docs]
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"]
)
[docs]
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
[docs]
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
[docs]
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
[docs]
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
[docs]
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,
)
[docs]
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
[docs]
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,
)