2. filter_functions package

2.1. Submodules

2.2. filter_functions.analytic module

This file provides functions for the analytical solutions to some of the dynamical decoupling sequences. Note that the filter functions given here differ by a factor of 1/omega**2 from those defined in this package due to different conventions. See for example [Cyw08]. Depending on the definition of the noise Hamiltonian one might also get different results. The functions here agree for

\[B_\alpha\equiv\sigma_z/2.\]

2.2.1. Functions

FID(): Free Induction Decay / Ramsey pulse
SE(): Spin Echo
PDD(): Periodic Dynamical Decoupling
CPMG(): Carr-Purcell-Meiboom-Gill Sequence
CDD(): Concatenated Dynamical Decoupling
UDD(): Uhrig Dynamical Decoupling

References

Cyw08: Cywiński, Ł., Lutchyn, R. M., Nave, C. P., & Das Sarma, S. (2008). How to enhance dephasing time in superconducting qubits. Physical Review B - Condensed Matter and Materials Physics, 77(17), 1–11. https://doi.org/10.1103/PhysRevB.77.174509

CDD(z, g)[source]

CPMG(z, n)[source]

FID(z)[source]

PDD(z, n)[source]

SE(z)[source]

UDD(z, n)[source]

2.3. filter_functions.basis module

This module defines the Basis class, a subclass of NumPy’s ndarray, to represent operator bases.

2.3.1. Classes

Basis: The operator basis as an array of shape (d**2, d, d) with d the dimension of the Hilbert space

2.3.2. Functions

normalize(): Function to normalize a Basis instance
expand(): Function to expand an array of operators in a given basis
ggm_expand(): Fast function to expand an array of operators in a Generalized Gell-Mann basis

class Basis(basis_array: Sequence, traceless: Optional[bool] = None, btype: Optional[str] = None, labels: Optional[Sequence[str]] = None)[source]

Bases: ndarray

Class for operator bases. There are several ways to instantiate a Basis object:

by just calling this constructor with a (possibly incomplete) array of basis matrices. No checks regarding orthonormality or hermiticity are performed.

by calling one of the classes alternative constructors (classmethods):

pauli(): Pauli operator basis

ggm(): Generalized Gell-Mann basis

from_partial() Generate an complete basis from partial elements

These bases guarantee the following properties:

hermitian

orthonormal

[traceless] (can be controlled by a flag)

Since Basis is a subclass of NumPy’s ndarray, it inherits all of its attributes, e.g. shape. The following attributes behave slightly differently to a ndarray, however

A == B is True if all elements evaluate almost equal, i.e. equivalent to np.allclose(A, B).

basis.T transposes the last two axes of basis. For a full basis, this corresponds to transposing each element individually. For a basis element, it corresponds to normal transposition.

Parameters

basis_array: array_like, shape (n, d, d): An array or list of square matrices that are elements of an operator basis spanning \(\mathbb{C}^{d\times d}\). n should be smaller than or equal to d**2.
traceless: bool, optional (default: auto): Controls whether a traceless basis is forced. Here, traceless means that the first element of the basis is the identity and the remaining elements are matrices of trace zero. If an element of basis_array is neither traceless nor the identity and traceless == True, an exception will be raised. Defaults to True if basis_array is traceless and False if not.
btype: str, optional (default: ``’custom’``): A string describing the basis type. For example, a basis created by the factory method pauli() has btype ‘pauli’.
labels: sequence of str, optional: A list of labels for the individual basis elements. Defaults to ‘C_0’, ‘C_1’, …

Attributes

Other than the attributes inherited from ``ndarray``, a ``Basis``
instance has the following attributes:
btype: str: Basis type.
labels: sequence of str: The labels for the basis elements.
d: int: Dimension of the space spanned by the basis.
H: Basis: Hermitian conjugate.
isherm: bool: If the basis is hermitian.
isorthonorm: bool: If the basis is orthonormal.
istraceless: bool: If the basis is traceless except for an identity element
iscomplete: bool: If the basis is complete, ie spans the full space.
sparse: COO, shape (n, d, d): Representation in the COO format supplied by the sparse package.
four_element_traces: COO, shape (n, n, n, n): Traces over all possible combinations of four elements of self. This is required for the calculation of the error transfer matrix and thus cached in the Basis instance.
Most of the attributes above are properties which are lazily
evaluated and cached.

Methods

Other than the methods inherited from ``ndarray``, a ``Basis``
instance has the following methods:
normalize(b)	Normalizes the basis (used internally when creating a basis from elements)
tidyup(eps_scale=None)	Cleans up floating point errors in-place to make zeros actual zeros. `eps_scale` is an optional argument multiplied to the data type’s `eps` to get the absolute tolerance.

Constructor.

property H: Basis: Return the basis hermitian conjugated element-wise.

property T: Basis: Return the basis transposed element-wise.

property four_element_traces: COO: Return all traces of the form \(\mathrm{tr}(C_i C_j C_k C_l)\) as a sparse COO array for \(i,j,k,l > 0\) (i.e. excluding the identity).

classmethod from_partial(partial_basis_array: Sequence, traceless: Optional[bool] = None, btype: Optional[str] = None, labels: Optional[Sequence[str]] = None) → Basis[source]

Generate complete and orthonormal basis from a partial set.

The basis is completed using singular value decomposition to determine the null space of the expansion coefficients of the partial basis with respect to another complete basis.

Parameters

partial_basis_array: array_like: A sequence of basis elements.
traceless: bool, optional: If a traceless basis should be generated (i.e. the first element is the identity and all the others have trace zero).
btype: str, optional: A custom identifier.
labels: Sequence[str], optional: A list of custom labels for each element. If len(labels) == len(partial_basis_array), the newly created elements get labels ‘C_i’.

Returns

basis: Basis, shape (d**2, d, d): The orthonormal basis.

Raises

ValueError: If the given elements are not orthonormal.
ValueError: If the given elements are not traceless but traceless==True.
ValueError: If not len(partial_basis_array) or d**2 labels were given.

classmethod ggm(d: int) → Basis[source]

Returns a generalized Gell-Mann basis in \(d\) dimensions [Bert08] where the elements \(\Lambda_i\) are normalized with respect to the Hilbert-Schmidt inner product,

\[\begin{split}\langle\Lambda_i,\Lambda_j\rangle &= \mathrm{Tr}\,\Lambda_i^\dagger\Lambda_j \\ &= \delta_{ij}.\end{split}\]

Parameters

d: int: The dimensionality of the space spanned by the basis

Returns

basis: Basis: The Basis object representing the GGM.

References

Bert08: Bertlmann, R. A., & Krammer, P. (2008). Bloch vectors for qudits. Journal of Physics A: Mathematical and Theoretical, 41(23). https://doi.org/10.1088/1751-8113/41/23/235303

property iscomplete: bool: Returns True if basis is complete.

property isherm: bool: Returns True if all basis elements are hermitian.

property isorthonorm: bool: Returns True if basis is orthonormal.

property istraceless: bool: Returns True if basis is traceless except for possibly the identity.

normalize(copy: bool = False) → Union[None, Basis][source]: Normalize the basis.

classmethod pauli(n: int) → Basis[source]

Returns a Pauli basis for \(n\) qubits, i.e. the basis spans the space \(\mathbb{C}^{d\times d}\) with \(d = 2^n\):

\[\mathcal{P} = \{I, X, Y, Z\}^{\otimes n}.\]

The elements \(\sigma_i\) are normalized with respect to the Hilbert-Schmidt inner product,

\[\begin{split}\langle\sigma_i,\sigma_j\rangle &= \mathrm{Tr}\,\sigma_i^\dagger\sigma_j \\ &= \delta_{ij}.\end{split}\]

Parameters

n: int: The number of qubits.

Returns

basis: Basis: The Basis object representing the Pauli basis.

property sparse: COO: Return the basis as a sparse COO array

tidyup(eps_scale: Optional[float] = None) → None[source]: Wraps util.remove_float_errors.

expand(M: Union[ndarray, Basis], basis: Union[ndarray, Basis], normalized: bool = True, hermitian: bool = False, tidyup: bool = False) → ndarray[source]

Expand the array M in the basis given by basis.

Parameters

M: array_like: The square matrix (d, d) or array of square matrices (…, d, d) to be expanded in basis
basis: array_like: The basis of shape (m, d, d) in which to expand.
normalized: bool {True}: Wether the basis is normalized.
hermitian: bool (default: False): If M is hermitian along its last two axes, the result will be real.
tidyup: bool {False}: Whether to set values below the floating point eps to zero.

Returns

coefficients: ndarray: The coefficient array with shape (…, m) or (m,) if M was 2-d

Notes

For an orthogonal matrix basis \(\mathcal{C} = \big\{C_k\in \mathbb{C}^{d\times d}: \langle C_k,C_l\rangle_\mathrm{HS} \propto \delta_{kl}\big\}_{k=0}^{d^2-1}\) with the Hilbert-Schmidt inner product as implemented by dot_HS() and \(M\in\mathbb{C}^{d\times d}\), the expansion of \(M\) in terms of \(\mathcal{C}\) is given by

\[\begin{split}M &= \sum_j c_j C_j, \\ c_j &= \frac{\mathrm{tr}\big(M C_j\big)} {\mathrm{tr}\big(C_j^\dagger C_j\big)}.\end{split}\]

ggm_expand(M: Union[ndarray, Basis], traceless: bool = False, hermitian: bool = False) → ndarray[source]

Expand the matrix M in a Generalized Gell-Mann basis [Bert08]. This function makes use of the explicit construction prescription of the basis and thus makes do without computing the expansion coefficients as the overlap between the matrix and each basis element.

Parameters

M: array_like: The square matrix (d, d) or array of square matrices (…, d, d) to be expanded in a GGM basis.
traceless: bool (default: False): Include the basis element proportional to the identity in the expansion. If it is known beforehand that M is traceless, the corresponding coefficient is zero and thus doesn’t need to be computed.
hermitian: bool (default: False): If M is hermitian along its last two axes, the result will be real.

Returns

coefficients: ndarray: The real coefficient array with shape (d**2,) or (…, d**2)

References

Bert08: Bertlmann, R. A., & Krammer, P. (2008). Bloch vectors for qudits. Journal of Physics A: Mathematical and Theoretical, 41(23). https://doi.org/10.1088/1751-8113/41/23/235303

normalize(b: Basis) → Basis[source]

Return a copy of the basis b normalized with respect to the Frobenius norm [Gol85]:

\(||A||_F = \left[\sum_{i,j} |a_{i,j}|^2\right]^{1/2}\)

or equivalently, with respect to the Hilbert-Schmidt inner product as implemented by dot_HS().

References

Gol85: G. H. Golub and C. F. Van Loan, Matrix Computations, Baltimore, MD, Johns Hopkins University Press, 1985, pg. 15

2.4. filter_functions.numeric module

This module defines the functions to calculate everything related to filter functions.

2.4.1. Functions

calculate_control_matrix_from_atomic(): Calculate the control matrix from those of atomic pulse sequences
calculate_control_matrix_from_scratch(): Calculate the control matrix from scratch
calculate_control_matrix_periodic(): Calculate the control matrix for a periodic Hamiltonian
calculate_noise_operators_from_atomic(): Calculate the interaction picture noise operators from atomic segments. Same calculation as calculate_control_matrix_from_atomic() except in Hilbert space.
calculate_noise_operators_from_scratch(): Calculate the interaction picture noise operators from scratch. Same calculation as calculate_control_matrix_from_scratch() except in Hilbert space.
calculate_cumulant_function(): Calculate the cumulant function for a given PulseSequence object.
calculate_decay_amplitudes(): Calculate the decay amplitudes, corresponding to first order terms of the Magnus expansion
calculate_frequency_shifts(): Calculate the frequency shifts, corresponding to second order terms of the Magnus expansion
calculate_filter_function(): Calculate the filter function from the control matrix
calculate_second_order_filter_function(): Calculate the second order filter function used to compute the frequency shifts.
calculate_pulse_correlation_filter_function(): Calculate the pulse correlation filter function from the control matrix
diagonalize(): Diagonalize a Hamiltonian
error_transfer_matrix(): Calculate the error transfer matrix of a pulse up to a unitary rotation
infidelity(): Function to compute the infidelity of a pulse defined by a PulseSequence instance for a given noise spectral density and frequencies

calculate_control_matrix_from_atomic(phases: ndarray, control_matrix_atomic: ndarray, propagators_liouville: ndarray, show_progressbar: bool = False, which: str = 'total') → ndarray[source]

Calculate the control matrix from the control matrices of atomic segments.

Parameters

phases: array_like, shape (n_dt, n_omega): The phase factors for \(g\in\{0, 1, \dots, G-1\}\).
control_matrix_atomic: array_like, shape (n_dt, n_nops, d**2, n_omega): The pulse control matrices for \(g\in\{1, 2, \dots, G\}\).
propagators_liouville: array_like, shape (n_dt, n_nops, d**2, d**2): The transfer matrices of the cumulative propagators for \(g\in\{0, 1, \dots, G-1\}\).
show_progressbar: bool, optional: Show a progress bar for the calculation.
which: str, (‘total’, ‘correlations’): Compute the total control matrix (the sum of all time steps) or the correlation control matrix (first axis holds each pulses’ contribution)

Returns

control_matrix: ndarray, shape ([n_pls,] n_nops, d**2, n_omega): The control matrix \(\tilde{\mathcal{B}}(\omega)\).

See also

calculate_control_matrix_from_scratch: Control matrix from scratch.
liouville_representation: Liouville representation for a given basis.

Notes

The control matrix is calculated by evaluating the sum

\[\tilde{\mathcal{B}}(\omega) = \sum_{g=1}^G e^{i\omega t_{g-1}} \tilde{\mathcal{B}}^{(g)}(\omega)\mathcal{Q}^{(g-1)}.\]

calculate_control_matrix_from_scratch(eigvals: ndarray, eigvecs: ndarray, propagators: ndarray, omega: Sequence[float], basis: Basis, n_opers: Sequence[Union[ndarray, Qobj]], n_coeffs: Sequence[Sequence[float]], dt: Sequence[float], t: Optional[Sequence[float]] = None, show_progressbar: bool = False, cache_intermediates: bool = False, out: Optional[ndarray] = None) → Union[ndarray, Tuple[ndarray, Dict[str, ndarray]]][source]

Calculate the control matrix from scratch, i.e. without knowledge of the control matrices of more atomic pulse sequences.

Parameters

eigvals: array_like, shape (n_dt, d): Eigenvalue vectors for each time pulse segment g with the first axis counting the pulse segment, i.e. eigvals == array([D_0, D_1, ...]).
eigvecs: array_like, shape (n_dt, d, d): Eigenvector matrices for each time pulse segment g with the first axis counting the pulse segment, i.e. eigvecs == array([V_0, V_1, ...]).
propagators: array_like, shape (n_dt+1, d, d): The propagators \(Q_g = P_g P_{g-1}\cdots P_0\) as a (d, d) array with d the dimension of the Hilbert space.
omega: array_like, shape (n_omega,): Frequencies at which the pulse control matrix is to be evaluated.
basis: Basis, shape (d**2, d, d): The basis elements in which the pulse control matrix will be expanded.
n_opers: array_like, shape (n_nops, d, d): Noise operators \(B_\alpha\).
n_coeffs: array_like, shape (n_nops, n_dt): The sensitivities of the system to the noise operators given by n_opers at the given time step.
dt: array_like, shape (n_dt): Sequence duration, i.e. for the \(g\)-th pulse \(t_g - t_{g-1}\).
t: array_like, shape (n_dt+1), optional: The absolute times of the different segments. Can also be computed from dt.
show_progressbar: bool, optional: Show a progress bar for the calculation.
cache_intermediates: bool, optional: Keep and return intermediate terms of the calculation that can be reused in other computations (second order or gradients). Otherwise the sum is performed in-place. The default is False.
out: ndarray, optional: A location into which the result is stored. See numpy.ufunc().

Returns

control_matrix: ndarray, shape (n_nops, d**2, n_omega): The control matrix \(\tilde{\mathcal{B}}(\omega)\)
intermediates: dict[str, ndarray]: Intermediate results of the calculation. Only if cache_intermediates is True.

See also

calculate_control_matrix_from_atomic: Control matrix from concatenation.
calculate_control_matrix_periodic: Control matrix for periodic system.

Notes

The control matrix is calculated according to

\[\tilde{\mathcal{B}}_{\alpha k}(\omega) = \sum_{g=1}^G e^{i\omega t_{g-1}} s_\alpha^{(g)}\mathrm{tr}\left( [\bar{B}_\alpha^{(g)}\circ I(\omega)] \bar{C}_k^{(g)} \right)\]

where

\[\begin{split}I^{(g)}_{nm}(\omega) &= \int_0^{t_l - t_{g-1}}\mathrm{d}t\, e^{i(\omega+\omega_n-\omega_m)t} \\ &= \frac{e^{i(\omega+\omega_n-\omega_m) (t_l - t_{g-1})} - 1} {i(\omega+\omega_n-\omega_m)}, \\ \bar{B}_\alpha^{(g)} &= V^{(g)\dagger} B_\alpha V^{(g)}, \\ \bar{C}_k^{(g)} &= V^{(g)\dagger} Q_{g-1} C_k Q_{g-1}^\dagger V^{(g)},\end{split}\]

and \(V^{(g)}\) is the matrix of eigenvectors that diagonalizes \(\tilde{\mathcal{H}}_n^{(g)}\), \(B_\alpha\) the \(\alpha\)-th noise operator \(s_\alpha^{(g)}\) the noise sensitivity during interval \(g\), and \(C_k\) the \(k\)-th basis element.

calculate_control_matrix_periodic(phases: ndarray, control_matrix: ndarray, total_propagator_liouville: ndarray, repeats: int, check_invertible: bool = True) → ndarray[source]

Calculate the control matrix of a periodic pulse given the phase factors, control matrix and transfer matrix of the total propagator, total_propagator_liouville, of the atomic pulse.

Parameters

phases: ndarray, shape (n_omega,): The phase factors \(e^{i\omega T}\) of the atomic pulse.
control_matrix: ndarray, shape (n_nops, d**2, n_omega): The control matrix \(\tilde{\mathcal{B}}^{(1)}(\omega)\) of the atomic pulse.
total_propagator_liouville: ndarray, shape (d**2, d**2): The transfer matrix \(\mathcal{Q}^{(1)}\) of the atomic pulse.
repeats: int: The number of repetitions.
check_invertible: bool, optional: Check for all frequencies if the inverse \(\mathbb{I} - e^{i\omega T} \mathcal{Q}^{(1)}\) exists by calculating the determinant. If it does not exist, the sum is evaluated explicitly for those points. The default is True.

Returns

control_matrix: ndarray, shape (n_nops, d**2, n_omega): The control matrix \(\tilde{\mathcal{B}}(\omega)\) of the repeated pulse.

Notes

The control matrix is computed as

\[\begin{split}\tilde{\mathcal{B}}(\omega) &= \tilde{\mathcal{B}}^{(1)}(\omega) \sum_{g=0}^{G-1} \left(e^{i\omega T}\right)^g \\ &= \tilde{\mathcal{B}}^{(1)}(\omega)\bigl( \mathbb{I} - e^{i\omega T} \mathcal{Q}^{(1)}\bigr)^{-1}\bigl( \mathbb{I} - \bigl(e^{i\omega T} \mathcal{Q}^{(1)}\bigr)^G\bigr).\end{split}\]

with \(G\) the number of repetitions.

calculate_cumulant_function(pulse: PulseSequence, spectrum: Optional[ndarray] = None, omega: Optional[Sequence[float]] = None, n_oper_identifiers: Optional[Sequence[str]] = None, which: str = 'total', second_order: bool = False, decay_amplitudes: Optional[ndarray] = None, frequency_shifts: Optional[ndarray] = None, show_progressbar: bool = False, memory_parsimonious: bool = False, cache_intermediates: Optional[bool] = None) → ndarray[source]

Calculate the cumulant function \(\mathcal{K}(\tau)\).

The error transfer matrix is obtained from the cumulant function by exponentiation, \(\langle\tilde{\mathcal{U}}\rangle = \exp\mathcal{K}(\tau)\).

Parameters

pulse: PulseSequence: The PulseSequence instance for which to compute the cumulant function.
spectrum: array_like, shape ([[n_nops,] n_nops,] n_omega), optional: The noise power spectral density in units of inverse frequencies as an array of shape (n_omega,), (n_nops, n_omega), or (n_nops, n_nops, n_omega). In the first case, the same spectrum is taken for all noise operators, in the second, it is assumed that there are no correlations between different noise sources and thus there is one spectrum for each noise operator. In the third and most general case, there may be a spectrum for each pair of noise operators corresponding to the correlations between them. n_nops is the number of noise operators considered and should be equal to len(n_oper_identifiers).
omega: array_like, shape (n_omega,), optional: The frequencies at which to evaluate the filter functions.
n_oper_identifiers: array_like, optional: The identifiers of the noise operators for which to evaluate the cumulant function. The default is all.
which: str, optional: Which decay amplitudes should be calculated, may be either ‘total’ (default) or ‘correlations’. See infidelity() and Notes. Note that the latter is not available for the second order terms.
second_order: bool, optional: Also take into account the frequency shifts \(\Delta\) that correspond to second order Magnus expansion and constitute unitary terms. Default False.
decay_amplitudes, array_like, shape ([[n_pls, n_pls,] n_nops,] n_nops, d**2, d**2), optional: A precomputed cumulant function. If given, spectrum, omega are not required.
frequency_shifts, array_like, shape ([[n_pls, n_pls,] n_nops,] n_nops, d**2, d**2), optional: A precomputed frequency shift. If given, spectrum, omega are not required for second order terms.
show_progressbar: bool, optional: Show a progress bar for the calculation of the control matrix.
memory_parsimonious: bool, optional: Trade memory footprint for performance. See calculate_decay_amplitudes(). The default is False.
cache_intermediates: bool, optional: Keep and return intermediate terms of the calculation of the control matrix that can be reused in other computations (second order or gradients). Otherwise the sum is performed in-place. Default is True if second_order=True, else False.

Returns

cumulant_function: ndarray, shape ([[n_pls, n_pls,] n_nops,] n_nops, d**2, d**2): The cumulant function. The individual noise operator contributions chosen by n_oper_identifiers are on the third to last axis / axes, depending on whether the noise is cross-correlated or not. If which == 'correlations', the first two axes correspond to the contributions of the pulses in the sequence.

See also

calculate_decay_amplitudes: Calculate the \(\Gamma_{\alpha\beta,kl}\)
error_transfer_matrix: Calculate the error transfer matrix \(\exp\mathcal{K}\).
infidelity: Calculate only infidelity of a pulse.
pulse_sequence.concatenate: Concatenate PulseSequence objects.
calculate_pulse_correlation_filter_function

Notes

The cumulant function is given by

\[\begin{split}K_{\alpha\beta,ij}(\tau) = -\frac{1}{2} \sum_{kl}\biggl( &\Delta_{\alpha\beta,kl}\left( T_{klji} - T_{lkji} - T_{klij} + T_{lkij} \right) \\ + &\Gamma_{\alpha\beta,kl}\left( T_{klji} - T_{kjli} - T_{kilj} + T_{kijl} \right) \biggr)\end{split}\]

Here, \(T_{ijkl} = \mathrm{tr}(C_i C_j C_k C_l)\) is a trivial function of the basis elements \(C_i\), and \(\Gamma_{\alpha\beta,kl}\) and \(\Delta_{\alpha\beta,kl}\) are the decay amplitudes and frequency shifts which correspond to first and second order in the Magnus expansion, respectively. Since the latter induce coherent errors, we can approximately neglect them if we assume that the pulse has been experimentally calibrated.

For a single qubit and represented in the Pauli basis, the above reduces to

\[\begin{split}K_{\alpha\beta,ij}(\tau) = \begin{cases} - \sum_{k\neq i}\Gamma_{\alpha\beta,kk} &\quad\mathrm{if}\: i = j, \\ - \Delta_{\alpha\beta,ij} + \Delta_{\alpha\beta,ji} + \Gamma_{\alpha\beta,ij} &\quad\mathrm{if}\: i\neq j, \end{cases}\end{split}\]

for \(i\in\{1,2,3\}\) and \(K_{0j} = K_{i0} = 0\).

Lastly, the pulse correlation cumulant function resolves correlations in the cumulant function of a sequence of pulses \(g = 1,\dotsc,G\) such that the following holds:

\[K_{\alpha\beta,ij}(\tau) = \sum_{g,g'=1}^G K_{\alpha\beta,ij}^{(gg')}(\tau).\]

calculate_decay_amplitudes(pulse: PulseSequence, spectrum: ndarray, omega: Sequence[float], n_oper_identifiers: Optional[Sequence[str]] = None, which: str = 'total', show_progressbar: bool = False, cache_intermediates: bool = False, memory_parsimonious: bool = False) → ndarray[source]

Get the decay amplitudes \(\Gamma_{\alpha\beta, kl}\) for noise sources \(\alpha,\beta\) and basis elements \(k,l\).

Parameters

pulse: PulseSequence

The PulseSequence instance for which to compute the decay amplitudes.

spectrum: array_like, shape ([[n_nops,] n_nops,] n_omega)

The noise power spectral density. If 1-d, the same spectrum is used for all noise operators. If 2-d, one (self-) spectrum for each noise operator is expected. If 3-d, should be a matrix of cross-spectra such that spectrum[i, j] == spectrum[j, i].conj().

omega: array_like,

The frequencies at which to calculate the filter functions.

n_oper_identifiers: array_like, optional

The identifiers of the noise operators for which to calculate the decay amplitudes. The default is all.

which: str, optional

Which decay amplitudes should be calculated, may be either ‘total’ (default) or ‘correlations’. See infidelity() and Notes.

show_progressbar: bool, optional

Show a progress bar for the calculation.

cache_intermediates: bool, optional

Keep and return intermediate terms of the calculation that are useful in other places (if control matrix not already cached).

memory_parsimonious: bool, optional

For large dimensions, the integrand

\[\tilde{\mathcal{B}}^\ast_{\alpha k}(\omega) S_{\alpha\beta}(\omega)\tilde{\mathcal{B}}_{\beta l}(\omega)\]

can consume quite a large amount of memory if set up for all \(\alpha,\beta,k,l\) at once. If True, it is only set up and integrated for a single \(k\) at a time and looped over. This is slower but requires much less memory. The default is False.

Returns

decay_amplitudes: ndarray, shape ([[n_pls, n_pls,] n_nops,] n_nops, d**2, d**2): The decay amplitudes.

Raises

ValueError: If spectrum has incompatible shape.

See also

infidelity: Compute the infidelity directly.
pulse_sequence.concatenate: Concatenate PulseSequence objects.
calculate_frequency_shifts: Second order (unitary) terms.
calculate_pulse_correlation_filter_function

Notes

The total decay amplitudes are given by

\[\Gamma_{\alpha\beta, kl} = \int\frac{\mathrm{d}\omega}{2\pi} \tilde{\mathcal{B}}^\ast_{\alpha k}(\omega) S_{\alpha\beta}(\omega)\tilde{\mathcal{B}}_{\beta l}(\omega).\]

If pulse correlations are taken into account, they are given by

\[\Gamma_{\alpha\beta, kl}^{(gg')} = \int \frac{\mathrm{d}\omega}{2\pi} S_{\alpha\beta}(\omega) F_{\alpha\beta, kl}^{(gg')}(\omega).\]

calculate_filter_function(control_matrix: ndarray, which: str = 'fidelity') → ndarray[source]

Compute the filter function from the control matrix.

Parameters

control_matrix: array_like, shape (n_nops, d**2, n_omega): The control matrix.
whichstr, optional: Which filter function to return. Either ‘fidelity’ (default) or ‘generalized’ (see Notes).

Returns

filter_function: ndarray, shape (n_nops, n_nops, [d**2, d**2,] n_omega): The filter functions for each noise operator correlation. The diagonal corresponds to the filter functions for uncorrelated noise sources.

See also

calculate_control_matrix_from_scratch: Control matrix from scratch.
calculate_control_matrix_from_atomic: Control matrix from concatenation.
calculate_pulse_correlation_filter_function: Pulse correlations.

Notes

The generalized filter function is given by

\[F_{\alpha\beta,kl}(\omega) = \tilde{\mathcal{B}}_{\alpha k}^\ast(\omega) \tilde{\mathcal{B}}_{\beta l}(\omega),\]

where \(\alpha,\beta\) are indices counting the noise operators \(B_\alpha\) and \(k,l\) indices counting the basis elements \(C_k\).

The fidelity filter function is obtained by tracing over the basis indices:

\[F_{\alpha\beta}(\omega) = \sum_{k} F_{\alpha\beta,kk}(\omega).\]

calculate_frequency_shifts(pulse: PulseSequence, spectrum: ndarray, omega: Sequence[float], n_oper_identifiers: Optional[Sequence[str]] = None, show_progressbar: bool = False) → ndarray[source]

Get the frequency shifts \(\Delta_{\alpha\beta, kl}\) for noise sources \(\alpha,\beta\) and basis elements \(k,l\).

Parameters

pulse: PulseSequence: The PulseSequence instance for which to compute the frequency shifts.
spectrum: array_like, shape ([[n_nops,] n_nops,] n_omega): The two-sided noise power spectral density. If 1-d, the same spectrum is used for all noise operators. If 2-d, one (self-) spectrum for each noise operator is expected. If 3-d, should be a matrix of cross-spectra such that spectrum[i, j] == spectrum[j, i].conj().
omega: array_like,: The frequencies. Note that the frequencies are assumed to be symmetric about zero.
n_oper_identifiers: array_like, optional: The identifiers of the noise operators for which to calculate the frequency shifts. The default is all.
show_progressbar: bool, optional: Show a progress bar for the calculation.

Returns

Delta: ndarray, shape ([n_nops,] n_nops, d**2, d**2): The frequency shifts.

Raises

ValueError: If spectrum has incompatible shape.

See also

calculate_second_order_filter_function: Corresponding filter function.
calculate_decay_amplitudes: First order (dissipative) terms.
infidelity: Compute the infidelity directly.
pulse_sequence.concatenate: Concatenate PulseSequence objects.
calculate_pulse_correlation_filter_function

Notes

The total frequency shifts are given by

\[\Delta_{\alpha\beta, kl} = \int_{-\infty}^\infty \frac{\mathrm{d}{\omega}}{2\pi} S_{\alpha\beta}(\omega) F_{\alpha\beta,kl}^{(2)}(\omega)\]

with \(F_{\alpha\beta,kl}^{(2)}(\omega)\) the second order filter function.

calculate_noise_operators_from_atomic(phases: ndarray, noise_operators_atomic: ndarray, propagators: ndarray, show_progressbar: bool = False) → ndarray[source]

Calculate the interaction picutre noise operators from atomic segments.

Parameters

phases: array_like, shape (n_dt, n_omega): The phase factors for \(g\in\{0, 1, \dots, G-1\}\).
noise_operators_atomic: array_like, shape (n_dt, n_nops, d, d, n_omega): The noise operators in the interaction picture of the g-th pulse, i.e. for \(g\in\{1, 2, \dots, G\}\).
propagators: array_like, shape (n_dt, d, d): The cumulative propagators of the pulses \(g\in\{0, 1, \dots, G-1\}\).
show_progressbar: bool, optional: Show a progress bar for the calculation.

Returns

noise_operators: ndarray, shape (n_omega, n_nops, d, d): The interaction picture noise operators \(\tilde{B}_\alpha(\omega)\).

See also

calculate_noise_operators_from_scratch: Compute the operators from scratch.
calculate_control_matrix_from_atomic: Same calculation in Liouville space.

Notes

The noise operators are calculated by evaluating the sum

\[\tilde{B}_\alpha(\omega) = \sum_{g=1}^G e^{i\omega t_{g-1}} Q_{g-1}^\dagger\tilde{B}_\alpha^{(g)}(\omega) Q_{g-1}.\]

The control matrix then corresponds to the coefficients of expansion in an operator basis \(\{C_k\}_k\):

\[\tilde{\mathcal{B}}_{k\alpha}(\omega) = \mathrm{tr}(\tilde{B}_\alpha(\omega) C_k).\]

Due to differences in implementation (for performance reasons), the axes of the result are transposed compared to the control matrix:

>>> ctrlmat = calculate_control_matrix_from_atomic(...)
>>> ctrlmat.shape
(n_nops, d**2, n_omega)
>>> noiseops = calculate_noise_operators_from_atomic(...)
>>> noiseops.shape
(n_omega, n_nops, d, d)
>>> ctrlmat_from_noiseops = ff.basis.expand(noiseops, basis)
>>> np.allclose(ctrlmat, ctrlmat_from_noiseops.transpose(1, 2, 0))
True

calculate_noise_operators_from_scratch(eigvals: ndarray, eigvecs: ndarray, propagators: ndarray, omega: Sequence[float], n_opers: Sequence[Union[ndarray, Qobj]], n_coeffs: Sequence[Sequence[float]], dt: Sequence[float], t: Optional[Sequence[float]] = None, show_progressbar: bool = False, cache_intermediates: bool = False) → Union[ndarray, Tuple[ndarray, Dict[str, ndarray]]][source]

Calculate the noise operators in interaction picture from scratch.

Parameters

eigvals: array_like, shape (n_dt, d): Eigenvalue vectors for each time pulse segment g with the first axis counting the pulse segment, i.e. eigvals == array([D_0, D_1, ...]).
eigvecs: array_like, shape (n_dt, d, d): Eigenvector matrices for each time pulse segment g with the first axis counting the pulse segment, i.e. eigvecs == array([V_0, V_1, ...]).
propagators: array_like, shape (n_dt+1, d, d): The propagators \(Q_g = P_g P_{g-1}\cdots P_0\) as a (d, d) array with d the dimension of the Hilbert space.
omega: array_like, shape (n_omega,): Frequencies at which the pulse control matrix is to be evaluated.
n_opers: array_like, shape (n_nops, d, d): Noise operators \(B_\alpha\).
n_coeffs: array_like, shape (n_nops, n_dt): The sensitivities of the system to the noise operators given by n_opers at the given time step.
dt: array_like, shape (n_dt): Sequence duration, i.e. for the \(g\)-th pulse \(t_g - t_{g-1}\).
t: array_like, shape (n_dt+1), optional: The absolute times of the different segments. Can also be computed from dt.
show_progressbar: bool, optional: Show a progress bar for the calculation.
cache_intermediates: bool, optional: Keep and return intermediate terms of the calculation that can be reused in other computations (second order or gradients). Otherwise the sum is performed in-place. The default is False.

Returns

noise_operators: ndarray, shape (n_omega, n_nops, d, d): The interaction picture noise operators \(\tilde{B}_\alpha(\omega)\).
intermediates: dict[str, ndarray]: Intermediate results of the calculation. Only if cache_intermediates is True.

See also

calculate_noise_operators_from_atomic: Compute the operators from atomic segments.
calculate_control_matrix_from_scratch: Same calculation in Liouville space.

Notes

The interaction picture noise operators are calculated according to

\[\tilde{B}_\alpha(\omega) = \sum_{g=1}^G e^{i\omega t_{g-1}} s_\alpha^{(g)} P^{(g)\dagger}\left[ \bar{B}^{(g)}_\alpha \circ I^{(g)}(\omega) \right] P^{(g)}\]

where

\[\begin{split}I^{(g)}_{nm}(\omega) &= \int_0^{t_g - t_{g-1}}\mathrm{d}t\, e^{i(\omega+\omega_n-\omega_m)t} \\ &= \frac{e^{i(\omega+\omega_n-\omega_m) (t_g - t_{g-1})} - 1} {i(\omega+\omega_n-\omega_m)}, \\ \bar{B}_\alpha^{(g)} &= V^{(g)\dagger} B_\alpha V^{(g)}, \\ P^{(g)} &= V^{(g)\dagger} Q_{g-1},\end{split}\]

and \(V^{(g)}\) is the matrix of eigenvectors that diagonalizes \(\tilde{\mathcal{H}}_n^{(g)}\), \(B_\alpha\) the \(\alpha\)-th noise operator, and \(s_\alpha^{(g)}\) the noise sensitivity during interval \(g\).

The control matrix then corresponds to the coefficients of expansion in an operator basis \(\{C_k\}_k\):

\[\tilde{\mathcal{B}}_{k\alpha}(\omega) = \mathrm{tr}(\tilde{B}_\alpha(\omega) C_k).\]

Due to differences in implementation (for performance reasons), the axes of the result are transposed compared to the control matrix:

>>> ctrlmat = calculate_control_matrix_from_scratch(...)
>>> ctrlmat.shape
(n_nops, d**2, n_omega)
>>> noiseops = calculate_noise_operators_from_scratch(...)
>>> noiseops.shape
(n_omega, n_nops, d, d)
>>> ctrlmat_from_noiseops = ff.basis.expand(noiseops, basis)
>>> np.allclose(ctrlmat, ctrlmat_from_noiseops.transpose(1, 2, 0))
True

calculate_pulse_correlation_filter_function(control_matrix: ndarray, which: str = 'fidelity') → ndarray[source]

Compute pulse correlation filter function from control matrix.

Parameters

control_matrix: array_like, shape (n_pulses, n_nops, d**2, n_omega): The control matrix.
whichstr, optional: Which filter function to return. Either ‘fidelity’ (default) or ‘generalized’ (see Notes).

Returns

filter_function_pc: ndarray, shape (n_pls, n_pls, n_nops, n_nops, [d**2, d**2,] n_omega): The pulse correlation filter functions for each pulse and noise operator correlations. The first two axes hold the pulse correlations, the second two the noise correlations.
whichstr, optional: Which filter function to return. Either ‘fidelity’ (default) or ‘generalized’ (see Notes).

See also

calculate_control_matrix_from_scratch: Control matrix from scratch.
calculate_control_matrix_from_atomic: Control matrix from concatenation.
calculate_filter_function: Regular filter function.

Notes

The generalized pulse correlation filter function is given by

\[F_{\alpha\beta,kl}^{(gg')}(\omega) = \bigl[ \mathcal{Q}^{(g'-1)\dagger} \tilde{\mathcal{B}}^{(g')\dagger}(\omega) \bigr]_{k\alpha} \bigl[ \tilde{\mathcal{B}}^{(g)}(\omega)\mathcal{Q}^{(g-1)} \bigr]_{\beta l} e^{i\omega(t_{g-1} - t_{g'-1})},\]

with \(\tilde{\mathcal{B}}^{(g)}\) the control matrix of the \(g\)-th pulse. The fidelity pulse correlation function is obtained by tracing out the basis indices,

\[F_{\alpha\beta}^{(gg')}(\omega) = \sum_{k} F_{\alpha\beta,kk}^{(gg')}(\omega)\]

calculate_second_order_filter_function(eigvals: ndarray, eigvecs: ndarray, propagators: ndarray, omega: Sequence[float], basis: Basis, n_opers: Sequence[Union[ndarray, Qobj]], n_coeffs: Sequence[Sequence[float]], dt: Sequence[float], intermediates: Optional[Dict[str, ndarray]] = None, show_progressbar: bool = False) → ndarray[source]

Calculate the second order filter function for frequency shifts.

Parameters

eigvals: array_like, shape (n_dt, d): Eigenvalue vectors for each time pulse segment l with the first axis counting the pulse segment, i.e. eigvals == array([D_0, D_1, ...]).
eigvecs: array_like, shape (n_dt, d, d): Eigenvector matrices for each time pulse segment l with the first axis counting the pulse segment, i.e. eigvecs == array([V_0, V_1, ...]).
propagators: array_like, shape (n_dt+1, d, d): The propagators \(Q_l = P_l P_{l-1}\cdots P_0\) as a (d, d) array with d the dimension of the Hilbert space.
omega: array_like, shape (n_omega,): Frequencies at which the pulse control matrix is to be evaluated.
basis: Basis, shape (d**2, d, d): The basis elements in which the pulse control matrix will be expanded.
n_opers: array_like, shape (n_nops, d, d): Noise operators \(B_\alpha\).
n_coeffs: array_like, shape (n_nops, n_dt): The sensitivities of the system to the noise operators given by n_opers at the given time step.
dt: array_like, shape (n_dt): Sequence duration, i.e. for the \(l\)-th pulse \(t_l - t_{l-1}\).
intermediates: Dict[str, ndarray], optional: Intermediate terms of the calculation of the control matrix that can be reused here. If None (default), they are computed from scratch.
show_progressbar: bool, optional: Show a progress bar for the calculation.

Returns

second_order_filter_function: ndarray, shape (n_nops, n_nops, d**2, d**2, n_omega): The second order filter function.

See also

calculate_frequency_shifts: Integrate over filter function times spectrum.
calculate_decay_amplitudes: First order (dissipative) terms.
infidelity: Compute the infidelity directly.
pulse_sequence.concatenate: Concatenate PulseSequence objects.
calculate_pulse_correlation_filter_function

Notes

The second order filter function is given by

\[F_{\alpha\beta, kl}^{(2)} = \sum_{g=1}^G\left[ \mathcal{G}_{\alpha k}^{(g)\ast}(\omega) \sum_{g'=1}^{g-1}\mathcal{G}_{\beta l}^{(g')}(\omega) + s_\alpha^{(g)}\bar{B}_{\alpha,ij}^{(g)}\bar{C}_{k,ji}^{(g)} I_{ijmn}^{(g)}(\omega)\bar{C}_{l,nm}^{(g) \bar{B}_{\beta,mn}^{(g)}s_\beta^{(g)}} \right]\]

with

\[\begin{split}\mathcal{G}^{(g)}(\omega) &= e^{i\omega t_{g-1}}\mathcal{B}^{(g)}(\omega) \mathcal{Q}^{(g-1)}, \\ I_{ijmn}^{(g)}(\omega) &= \int_{t_{g-1}}^{t_g}\mathrm{d}{t} e^{i\Omega_{ij}^{(g)}(t - t_{g-1}) - i\omega t} \int_{t_{g-1}}^{t}\mathrm{d}{t'} e^{i\Omega_{mn}^{(g)}(t' - t_{g-1}) + i\omega t'}.\end{split}\]

diagonalize(hamiltonian: ndarray, dt: Sequence[float]) → Tuple[ndarray][source]

Diagonalize a Hamiltonian.

Diagonalize the Hamiltonian which is piecewise constant during the times given by dt and return eigenvalues, eigenvectors, and the cumulative propagators \(Q_l\). Note that we calculate in units where \(\hbar\equiv 1\) so that

\[U(t, t_0) = \mathcal{T}\exp\left( -i\int_{t_0}^t\mathrm{d}t'\mathcal{H}(t') \right).\]

Parameters

hamiltonian: array_like, shape (n_dt, d, d): Hamiltonian of shape (n_dt, d, d) with d the dimensionality of the system
dt: array_like: The time differences

Returns

eigvals: ndarray: Array of eigenvalues of shape (n_dt, d)
eigvecs: ndarray: Array of eigenvectors of shape (n_dt, d, d)
propagators: ndarray: Array of cumulative propagators of shape (n_dt+1, d, d)

error_transfer_matrix(pulse: Optional[PulseSequence] = None, spectrum: Optional[ndarray] = None, omega: Optional[Sequence[float]] = None, n_oper_identifiers: Optional[Sequence[str]] = None, second_order: bool = False, cumulant_function: Optional[ndarray] = None, show_progressbar: bool = False, memory_parsimonious: bool = False, cache_intermediates: bool = False) → ndarray[source]

Compute the error transfer matrix up to unitary rotations.

Parameters

pulse: PulseSequence: The PulseSequence instance for which to compute the error transfer matrix.
spectrum: array_like, shape ([[n_nops,] n_nops,] n_omega): The two-sided noise power spectral density in units of inverse frequencies as an array of shape (n_omega,), (n_nops, n_omega), or (n_nops, n_nops, n_omega). In the first case, the same spectrum is taken for all noise operators, in the second, it is assumed that there are no correlations between different noise sources and thus there is one spectrum for each noise operator. In the third and most general case, there may be a spectrum for each pair of noise operators corresponding to the correlations between them. n_nops is the number of noise operators considered and should be equal to len(n_oper_identifiers).
omega: array_like, shape (n_omega,): The frequencies at which to calculate the filter functions.
n_oper_identifiers: array_like, optional: The identifiers of the noise operators for which to evaluate the error transfer matrix. The default is all. Note that, since in general contributions from different noise operators won’t commute, not selecting all noise operators results in neglecting terms of order \(\xi^4\).
second_order: bool, optional: Also take into account the frequency shifts \(\Delta\) that correspond to second order Magnus expansion and constitute unitary terms. Default False.
cumulant_function: ndarray, shape ([[n_pls, n_pls,] n_nops,] n_nops, d**2, d**2): A precomputed cumulant function. If given, pulse, spectrum, omega are not required.
show_progressbar: bool, optional: Show a progress bar for the calculation of the control matrix.
memory_parsimonious: bool, optional: Trade memory footprint for performance. See calculate_decay_amplitudes(). The default is False.
cache_intermediates: bool, optional: Keep and return intermediate terms of the calculation of the control matrix (if it is not already cached) that can be reused for second order or gradients. Can consume large amount of memory, but speed up the calculation.

Returns

error_transfer_matrix: ndarray, shape (d**2, d**2): The error transfer matrix. The individual noise operator contributions are summed up before exponentiating as they might not commute.

See also

calculate_cumulant_function: Calculate the cumulant function \(\mathcal{K}\)
calculate_decay_amplitudes: Calculate the \(\Gamma_{\alpha\beta,kl}\)
infidelity: Calculate only infidelity of a pulse.

Notes

The error transfer matrix is given by

\[\tilde{\mathcal{U}} = \exp\mathcal{K}(\tau)\]

with \(\mathcal{K}(\tau)\) the cumulant function (see calculate_cumulant_function()). For Gaussian noise this expression is exact when taking into account the decay amplitudes \(\Gamma\) and frequency shifts \(\Delta\). As the latter effects coherent errors it can be neglected if we assume that the experimenter has calibrated their pulse.

For non-Gaussian noise the expression above is perturbative and includes noise up to order \(\xi^2\) and hence \(\tilde{\mathcal{U}} = \mathbb{1} + \mathcal{K}(\tau) + \mathcal{O}(\xi^4)\) (although it is evaluated as a matrix exponential in any case).

Given the above expression of the error transfer matrix, the entanglement fidelity is given by

\[\mathcal{F}_\mathrm{e} = \frac{1}{d^2}\mathrm{tr}\,\tilde{\mathcal{U}}.\]

infidelity(pulse: PulseSequence, spectrum: Union[Sequence[float], Callable], omega: Union[Sequence[float], Dict[str, Union[int, str]]], n_oper_identifiers: Optional[Sequence[str]] = None, which: str = 'total', show_progressbar: bool = False, cache_intermediates: bool = False, return_smallness: bool = False, test_convergence: bool = False) → Union[ndarray, Any][source]

Calculate the leading order entanglement infidelity.

This function calculates the infidelity approximately from the leading peturbation (see Notes). To compute it exactly for Gaussian noise and vanishing coherent errors (second order Magnus terms), use error_transfer_matrix() to obtain it from the full process matrix.

Parameters

pulse: PulseSequence

The PulseSequence instance for which to calculate the infidelity for.

spectrum: array_like, shape ([[n_nops,] n_nops,] omega) or callable

The two-sided noise power spectral density in units of inverse frequencies as an array of shape (n_omega,), (n_nops, n_omega), or (n_nops, n_nops, n_omega). In the first case, the same spectrum is taken for all noise operators, in the second, it is assumed that there are no correlations between different noise sources and thus there is one spectrum for each noise operator. In the third and most general case, there may be a spectrum for each pair of noise operators corresponding to the correlations between them. n_nops is the number of noise operators considered and should be equal to len(n_oper_identifiers).

If test_convergence is True, a function handle to compute the power spectral density from a sequence of frequencies is expected.

omega: array_like or dict

The frequencies at which the integration is to be carried out. If test_convergence is True, a dict with possible keys (‘omega_IR’, ‘omega_UV’, ‘spacing’, ‘n_min’, ‘n_max’, ‘n_points’), where all entries are integers except for spacing which should be a string, either ‘linear’ or ‘log’. ‘n_points’ controls how many steps are taken.

n_oper_identifiers: array_like, optional

The identifiers of the noise operators for which to calculate the infidelity contribution. If given, the infidelities for each noise operator will be returned. Otherwise, all noise operators will be taken into account.

which: str, optional

Which infidelities should be calculated, may be either ‘total’ (default) or ‘correlations’. In the former case, one value is returned for each noise operator, corresponding to the total infidelity of the pulse (or pulse sequence). In the latter, an array of infidelities is returned where element (i,j) corresponds to the infidelity contribution of the correlations between pulses i and j (see Notes). Note that this option is only available if the pulse correlation filter functions have been computed during concatenation (see calculate_pulse_correlation_filter_function() and concatenate()).

show_progressbar: bool, optional

Show a progressbar for the calculation of the control matrix.

cache_intermediates: bool, optional

Keep and return intermediate terms of the calculation of the control matrix (if it is not already cached) that can be reused for second order or gradients. The default is False.

return_smallness: bool, optional

Return the smallness parameter \(\xi\) for the given spectrum.

test_convergence: bool, optional

Test the convergence of the integral with respect to the number of frequency samples. Returns the number of frequency samples and the corresponding fidelities. See spectrum and omega for more information.

Returns

infid: ndarray, shape ([[n_pls, n_pls,], n_nops,] n_nops): Array with the infidelity contributions for each spectrum spectrum on the last axis or axes, depending on the shape of spectrum and which. If which is correlations, the first two axes are the individual pulse contributions. If spectrum is 2-d (3-d), the last axis (two axes) are the individual spectral contributions. Only if test_convergence is False.
n_samples: array_like: Array with number of frequency samples used for convergence test. Only if test_convergence is True.
convergence_infids: array_like: Array with infidelities calculated in convergence test. Only if test_convergence is True.

See also

calculate_decay_amplitudes: Calculate the full matrix of first order terms.
error_transfer_matrix: Calculate the full process matrix.
plotting.plot_infidelity_convergence: Convenience function to plot results.
pulse_sequence.concatenate: Concatenate PulseSequence objects.
calculate_pulse_correlation_filter_function.

Notes

The infidelity is given by

\[\begin{split}\mathcal{I}_{\alpha\beta} &= 1 - \frac{1}{d^2}\mathrm{tr}\:\tilde{\mathcal{U}} \\ &= \frac{1}{d}\int_{-\infty}^{\infty} \frac{\mathrm{d}\omega}{2\pi}S_{\alpha\beta}(\omega) F_{\alpha\beta}(\omega) + \mathcal{O}\big(\xi^4\big) \\ &= \sum_{g,g'=1}^G \mathcal{I}_{\alpha\beta}^{(gg')}\end{split}\]

with \(S_{\alpha\beta}(\omega)\) the two-sided noise spectral density and \(F_{\alpha\beta}(\omega)\) the first-order filter function for noise sources \(\alpha,\beta\). The noise spectrum may include correlated noise sources, that is, its entry at \((\alpha,\beta)\) corresponds to the correlations between sources \(\alpha\) and \(\beta\). \(\mathcal{I}_{\alpha\beta}^{(gg')}\) are the correlation infidelities that can be computed by setting which='correlations'.

To convert to the average gate infidelity, use the following relation given by Horodecki et al. [Hor99] and Nielsen [Nie02]:

\[\mathcal{I}_\mathrm{avg} = \frac{d}{d+1}\mathcal{I}.\]

The smallness parameter is given by

\[\xi^2 = \sum_\alpha\left[ \lvert\lvert B_\alpha\rvert\rvert^2 \int_{-\infty}^\infty\frac{\mathrm{d}\omega}{2\pi} S_\alpha(\omega)\left(\sum_gs_\alpha^{(g)}\Delta t_g \right)^2 \right].\]

Note that in practice, the integral is only evaluated on the interval \(\omega\in[\omega_\mathrm{min},\omega_\mathrm{max}]\).

References

Hor99: Horodecki, M., Horodecki, P., & Horodecki, R. (1999). General teleportation channel, singlet fraction, and quasidistillation. Physical Review A - Atomic, Molecular, and Optical Physics, 60(3), 1888–1898. https://doi.org/10.1103/PhysRevA.60.1888
Nie02: Nielsen, M. A. (2002). A simple formula for the average gate fidelity of a quantum dynamical operation. Physics Letters, Section A: General, Atomic and Solid State Physics, 303(4), 249–252. https://doi.org/10.1016/S0375-9601(02)01272-0

2.5. filter_functions.plotting module

This module provides various plotting functions.

2.5.1. Functions

plot_bloch_vector_evolution(): Plot the evolution of the Bloch vector on a QuTiP-generated Bloch sphere
plot_filter_function(): Plot the filter function of a given PulseSequence
plot_infidelity_convergence(): Helper function called by infidelity() to plot the convergence of the infidelity
plot_pulse_correlation_filter_function(): Plot the pulse correlation filter function of a given PulseSequence
plot_pulse_train(): Plot the pulse train of a given PulseSequence
plot_cumulant_function(): Plot the cumulant function of a PulseSequence for a given spectrum as an image.

plot_bloch_vector_evolution(pulse: PulseSequence, psi0: Optional[Union[ndarray, Qobj]] = None, b: Optional[Bloch] = None, n_samples: Optional[int] = None, cmap: Union[Colormap, str] = 'winter', add_cbar: bool = False, show: bool = True, return_Bloch: bool = False, cbar_kwargs: Optional[dict] = None, **bloch_kwargs) → Union[None, Bloch][source]

Plot the evolution of the Bloch vector under the given pulse sequence.

Parameters

pulse: PulseSequence: The PulseSequence instance whose control Hamiltonian determines the time evolution of the Bloch vector.
psi0: Qobj or array_like, optional: The initial state before the pulse is applied. Defaults to \(|0\rangle\).
b: qutip.Bloch, optional: If given, the QuTiP Bloch instance on which to plot the time evolution.
n_samples: int, optional: The number of time points to be sampled.
cmap: matplotlib colormap, optional: The colormap for the trajectory. Requires matplotlib >= 3.3.0.
add_cbar: bool, optional: Add a colorbar encoding the time evolution to the figure. Default is false.
show: bool, optional: Whether to show the sphere (by calling b.make_sphere()).
return_Bloch: bool, optional: Whether to return the qutip.bloch.Bloch instance
cbar_kwargs: dict, optional: A dictionary with keyword arguments to be fed into the colorbar constructor (if add_cbar == True).
**bloch_kwargs: optional: Keyword arguments to be fed into the qutip.Bloch constructor (if b not given).

Returns

b: qutip.Bloch: The qutip.Bloch instance

Raises

ValueError: If the pulse is for more than one qubit

See also

qutip.bloch.Bloch: Qutip’s Bloch sphere implementation.

plot_cumulant_function(pulse: Optional[PulseSequence] = None, spectrum: Optional[ndarray] = None, omega: Optional[Sequence[float]] = None, cumulant_function: Optional[ndarray] = None, n_oper_identifiers: Optional[Sequence[int]] = None, second_order: bool = False, colorscale: str = 'linear', linthresh: Optional[float] = None, basis_labels: Optional[Sequence[str]] = None, basis_labelsize: Optional[int] = None, cbar_label: str = 'Cumulant Function', cbar_labelsize: Optional[int] = None, fig: Optional[Figure] = None, grid: Optional[ImageGrid] = None, cmap: Optional[Union[Colormap, str]] = None, grid_kw: Optional[dict] = None, cbar_kw: Optional[dict] = None, imshow_kw: Optional[dict] = None, **figure_kw) → Tuple[Figure, ImageGrid][source]

Plot the cumulant function for a given noise spectrum as an image.

The cumulant function generates the error transfer matrix \(\tilde{\mathcal{U}}\) exactly for Gaussian noise and to second order for non-Gaussian noise.

The function may be called with either a PulseSequence, a spectrum, and a list of frequencies in which case the cumulant function is calculated for those parameters, or with a precomputed cumulant function.

As of now, only auto-correlated spectra are implemented.

Parameters

pulse: PulseSequence: The pulse sequence.
spectrum: ndarray: The two-sided noise spectrum.
omega: array_like: The frequencies for which to evaluate the error transfer matrix.
cumulant_function: ndarray, shape (n_nops, d**2, d**2): A precomputed cumulant function. If given, pulse, spectrum, omega are not required.
n_oper_identifiers: array_like, optional: The identifiers of the noise operators for which the cumulant function should be plotted. All identifiers can be accessed via pulse.n_oper_identifiers. Defaults to all.
second_order: bool, optional: Also take into account the frequency shifts \(\Delta\) that correspond to second order Magnus expansion and constitute unitary terms. Default False.
colorscale: str, optional: The scale of the color code (‘linear’ or ‘log’ (default))
linthresh: float, optional: The threshold below which the colorscale will be linear (only for ‘log’) colorscale
basis_labels: array_like (str), optional: Custom labels for the elements of the cumulant function (the basis elements). Grabbed from the basis by default.
basis_labelsize: int, optional: The size in points for the basis labels.
cbar_label: str, optional: The label for the colorbar. Default: ‘Cumulant Function’.
cbar_labelsize: int, optional: The size in points for the colorbar label.
fig: matplotlib figure, optional: A matplotlib figure instance to plot in
grid: matplotlib ImageGrid, optional: An ImageGrid instance to use for plotting.
cmap: matplotlib colormap, optional: The colormap for the matrix plot.
grid_kw: dict, optional: Dictionary with keyword arguments passed to the ImageGrid constructor.
cbar_kw: dict, optional: Dictionary with keyword arguments passed to the colorbar constructor.
imshow_kw: dict, optional: Dictionary with keyword arguments passed to imshow.
figure_kw: optional: Keyword argument dictionaries that are fed into the matplotlib.pyplot.figure() function if no fig instance is specified.

Returns

fig: matplotlib figure: The matplotlib figure instance used for plotting.
grid: matplotlib ImageGrid: The ImageGrid instance used for plotting.

plot_filter_function(pulse: PulseSequence, omega: Optional[Sequence[float]] = None, n_oper_identifiers: Optional[Sequence[int]] = None, fig: Optional[Figure] = None, axes: Optional[Axes] = None, xscale: str = 'log', yscale: str = 'linear', omega_in_units_of_tau: bool = True, cycler: Optional[Cycler] = None, plot_kw: dict = {}, subplot_kw: Optional[dict] = None, gridspec_kw: Optional[dict] = None, **figure_kw) → Tuple[Figure, Axes, Legend][source]

Plot the fidelity filter function(s) of the given PulseSequence for positive frequencies. As of now only the diagonal elements of \(F_{\alpha\beta}\) are implemented, i.e. the filter functions corresponding to uncorrelated noise sources.

Parameters

pulse: PulseSequence: The pulse sequence whose filter function to plot.
omega: array_like, optional: The frequencies at which to evaluate the filter function. If not given, the pulse sequence’s omega attribute is used (if set) or sensible values are chosen automatically (if None)
n_oper_identifiers: array_like, optional: The identifiers of the noise operators for which the filter function should be plotted. All identifiers can be accessed via pulse.n_oper_identifiers. Defaults to all.
fig: matplotlib figure, optional: A matplotlib figure instance to plot in
axes: matplotlib axes, optional: A matplotlib axes instance to use for plotting.
xscale: str, optional: x-axis scaling. One of (‘linear’, ‘log’).
yscale: str, optional: y-axis scaling. One of (‘linear’, ‘log’).
omega_in_units_of_tau: bool, optional: Plot \(\omega\tau\) or just \(\omega\) on x-axis.
cycler: cycler.Cycler, optional: A Cycler instance used to set the style cycle if multiple lines are to be drawn
plot_kw: dict, optional: Dictionary with keyword arguments passed to the plot function
subplot_kw: dict, optional: Dictionary with keyword arguments passed to the subplots constructor
gridspec_kw: dict, optional: Dictionary with keyword arguments passed to the gridspec constructor
figure_kw: optional: Keyword argument dictionaries that are fed into the matplotlib.pyplot.subplots() function if no fig instance is specified.

Returns

fig: matplotlib figure: The matplotlib figure instance used for plotting.
axes: matplotlib axes: The matplotlib axes instance used for plotting.
legend: matplotlib legend: The matplotlib legend instance in the plot.

Raises

ValueError: If an invalid number of n_oper_labels were given

plot_infidelity_convergence(n_samples: Sequence[int], infids: Sequence[float], axes: Optional[Axes] = None) → Tuple[Figure, Axes][source]

Plot the convergence of the infidelity integral. The function arguments are those returned by infidelity() with the test_convergence flag set to True.

Parameters

n_samples: array_like: Array with the number of samples at which the integral was evaluated
infids: array_like, shape (n_samples, [n_oper_inds, optional]): Array with the calculated infidelities for each noise operator on the second axis or the second axis already traced out.
axes: sequence of two matplotlib axes, optional: Two axes that the result is plotted in.

Returns

fig: matplotlib figure: The matplotlib figure instance used for plotting.
axes: matplotlib axes: The matplotlib axes instances used for plotting.

plot_pulse_correlation_filter_function(pulse: PulseSequence, n_oper_identifiers: Optional[Sequence[int]] = None, fig: Optional[Figure] = None, xscale: str = 'log', yscale: str = 'linear', omega_in_units_of_tau: bool = True, cycler: Optional[Cycler] = None, plot_kw: dict = {}, subplot_kw: Optional[dict] = None, gridspec_kw: Optional[dict] = None, **figure_kw) → Tuple[Figure, Axes, Legend][source]

Plot the fidelity pulse correlation filter functions of the given PulseSequence if they were computed during concatenation for positive frequencies.

Returns a figure with n by n subplots where n is the number of pulses that were concatenated. As of now only the diagonal elements of \(F_{\alpha\beta}\) are implemented, i.e. the filter functions corresponding to uncorrelated noise sources.

Parameters

pulse: PulseSequence: The pulse sequence whose filter function to plot.
n_oper_identifiers: array_like, optional: The identifiers of the noise operators for which the filter function should be plotted. All identifiers can be accessed via pulse.n_oper_identifiers. Defaults to all.
fig: matplotlib figure, optional: A matplotlib figure instance to plot in
xscale: str, optional: x-axis scaling. One of (‘linear’, ‘log’).
yscale: str, optional: y-axis scaling. One of (‘linear’, ‘log’).
omega_in_units_of_tau: bool, optional: Plot \(\omega\tau\) or just \(\omega\) on x-axis.
cycler: cycler.Cycler, optional: A Cycler instance used to set the style cycle if multiple lines are to be drawn in one subplot. Used for all subplots.
plot_kw: dict, optional: Dictionary with keyword arguments passed to the plot function
subplot_kw: dict, optional: Dictionary with keyword arguments passed to the subplots constructor
gridspec_kw: dict, optional: Dictionary with keyword arguments passed to the gridspec constructor
figure_kw: optional: Keyword argument dictionaries that are fed into the matplotlib.pyplot.subplots() function if no fig instance if specified.

Returns

fig: matplotlib figure: The matplotlib figure instance used for plotting.
axes: matplotlib axes: The matplotlib axes instances used for plotting.
legend: matplotlib legend: The matplotlib legend instance in the plot.

Raises

CalculationError: If the pulse correlation filter function was not computed during concatenation.

plot_pulse_train(pulse: PulseSequence, c_oper_identifiers: Optional[Sequence[int]] = None, fig: Optional[Figure] = None, axes: Optional[Axes] = None, cycler: Optional[Cycler] = None, plot_kw: Optional[dict] = {}, subplot_kw: Optional[dict] = None, gridspec_kw: Optional[dict] = None, **figure_kw) → Tuple[Figure, Axes, Legend][source]

Plot the pulsetrain of the PulseSequence pulse.

Parameters

pulse: PulseSequence: The pulse sequence whose pulse train to plot.
c_oper_identifiers: array_like, optional: The identifiers of the control operators for which the pulse train should be plotted. All identifiers can be accessed via pulse.c_oper_identifiers. Defaults to all.
fig: matplotlib figure, optional: A matplotlib figure instance to plot in
axes: matplotlib axes, optional: A matplotlib axes instance to use for plotting.
cycler: cycler.Cycler, optional: A Cycler instance used to set the style cycle if multiple lines are to be drawn
plot_kw: dict, optional: Dictionary with keyword arguments passed to the plot function
subplot_kw: dict, optional: Dictionary with keyword arguments passed to the subplots constructor
gridspec_kw: dict, optional: Dictionary with keyword arguments passed to the gridspec constructor
figure_kw: optional: Keyword argument dictionaries that are fed into the matplotlib.pyplot.subplots() function if no fig instance is specified.

Returns

fig: matplotlib figure: The matplotlib figure instance used for plotting.
axes: matplotlib axes: The matplotlib axes instance used for plotting.
legend: matplotlib legend: The matplotlib legend instance in the plot.

Raises

ValueError: If an invalid number of c_oper_labels were given

2.6. filter_functions.pulse_sequence module

This module defines the PulseSequence class, the package’s central object.

2.6.1. Classes

PulseSequence: The pulse sequence defined by a Hamiltonian

2.6.2. Functions

concatenate(): Function to concatenate different PulseSequence instances and efficiently compute their joint filter function
concatenate_periodic(): Function to more efficiently concatenate many versions of the same PulseSequence instances and compute their joint filter function
extend(): Function to map several PulseSequence instances to different qubits, efficiently scaling up cached attributes.

class PulseSequence(*args, **kwargs)[source]

Bases: object

A class for pulse sequences and their filter functions.

The Hamiltonian is separated into a control and a noise part with

\[\begin{split}\mathcal{H}_c &= \sum_i a_i(t) A_i \\ \mathcal{H}_n &= \sum_j s_j(t) b_j(t) B_j\end{split}\]

where \(A_i\) and \(B_j\) are hermitian operators and \(b_j(t)\) are classically fluctuating noise variables captured in a power spectral density and not needed at instantiation of a PulseSequence.

Parameters

H_c: list of lists

A nested list of n_cops nested lists as taken by QuTiP functions (see for example qutip.propagator.propagator()) describing the control part of the Hamiltonian. The i-th entry of the list should be a list consisting of the i-th operator \(A_i\) making up the control Hamiltonian and a list or array \(a_i(t)\) describing the magnitude of that operator during the time intervals dt. Optionally, the list may also include operator identifiers. That is, H_c should look something like this:

H = [[c_oper1, c_coeff1, c_oper_identifier1],
     [c_oper2, c_coeff2, c_oper_identifier2], ...]

The operators may be given either as NumPy arrays or QuTiP Qobjs and each coefficient array should have the same number of elements as dt, and should be given in units of \(\hbar\). If not every sublist (read: operator) was given a identifier, they are automatically filled up with ‘A_i’ where i is the position of the operator.

H_n: list of lists

A nested list of n_nops nested lists as taken by QuTiP functions (see for example qutip.propagator.propagator()) describing the noise part of the Hamiltonian. The j-th entry of the list should be a list consisting of the j-th operator \(B_j\) making up the noise Hamiltonian and a list or array describing the sensitivity \(s_j(t)\) of the system to the noise operator during the time intervals dt. Optionally, the list may also include operator identifiers. That is, H_n should look something like this:

H = [[n_oper1, n_coeff1, n_oper_identifier1],
     [n_oper2, n_coeff2, n_oper_identifier2], ...]

The operators may be given either as NumPy arrays or QuTiP Qobjs and each coefficient array should have the same number of elements as dt, and should be given in units of \(\hbar\). If not every sublist (read: operator) was given a identifier, they are automatically filled up with ‘B_i’ where i is the position of the operator.

dt: array_like, shape (n_dt,)

The segment durations of the Hamiltonian (i.e. durations of constant control). Internally, the control operation is taken to start at \(t_0\equiv 0\), i.e. the edges of the constant control segments are at times t = [0, *np.cumsum(dt)].

basis: Basis, shape (d**2, d, d), optional

The operator basis in which to calculate. If a Generalized Gell-Mann basis (see ggm()) is chosen, some calculations will be faster for large dimensions due to a simpler basis expansion. However, when extending the pulse sequence to larger qubit registers, cached filter functions cannot be retained since the GGM basis does not factor into tensor products. In this case a Pauli basis is preferable.

Notes

Due to the heavy use of NumPy’s einsum() function, results have a floating point error of ~1e-13.

Examples

A rotation by \(\pi\) around the axis between x and y preceeded and followed by a period of free evolution with the system subject to dephasing noise.

>>> import qutip as qt; import numpy as np
>>> H_c = [[qt.sigmax(), [0, np.pi, 0]],
           [qt.sigmay(), [0, np.pi, 0]]]
>>> # Equivalent pulse:
>>> # H_c = [[qt.sigmax() + qt.sigmay(), [0, np.pi, 0]]]
>>> # The noise sensitivity is constant
>>> H_n = [[qt.sigmaz()/np.sqrt(2), [1, 1, 1], 'Z']]
>>> dt = [1, 1, 1]
>>> # Free evolution between t=0 and t=1, rotation between t=1 and t=2,
>>> # and free evolution again from t=2 to t=3.
>>> pulse = PulseSequence(H_c, H_n, dt)
>>> pulse.c_oper_identifiers
['A_0', 'A_1']
>>> pulse.n_oper_identifiers
['Z']
>>> omega = np.logspace(-1, 2, 500)
>>> F = pulse.get_filter_function(omega)    # shape (1, 500)
>>> # Plot the resulting filter function:
>>> from filter_functions import plotting
>>> fig, ax, leg = plotting.plot_filter_function(pulse)

Attributes

c_opers: ndarray, shape (n_cops, d, d): Control operators. Note that they are stored sorted by their corresponding identifiers.
n_opers: ndarray, shape (n_nops, d, d): Noise operators. Note that they are stored sorted by their corresponding identifiers.
c_oper_identifers: sequence of str: Identifiers for the control operators of the system. Stored sorted.
n_oper_identifers: sequence of str: Identifiers for the noise operators of the system. Stored sorted.
c_coeffs: ndarray, shape (n_cops, n_dt): Control parameters in units of \(\hbar\). Note that they are stored sorted by their corresponding identifiers.
n_coeffs: ndarray, shape (n_nops, n_dt): Noise sensitivities in units of \(\hbar\). Note that they are stored sorted by their corresponding identifiers.
dt: ndarray, shape (n_dt,): Time steps
t: ndarray, shape (n_dt + 1,): Absolute times taken to start at \(t_0\equiv 0\)
tau: float: Total duration. Equal to t[-1].
d: int: Dimension of the Hamiltonian
basis: Basis, shape (d**2, d, d): The operator basis used for calculation
nbytes: int: An estimate of the memory consumed by the PulseSequence instance and its attributes
If the Hamiltonian was diagonalized, the eigenvalues and -vectors as
well as the cumulative propagators are cached:
eigvals: ndarray, shape (n_dt, d): Eigenvalues \(D^{(g)}\)
eigvecs: ndarray, shape (n_dt, d, d): Eigenvectors \(V^{(g)}\)
propagators: ndarray, shape (n_dt+1, d, d): Cumulative propagators \(Q_g\)
total_propagator: ndarray, shape (d, d): The total propagator \(Q\) of the pulse alone. That is, \(|\psi(\tau)\rangle = Q|\psi(0)\rangle\).
total_propagator_liouville: array_like, shape (d**2, d**2): The transfer matrix for the total propagator of the pulse. Given by liouville_representation(pulse.total_propagator, pulse.basis).
Furthermore, when the filter function is calculated, the frequencies
are cached as well as other relevant quantities.

Methods

get_filter_function_derivative(omega[, ...])

Calculate the pulse sequence's filter function derivative.

cleanup(method=’conservative’)	Delete cached attributes
is_cached(attr)	Checks if a given attribute of the `PulseSequence` is cached
diagonalize()	Diagonalize the Hamiltonian of the pulse sequence, computing eigenvalues and -vectors as well as cumulative propagators
get_control_matrix(omega, show_progressbar=False)	Calculate the control matrix for frequencies omega
get_filter_function(omega, which=’fidelity’, show_progressbar=False)	Calculate the filter function for frequencies omega
get_pulse_correlation_filter_function(which=’fidelity’)	Get the pulse correlation filter function (only possible if computed during concatenation)
propagator_at_arb_t(t)	Calculate the propagator at arbitrary times

Initialize a PulseSequence instance.

cache_control_matrix(omega: Sequence[float], control_matrix: Optional[ndarray] = None, show_progressbar: bool = False, cache_intermediates: bool = False) → None[source]

Cache the control matrix and the frequencies it was calculated for.

Parameters

omega: array_like, shape (n_omega,): The frequencies for which to cache the filter function.
control_matrix: array_like, shape ([n_nops,] n_nops, d**2, n_omega), optional: The control matrix for the frequencies omega. If None, it is computed.
show_progressbar: bool: Show a progress bar for the calculation of the control matrix.
cache_intermediates: bool, optional: Keep intermediate terms of the calculation that are also required by other computations. Only applies if control_matrix is not supplied.

cache_filter_function(omega: Sequence[float], control_matrix: Optional[ndarray] = None, filter_function: Optional[ndarray] = None, which: str = 'fidelity', order: int = 1, show_progressbar: bool = False, cache_intermediates: bool = False) → None[source]

Cache the filter function. If control_matrix.ndim == 4, it is taken to be the ‘pulse correlation control matrix’ and summed along the first axis. In that case, also the pulse correlation filter function is calculated and cached. Total phase factors and transfer matrices of the the cumulative propagator are also cached so they can be reused during concatenation.

Parameters

omega: array_like, shape (n_omega,): The frequencies for which to cache the filter function.
control_matrix: array_like, shape ([n_nops,] n_nops, d**2, n_omega), optional: The control matrix for the frequencies omega. If None, it is computed and the filter function derived from it.
filter_function: array_like, shape (n_nops, n_nops, [d**2, d**2,] n_omega), optional: The filter function for the frequencies omega. If None, it is computed from R in the case order == 1 and from scratch else.
which: str, optional: Which filter function to cache. Either ‘fidelity’ (default) or ‘generalized’.
order: int, optional: First or second order filter function.
show_progressbar: bool, optional: Show a progress bar for the calculation of the control matrix.
cache_intermediates: bool, optional: Keep intermediate terms of the calculation that are also required by other computations.

See also

PulseSequence.get_filter_function: Getter method

cache_total_phases(omega: Sequence[float], total_phases: Optional[ndarray] = None) → None[source]

Cache the total phase factors for this pulse and omega.

Parameters

omega: array_like, shape (n_omega,): The frequencies for which to cache the phase factors.
total_phases: array_like, shape (n_omega,), optional: The total phase factors for the frequencies omega. If None, they are computed.

cleanup(method: str = 'conservative') → None[source]

Delete cached byproducts of the calculation of the filter function that are not necessarily needed anymore in order to free up memory.

Parameters

method: optional

If set to ‘conservative’ (the default), only the following attributes are deleted:

_eigvals

_eigvecs

_propagators

If set to ‘greedy’, all of the above as well as the following attributes are deleted:

_total_propagator

_total_propagator_liouville

_total_phases

_control_matrix

_control_matrix_pc

If set to ‘all’, all of the above as well as the following attributes are deleted:

omega

_filter_function

_filter_function_gen

_filter_function_pc

_filter_function_pc_gen

_filter_function_2

_intermediates[‘control_matrix_step’]

If set to ‘frequency dependent’ only attributes that are functions of frequency are initalized to None.

Note that if this PulseSequence is concatenated with another one, some of the attributes might need to be calculated again, resulting in slower execution of the concatenation.

diagonalize() → None[source]: Diagonalize the Hamiltonian defining the pulse sequence.

property duration: float: The duration of the pulse. Alias of tau.

property eigvals: ndarray: Get the eigenvalues of the pulse’s Hamiltonian.

property eigvecs: ndarray: Get the eigenvectors of the pulse’s Hamiltonian.

get_control_matrix(omega: Sequence[float], show_progressbar: bool = False, cache_intermediates: bool = False) → ndarray[source]

Get the control matrix for the frequencies omega. If it has been cached for the same frequencies, the cached version is returned, otherwise it is calculated from scratch.

Parameters

omega: array_like, shape (n_omega,): The frequencies at which to evaluate the control matrix.
show_progressbar: bool: Show a progress bar for the calculation of the control matrix.
cache_intermediates: bool, optional: Keep intermediate terms of the calculation that are also required by other computations.

Returns

control_matrix: ndarray, shape (n_nops, d**2, n_omega): The control matrix for the noise operators.

get_filter_function(omega: Sequence[float], which: str = 'fidelity', order: int = 1, show_progressbar: bool = False, cache_intermediates: bool = False) → ndarray[source]

Get the first or second order filter function.

The filter function is cached so it doesn’t need to be calculated twice for the same frequencies.

Parameters

omega: array_like, shape (n_omega,): The frequencies at which to evaluate the filter function.
which: str, optional: Which filter function to return. Either ‘fidelity’ (default) or ‘generalized’ (see Notes). Only if order == 1.
order: int, optional: First or second order filter function.
show_progressbar: bool, optional: Show a progress bar for the calculation of the control matrix.
cache_intermediates: bool, optional: Keep intermediate terms of the calculation that are also required by other computations.

Returns

filter_function: ndarray, shape (n_nops, n_nops, [d**2, d**2,] n_omega): The filter function for each combination of noise operators as a function of omega.

Notes

The first-order generalized filter function is given by

\[F_{\alpha\beta,kl}(\omega) = \tilde{\mathcal{B}}_{\alpha k}^\ast(\omega) \tilde{\mathcal{B}}_{\beta l}(\omega),\]

where \(\alpha,\beta\) are indices counting the noise operators \(B_\alpha\) and \(k,l\) indices counting the basis elements \(C_k\).

The fidelity filter function is obtained by tracing over the basis indices:

\[F_{\alpha\beta}(\omega) = \sum_{k} F_{\alpha\beta,kk}(\omega).\]

get_filter_function_derivative(omega: Sequence[float], control_identifiers: Optional[Sequence[str]] = None, n_oper_identifiers: Optional[Sequence[str]] = None, n_coeffs_deriv: Optional[Sequence[Sequence[float]]] = None) → ndarray[source]

Calculate the pulse sequence’s filter function derivative.

Parameters

omega: array_like, shape (n_omega,): Frequencies at which the pulse control matrix is to be evaluated.
control_identifiers: Sequence[str], shape (n_ctrl,): Sequence of strings with the control identifiers to distinguish between control and drift Hamiltonian. The default is None, in which case the derivative is computed for all known non-noise operators.
n_oper_identifiers: Sequence[str], shape (n_nops,): Sequence of strings with the noise identifiers for which to compute the derivative contribution. The default is None, in which case it is computed for all known noise operators.
n_coeffs_deriv: array_like, shape (n_nops, n_ctrl, n_dt): The derivatives of the noise susceptibilities by the control amplitudes. The rows and columns should be in the same order as the corresponding identifiers above. Defaults to None, in which case the coefficients are assumed to be constant and hence their derivative vanishing.

Warning

Internally, control and noise terms of the Hamiltonian are stored alphanumerically sorted by their identifiers. If the noise and/or control identifiers above are not explicitly given, the rows and/or columns of this parameter need to be sorted in the same fashion.

Returns

filter_function_deriv: ndarray, shape (n_nops, n_t, n_ctrl, n_omega): The regular filter functions’ derivatives for variation in each control contribution. Sorted in the same fashion as n_coeffs_deriv or, if not given, alphanumerically by the identifiers.

get_pulse_correlation_control_matrix() → ndarray[source]: Get the pulse correlation control matrix if it was cached.

get_pulse_correlation_filter_function(which: str = 'fidelity') → ndarray[source]

Get the pulse correlation filter function given by

\[F_{\alpha\beta}^{(gg')}(\omega) = e^{i\omega(t_{g-1} - t_{g'-1})} \tilde{\mathcal{B}}^{(g)}(\omega)\mathcal{Q}^{(g-1)} \mathcal{Q}^{(g'-1)\dagger} \tilde{\mathcal{B}}^{(g')\dagger}(\omega),\]

where \(g,g'\) index the pulse in the sequence and \(\alpha,\beta\) index the noise operators, if it was computed during concatenation. Since the calculation requires the individual pulse’s control matrices and phase factors, which are not retained after concatenation, the pulse correlation filter function cannot be computed afterwards.

Note that the frequencies for which the filter function was calculated are not stored.

Returns

filter_function_pc: ndarray, shape (n_pls, n_pls, n_nops, n_nops, n_omega)

The pulse correlation filter function for each noise operator as a function of omega. The first two axes correspond to the pulses in the sequence, i.e. if the concatenated pulse sequence is \(C\circ B\circ A\), the first two axes are arranged like

\[\begin{split}F_{\alpha\beta}^{(gg')} &= \begin{pmatrix} F_{\alpha\beta}^{(AA)} & F_{\alpha\beta}^{(AB)} & F_{\alpha\beta}^{(AC)} \\ F_{\alpha\beta}^{(BA)} & F_{\alpha\beta}^{(BB)} & F_{\alpha\beta}^{(BC)} \\ F_{\alpha\beta}^{(CA)} & F_{\alpha\beta}^{(CB)} & F_{\alpha\beta}^{(CC)} \end{pmatrix}\end{split}\]

for \(g,g'\in\{A, B, C\}\).

get_total_phases(omega: Sequence[float]) → ndarray[source]: Get the (cached) total phase factors for this pulse and omega.

is_cached(attr: str) → bool[source]: Returns True if the attribute is cached

property nbytes: int: Return an estimate of the amount of memory consumed by this object (or, more precisely, the array attributes of this object).

property omega: ndarray: Cached frequencies

propagator_at_arb_t(t: Sequence[float]) → ndarray[source]: Calculate the cumulative propagator Q(t) at times t by making use of the fact that we assume piecewise-constant control.

property propagators: ndarray: Get the eigenvectors of the pulse’s Hamiltonian.

property t: ndarray: The times of the pulse.

property tau: Union[float, int]: The duration of the pulse.

property total_propagator: ndarray: Get total propagator of the pulse.

property total_propagator_liouville: ndarray: Get the transfer matrix for the total propagator of the pulse.

concatenate(pulses: Iterable[PulseSequence], calc_pulse_correlation_FF: bool = False, calc_filter_function: Optional[bool] = None, which: str = 'fidelity', omega: Optional[Sequence[float]] = None, show_progressbar: bool = False) → PulseSequence[source]

Concatenate an arbitrary number of pulses. Note that pulses are concatenated left-to-right, that is,

\[\mathtt{concatenate((A, B))} \equiv B \circ A\]

so that \(A\) is executed before \(B\) when applying the concatenated pulse.

Parameters

pulses: sequence of PulseSequences: The PulseSequence instances to be concatenated. If any of the instances have a cached filter function, the filter function for the composite pulse will also be calculated in order to make use of the speedup gained from concatenating the filter functions. If omega is given, calculation of the composite filter function is forced.
calc_pulse_correlation_FF: bool, optional: Switch to control whether the pulse correlation filter function (see PulseSequence.get_pulse_correlation_filter_function()) is calculated. If omega is not given, the cached frequencies of all pulses need to be equal.
calc_filter_function: bool, optional: Switch to force the calculation of the filter function to be carried out or not. Overrides the automatic behavior of calculating it if at least one pulse has a cached control matrix. If True and no pulse has a cached control matrix, a list of frequencies must be supplied as omega.
which: str, optional: Which filter function to compute. Either ‘fidelity’ (default) or ‘generalized’ (see PulseSequence.get_filter_function() and PulseSequence.get_pulse_correlation_filter_function()).
omega: array_like, optional: Frequencies at which to evaluate the (pulse correlation) filter functions. If None, an attempt is made to use cached frequencies.
show_progressbar: bool: Show a progress bar for the calculation of the control matrix.

Returns

pulse: PulseSequence: The concatenated pulse.

concatenate_periodic(pulse: PulseSequence, repeats: int) → PulseSequence[source]

Concatenate a pulse sequence pulse whose Hamiltonian is periodic repeats times. Although performing the same task, this function is much faster for concatenating many identical pulses with filter functions than concatenate().

Note that for large dimensions, the calculation of the control matrix using this function might be very memory intensive.

Parameters

pulse: PulseSequence: The PulseSequence instance to be repeated. If it has a cached filter function, the filter function for the new pulse will also be computed.
repeats: int: The number of repetitions

Returns

newpulse: PulseSequence: The concatenated PulseSequence

See also

concatenate: Concatenate arbitrary PulseSequences.

Notes

The total control matrix is given by

\[\begin{split}\tilde{\mathcal{B}}(\omega) &= \tilde{\mathcal{B}}^{(1)}(\omega) \sum_{g=0}^{G-1} \left(e^{i\omega T}\right)^g \\ &= \tilde{\mathcal{B}}^{(1)}(\omega)\bigl( \mathbb{I} - e^{i\omega T} \mathcal{Q}^{(1)}\bigr)^{-1}\bigl( \mathbb{I} - \bigl(e^{i\omega T} \mathcal{Q}^{(1)}\bigr)^G\bigr).\end{split}\]

with \(T\) the period of the control Hamiltonian and \(G\) the number of periods. The last equality is valid only if \(\mathbb{I} - e^{i\omega T}\mathcal{Q}^{(1)}\) is invertible.

extend(pulse_to_qubit_mapping: Sequence[Sequence[Optional[Union[PulseSequence, Sequence[int], int, Mapping[str, str]]]]], N: Optional[int] = None, d_per_qubit: int = 2, additional_noise_Hamiltonian: Optional[Sequence[Sequence[Union[ndarray, Qobj, Sequence[float]]]]] = None, cache_diagonalization: Optional[bool] = None, cache_filter_function: Optional[bool] = None, omega: Optional[Sequence[float]] = None, show_progressbar: bool = False) → PulseSequence[source]

Map one or more pulse sequences to different qubits.

Parameters

pulse_to_qubit_mapping: sequence of mapping tuples

A sequence of tuples with the first entry a PulseSequence instance and the second an int or tuple of ints indicating the qubits that the PulseSequence should be mapped to. A mapping of operator identifiers may optionally be given as a third element of each tuple. By default, the index of the qubit the operator is mapped to is appended to its identifier.

Pulse sequences defined for multiple qubits may also be extended to non-neighboring qubits. Note that for multi-qubit pulses the order of the qubits is respected, i.e. mapping a pulse to (1, 0) is different from mapping it to (0, 1).

N: int

The total number of qubits the new PulseSequence should be defined for. By default, this is inferred from pulse_to_qubit_mapping.

d_per_qubit: int

The size of the Hilbert space a single qubit requires.

additional_noise_Hamiltonian: list of lists

Additional noise operators and corresponding sensitivities for the new pulse sequence.

cache_diagonalization: bool

Force diagonalizing the new pulse sequence. By default, diagonalization is cached if all pulses in pulse_to_qubit_mapping have been diagonalized since it is much cheaper to get the relevant quantities as tensor products from the mapped pulses instead of diagonalizing the new pulse.

cache_filter_function: bool

Force computing the filter functions for the new pulse sequence. Noise operators of individual pulses will be extended to the new Hilbert space. By default, this is done if all pulses in pulse_to_qubit_mapping have their filter functions cached.

Note that extending the filter functions is only possible if they the mapped pulses are using a separable basis like the Pauli basis.

omega: array_like

Frequencies for which to compute the filter functions if cache_filter_function == True. Defaults to None, in which case the cached frequencies of the individual pulses need to be the same.

show_progressbar: bool

Show a progress bar for the calculation of the control matrix.

Returns

newpulse: PulseSequence: The new pulse sequence on the larger qubit register. The noise operators (and possibly filter functions) are stored in the following order: first those of the multi-qubit pulses in the order they appeared in pulse_to_qubit_mapping, then those of the single-qubit pulses, and lastly any additional ones that may be given by additional_noise_Hamiltonian.

See also

remap: Map PulseSequence to a different qubit.
concatenate: Concatenate PulseSequences (in time).
concatenate_periodic: Periodically concatenate a PulseSequence.

Examples

>>> import filter_functions as ff
>>> I, X, Y, Z = ff.util.paulis
>>> X_pulse = ff.PulseSequence([[X, [np.pi/2], 'X']],
...                            [[X, [1], 'X'], [Z, [1], 'Z']],
...                            [1], basis=ff.Basis.pauli(1))
>>> XX_pulse = ff.extend([(X_pulse, 0), (X_pulse, 1)])
>>> XX_pulse.d
4
>>> XIX_pulse_1 = ff.extend([(X_pulse, 0), (X_pulse, 2)])
>>> XIX_pulse_1.d
8
>>> XXI_pulse = ff.extend([(X_pulse, 0), (X_pulse, 1)], N=3)
>>> XXI_pulse.d
8

Filter functions are automatically cached if they are for mapped pulses:

>>> omega = ff.util.get_sample_frequencies(X_pulse)
>>> X_pulse.cache_filter_function(omega)
>>> XX_pulse = ff.extend([(X_pulse, 0), (X_pulse, 1)])
>>> XX_pulse.is_cached('filter_function')
True

This behavior can also be overriden manually:

>>> XX_pulse = ff.extend([(X_pulse, 0), (X_pulse, 1)],
...                      cache_filter_function=False)
>>> XX_pulse.is_cached('filter_function')
False

Mapping pulses to non-neighboring qubits is also possible:

>>> Y_pulse = ff.PulseSequence([[Y, [np.pi/2], 'Y']],
...                            [[Y, [1], 'Y'], [Z, [1], 'Z']],
...                            [1], basis=ff.Basis.pauli(1))
>>> XXY_pulse = ff.extend([(XX_pulse, (0, 1)), (Y_pulse, 2)])
>>> XYX_pulse = ff.extend([(XX_pulse, (0, 2)), (Y_pulse, 1)])

Additionally, pulses can have the order of the qubits they are defined for permuted (see remap()):

>>> Z_pulse = ff.PulseSequence([[Z, [np.pi/2], 'Z']], [[Z, [1], 'Z']],
...                            [1], basis=ff.Basis.pauli(1))
>>> XY_pulse = ff.extend([(X_pulse, 0), (Y_pulse, 1)])
>>> YZX_pulse = ff.extend([(XY_pulse, (2, 0)), (Z_pulse, 1)])

Control and noise operator identifiers can be mapped according to a specified mapping:

>>> YX_pulse = ff.extend([(X_pulse, 1, {'X': 'IX', 'Z': 'IZ'}),
...                       (Y_pulse, 0, {'Y': 'YI', 'Z': 'ZI'})])
>>> YX_pulse.c_oper_identifiers
array(['IX', 'YI'], dtype='<U2')
>>> YX_pulse.n_oper_identifiers
array(['IX', 'IZ', 'YI', 'ZI'], dtype='<U2')

We can also add an additional noise Hamiltonian:

>>> H_n = [[ff.util.tensor(Z, Z, Z), [1], 'ZZZ']]
>>> XYX_pulse = ff.extend([(XX_pulse, (0, 2)), (Y_pulse, 1)],
...                       additional_noise_Hamiltonian=H_n)
>>> 'ZZZ' in XYX_pulse.n_oper_identifiers
True

remap(pulse: PulseSequence, order: Sequence[int], d_per_qubit: int = 2, oper_identifier_mapping: Optional[Mapping[str, str]] = None) → PulseSequence[source]

Remap a PulseSequence by changing the order of qubits in the register. Cached attributes are automatically attempted to be retained.

Caution

This function simply permutes the order of the tensor product elements of control and noise operators. Thus, the resultant pulse will have its filter functions defined for different noise operators than the original one.

Parameters

pulse: PulseSequence: The pulse whose qubit order should be permuted.
order: sequence of ints: A list of permutation indices. E.g., if pulse is defined for two qubits, order == [1, 0] will reverse the order of qubits.
d_per_qubit: int (default: 2): The size of the Hilbert space a single qubit inhabitates.
oper_identifier_mapping: dict_like: A mapping that maps operator identifiers from the old pulse to the remapped pulse. The default is the identity mapping.

Returns

remapped_pulse: PulseSequence: A new PulseSequence instance with the order of the qubits permuted according to order.

See also

extend: Map PulseSequences to composite Hilbert spaces.
util.tensor_transpose: Transpose the order of a tensor product.

Examples

>>> X, Y = util.paulis[1:3]
>>> XY, YX = util.tensor(X, Y), util.tensor(Y, X)
>>> pulse = PulseSequence([[XY, [np.pi/2], 'XY']], [[YX, [1], 'YX']], [1],
...                       Basis.pauli(2))
>>> mapping = {'XY': 'YX', 'YX': 'XY'}
>>> remapped_pulse = remap(pulse, (1, 0), oper_identifier_mapping=mapping)
>>> target_pulse = PulseSequence([[YX, [np.pi/2], 'YX']],
...                              [[XY, [1], 'XY']], [1], Basis.pauli(2))
>>> remapped_pulse == target_pulse
True

Caching of attributes is automatically handled >>> remapped_pulse.is_cached(‘filter_function’) False >>> pulse.cache_filter_function(util.get_sample_frequencies(pulse)) >>> remapped_pulse = remap(pulse, (1, 0)) >>> remapped_pulse.is_cached(‘filter_function’) True

2.7. filter_functions.gradient module

This module implements functions to calculate filter function and infidelity derivatives. Currently only auto-correlated noise (i.e. no cross-correlations) is implemented.

Throughout this documentation the following notation will be used:

n_dt denotes the number of time steps,
n_cops the number of all control operators,
n_ctrl the number of accessible control operators (if identifiers are provided, otherwise n_ctrl=n_cops),
n_nops the number of noise operators,
n_omega the number of frequency samples, and
d the dimension of the system.

2.7.1. Functions

calculate_derivative_of_control_matrix_from_scratch(): Calculate the derivative of the control matrix from scratch.
calculate_canonical_filter_function_derivative(): Compute the filter function derivative from the control matrix.
infidelity_derivative(): Calculate the infidelity derivative.

calculate_derivative_of_control_matrix_from_scratch(omega: Sequence[float], propagators: ndarray, eigvals: ndarray, eigvecs: ndarray, basis: Basis, t: Sequence[float], dt: Sequence[float], n_opers: Sequence[Union[ndarray, Qobj]], n_coeffs: Sequence[Sequence[float]], c_opers: Sequence[Union[ndarray, Qobj]], n_coeffs_deriv: Optional[Sequence[Sequence[float]]] = None, intermediates: Optional[Dict[str, ndarray]] = None) → ndarray[source]

Calculate the derivative of the control matrix from scratch.

Parameters

omega: array_like, shape (n_omega,): Frequencies, at which the pulse control matrix is to be evaluated.
propagators: array_like, shape (n_dt+1, d, d): The propagators \(Q_g = P_g P_{g-1}\cdots P_0\) as a (d, d) array with d the dimension of the Hilbert space.
eigvals: array_like, shape (n_dt, d): Eigenvalue vectors for each time pulse segment g with the first axis counting the pulse segment, i.e. HD == array([D_0, D_1, ...]).
eigvecs: array_like, shape (n_dt, d, d): Eigenvector matrices for each time pulse segment g with the first axis counting the pulse segment, i.e. HV == array([V_0, V_1, ...]).
basis: Basis, shape (d**2, d, d): The basis elements, in which the pulse control matrix will be expanded.
t: array_like, shape (n_dt+1), optional: The absolute times of the different segments.
dt: array_like, shape (n_dt): Sequence duration, i.e. for the \(g\)-th pulse \(t_g - t_{g-1}\).
n_opers: array_like, shape (n_nops, d, d): Noise operators \(B_\alpha\).
n_coeffs: array_like, shape (n_nops, n_dt): The sensitivities of the system to the noise operators given by n_opers at the given time step.
c_opers: array_like, shape (n_ctrl, d, d): Control operators \(H_k\) with respect to which the derivative is computed.
n_coeffs_deriv: array_like, shape (n_nops, n_ctrl, n_dt): The derivatives of the noise susceptibilities by the control amplitudes. Defaults to None.
intermediates: Dict[str, ndarray], optional: Optional dictionary containing intermediate results of the calculation of the control matrix.

Returns

ctrlmat_deriv: ndarray, shape (n_ctrl, n_omega, n_dt, n_nops, d**2): Partial derivatives of the total control matrix with respect to each control direction \(\frac{\partial\mathcal{B}_{\alpha k}(\omega)}{\partial u_h(t_{g'})}\).

See also

_liouville_derivative
_control_matrix_at_timestep_derivative

Notes

The derivative of the control matrix is calculated according to

\[\frac{\partial\mathcal{B}_{\alpha k}(\omega)}{\partial u_h(t_{g'})} = \sum_{g=1}^G \mathrm{e}^{\mathrm{i}\omega t_{g-1}}\cdot\left(\sum_j \left[\frac{\partial\mathcal{B}_{\alpha j}^{(g)}(\omega)} {\partial u_h(t_{g'})} \cdot \mathcal{Q}_{jk}^{(g-1)} + \mathcal{B}_{\alpha j}^{(g)}(\omega) \cdot\frac{\partial \mathcal{Q}_{jk}^{(g-1)}} {\partial u_h(t_{g'})} \right] \right)\]

calculate_filter_function_derivative(ctrlmat: ndarray, ctrlmat_deriv: ndarray) → ndarray[source]

Compute the filter function derivative from the control matrix.

Parameters

ctrlmat: array_like, shape (n_nops, d**2, n_omega): The control matrix.
ctrlmat_deriv: array_like, shape (n_nops, d**2, n_t, n_ctrl, n_omega): The derivative of the control matrix.

Returns

filter_function_deriv: ndarray, shape (n_nops, n_dt, n_ctrl, n_omega): The regular filter functions’ derivatives for variation in each control contribution \(\frac{\partial F_\alpha(\omega)}{\partial u_h(t_{g'})}\).

Notes

The filter function derivative is calculated according to

\[\frac{\partial F_\alpha(\omega)}{\partial u_h(t_{g'})} = 2 \mathrm{Re}\left(\sum_k \mathcal{B}_{\alpha k}^\ast(\omega) \frac{\partial\mathcal{B}_{\alpha k}(\omega)} {\partial u_h(t_{g'})} \right)\]

infidelity_derivative(pulse: PulseSequence, spectrum: Sequence[float], omega: Sequence[float], control_identifiers: Optional[Sequence[str]] = None, n_oper_identifiers: Optional[Sequence[str]] = None, n_coeffs_deriv: Optional[Sequence[Sequence[float]]] = None) → ndarray[source]

Calculate the entanglement infidelity derivative of the PulseSequence pulse.

Parameters

pulse: PulseSequence: The PulseSequence instance for which to calculate the infidelity.
spectrum: array_like, shape ([[n_nops,] n_nops,] omega): The two-sided noise power spectral density in units of inverse frequencies as an array of shape (n_omega,), (n_nops, n_omega), or (n_nops, n_nops, n_omega). In the first case, the same spectrum is taken for all noise operators, in the second, it is assumed that there are no correlations between different noise sources and thus there is one spectrum for each noise operator. In the third and most general case, there may be a spectrum for each pair of noise operators corresponding to the correlations between them. n_nops is the number of noise operators considered and should be equal to len(n_oper_identifiers).
omega: array_like, shape (n_omega,): The frequencies at which the integration is to be carried out.
control_identifiers: Sequence[str], shape (n_ctrl,): Sequence of strings with the control identifiers to distinguish between control and drift Hamiltonian. The default is None, in which case the derivative is computed for all known non-noise operators.
n_oper_identifiers: Sequence[str], shape (n_nops,): Sequence of strings with the noise identifiers for which to compute the derivative contribution. The default is None, in which case it is computed for all known noise operators.
n_coeffs_deriv: array_like, shape (n_nops, n_ctrl, n_dt): The derivatives of the noise susceptibilities by the control amplitudes. The rows and columns should be in the same order as the corresponding identifiers above. Defaults to None, in which case the coefficients are assumed to be constant and hence their derivative vanishing.

Warning

Internally, control and noise terms of the Hamiltonian are stored alphanumerically sorted by their identifiers. If the noise and/or control identifiers above are not explicitly given, the rows and/or columns of this parameter need to be sorted in the same fashion.

Returns

infid_deriv: ndarray, shape (n_nops, n_dt, n_ctrl): Array with the derivative of the infidelity for each noise source taken for each control direction at each time step \(\frac{\partial I_e}{\partial u_h(t_{g'})}\). Sorted in the same fashion as n_coeffs_deriv or, if not given, alphanumerically by the identifiers.

Raises

ValueError: If the provided noise spectral density does not fit expected shape.

Notes

The infidelity’s derivative is given by

\[\frac{\partial I_e}{\partial u_h(t_{g'})} = \frac{1}{d} \int_{-\infty}^\infty \frac{d\omega}{2\pi} S_\alpha(\omega) \frac{\partial F_\alpha(\omega)} {\partial u_h(t_{g'})}\]

with \(S_{\alpha}(\omega)\) the noise spectral density and \(F_{\alpha}(\omega)\) the canonical filter function for noise source \(\alpha\).

To convert to the average gate infidelity, use the following relation given by Horodecki et al. [Hor99] and Nielsen [Nie02]:

\[\big\langle\mathcal{I}_\mathrm{avg}\big\rangle = \frac{d}{d+1} \big\langle\mathcal{I}_\mathrm{e}\big\rangle.\]

References

Hor99: Horodecki, M., Horodecki, P., & Horodecki, R. (1999). General teleportation channel, singlet fraction, and quasidistillation. Physical Review A - Atomic, Molecular, and Optical Physics, 60(3), 1888–1898. https://doi.org/10.1103/PhysRevA.60.1888
Nie02: Nielsen, M. A. (2002). A simple formula for the average gate fidelity of a quantum dynamical operation. Physics Letters, Section A: General, Atomic and Solid State Physics, 303(4), 249–252. https://doi.org/10.1016/S0375-9601(02)01272-0

2.8. filter_functions.superoperator module

This module provides some functions related to superoperators and quantum maps.

2.8.1. Functions

liouville_representation(): Calculate the Liouville representation of a unitary with respect to a basis
liouville_to_choi(): Convert from Liouville to Choi matrix representation.
liouville_is_CP(): Check if superoperator in Liouville representation is completely positive.
liouville_is_cCP(): Check if superoperator in Liouville representation is conditional CP.

liouville_is_CP(superoperator: ndarray, basis: Basis, return_eig: Optional[bool] = False, atol: Optional[float] = None) → Union[bool, Tuple[bool, Tuple[ndarray, ndarray]]][source]

Test if a Liouville superoperator is completely positive (CP).

Parameters

superoperator: ndarray, shape (…, d**2, d**2): The superoperator in Liouville representation to be checked for CPness.
basis: Basis, shape (d**2, d, d): The operator basis defining the Liouville representation.
return_eig: bool, optional: Return the tuple of eigenvalues and eigenvectors of the Choi matrix. The default is False.
atol: float, optional: Absolute tolerance for the complete positivity.

Returns

CP: bool, (shape (…,)): The (array, if broadcasted) of bools indicating if superoperator is CP.
(D, V): Tuple[ndarray, ndarray]: The eigenvalues and eigenvectors of the Choi matrix (only if return_eig is True).

See also

liouville_representation: Calculate Liouville representation of a unitary.
Liouville_to_choi: Convert from Liouville to Choi matrix representation.
liouville_is_cCP: Test if a superoperator is conditional CP.

Notes

A superoperator \(\mathcal{S}\) is completely positive (CP) if and only if its Choi matrix representation is positive semidefinite:

\[\mathcal{S}\text{ is CP }\Leftrightarrow \mathrm{choi}(\mathcal{S})\geq 0.\]

liouville_is_cCP(superoperator: ndarray, basis: Basis, return_eig: Optional[bool] = False, atol: Optional[float] = None) → Union[bool, Tuple[bool, Tuple[ndarray, ndarray]]][source]

Test if a Liouville superoperator is conditional completely positive.

Parameters

superoperator: ndarray, shape (…, d**2, d**2): The superoperator in Liouville representation to be checked for cCPness
basis: Basis, shape (d**2, d, d): The operator basis defining the Liouville representation.
return_eig: bool, optional: Return the tuple of eigenvalues and eigenvectors of the Choi matrix projected on the complement of the maximally entangled state. The default is False.
atol: float, optional: Absolute tolerance for the complete positivity.

Returns

cCP: bool, (shape (…,)): The (array, if broadcasted) of bools indicating if superoperator is cCP
(D, V): Tuple[ndarray, ndarray]: The eigenvalues and eigenvectors of the projected Choi matrix (only if return_eig is True).

See also

liouville_representation: Calculate Liouville representation of a unitary.
Liouville_to_choi: Convert from Liouville to Choi matrix representation.
liouville_is_CP: Test if a superoperator is CP.

Notes

A superoperator \(\mathcal{S}\) is conditional completely positive (cCP) if and only if its Choi matrix projected on the complement of the maximally entangled state is positive semidefinite:

\[\mathcal{S}\text{ is cCP }\Leftrightarrow Q\mathrm{choi}(\mathcal{S})Q\geq 0\]

with \(Q = \mathbb{I} - |\Omega\rangle\langle\Omega|\).

liouville_representation(U: ndarray, basis: Basis) → ndarray[source]

Get the Liouville representaion of the unitary U with respect to the basis.

Parameters

U: ndarray, shape (…, d, d): The unitary.
basis: Basis, shape (d**2, d, d): The basis used for the representation, e.g. a Pauli basis.

Returns

R: ndarray, shape (…, d**2, d**2): The Liouville representation of U.

Notes

The Liouville representation of a unitary quantum operation \(\mathcal{U}:\rho\rightarrow U\rho U^\dagger\) is given by

\[\mathcal{U}_{ij} = \mathrm{tr}(C_i U C_j U^\dagger)\]

with \(C_i\) elements of the basis spanning \(\mathbb{C}^{d\times d}\) with \(d\) the dimension of the Hilbert space.

liouville_to_choi(superoperator: ndarray, basis: Basis) → ndarray[source]

Convert from Liouville to Choi matrix representation.

Parameters

superoperator: ndarray, shape (…, d**2, d**2): The Liouville representation of a superoperator.
basis: Basis, shape (d**2, d, d): The operator basis defining the Liouville representation.

Returns

choi: ndarray, shape (…, d**2, d**2): The Choi matrix representation of the superoperator.

See also

liouville_representation: Calculate Liouville representation of a unitary.
liouville_is_CP: Test if a superoperator is completely positive (CP).
liouville_is_cCP: Test if a superoperator is conditional CP.

Notes

The Choi matrix is given by

\[\begin{split}\mathrm{choi}(\mathcal{S}) &= (\mathbb{I}\otimes\mathcal{S}) (|\Omega\rangle\langle\Omega|) \\ &= \sum_{ij} E_{ij}\otimes\mathcal{S}(E_{ij}) \\ &= \sum_{ij}\mathcal{S}_{ij} C_j^T\otimes C_i\end{split}\]

where \(|\Omega\rangle\) is a maximally entangled state, \(E_{ij} = |i\rangle\langle j|\), and \(C_i\) are the basis elements that define the Liouville representation \(\mathcal{S}_{ij}\) [Mer13].

References

Mer13: Merkel, S. T. et al. Self-consistent quantum process tomography. Physical Review A - Atomic, Molecular, and Optical Physics, 87, 062119 (2013). https://doi.org/10.1103/PhysRevA.87.062119

2.9. filter_functions.types module

Defines custom types for the package.

2.10. filter_functions.util module

This module provides various helper functions.

2.10.1. Functions

abs2(): Absolute value squared
get_indices_from_identifiers(): The the indices of a subset of identifiers within a list of identifiers.
tensor(): Fast, flexible tensor product of an arbitrary number of inputs using einsum()
tensor_insert(): For an array that is known to be a tensor product, insert arrays at a given position in the product chain
tensor_merge(): For two arrays that are tensor products of known dimensions, merge them at arbitary positions in the product chain
tensor_transpose(): For a tensor product, transpose the order of the constituents in the product chain
mdot(): Multiple matrix product
remove_float_errors(): Set entries whose absolute value is below a certain threshold to zero
oper_equiv(): Determine if two vectors or operators are equal up to a global phase
dot_HS(): Hilbert-Schmidt inner product
get_sample_frequencies(): Get frequencies with typical infrared and ultraviolet cutoffs for a PulseSequence
progressbar(): A progress bar for loops. Uses tqdm if available and a simple custom one if not.
hash_array_along_axis(): Return a list of hashes along a given axis
all_array_equal(): Check if all arrays in an iterable are equal

2.10.2. Exceptions

CalculationError: Exception raised if trying to fetch the pulse correlation function when it was not computed during concatenation

abs2(x: ndarray) → ndarray[source]

Fast function to calculate the absolute value squared,

\[|\cdot|^2 := \Re(\cdot)^2 + \Im(\cdot)^2\]

Equivalent to:

np.abs(x)**2

all_array_equal(it: Iterable) → bool[source]: Return True if all array elements of it are equal by hashing the bytes representation of each array. Note that this is not thread-proof.

dot_HS(U: Union[ndarray, Qobj], V: Union[ndarray, Qobj], eps: Optional[float] = None) → Union[float, complex, ndarray][source]

Return the Hilbert-Schmidt inner product of U and V,

\[\langle U, V\rangle_\mathrm{HS} := \mathrm{tr}(U^\dagger V).\]

Parameters

U, V: qutip.Qobj or ndarray: Objects to compute the inner product of.
eps: float: The floating point precision. The result is rounded to abs(int(np.log10(eps))) decimals if eps > 0.

Returns

result: float, complex: The result rounded to precision eps.

Examples

>>> U, V = paulis[1:3]
>>> dot_HS(U, V)
0.0
>>> dot_HS(U, U)
2.0

get_sample_frequencies(pulse: PulseSequence, n_samples: int = 300, spacing: str = 'log', include_quasistatic: bool = False, omega_min: Optional[float] = None, omega_max: Optional[float] = None) → ndarray[source]

Get n_samples sample frequencies spaced ‘linear’ or ‘log’.

The ultraviolet cutoff is taken to be one order of magnitude larger than the timescale of the pulse tau. In the case of log spacing, the values are clipped in the infrared at two orders of magnitude below the timescale of the pulse.

Parameters

pulse: PulseSequence: The pulse to get frequencies for.
n_samples: int, optional: The number of frequency samples. Default is 300.
spacing: str, optional: The spacing of the frequencies. Either ‘log’ or ‘linear’, default is ‘log’.
include_quasistatic: bool, optional: Include zero frequency. Default is False.
omega_min, omega_max: float, optional: Minimum and maximum angular frequencies included (DC notwithstanding). Default to \(2\pi\times 10^{-2}/\tau\) and \(2\pi\times 10^{+1}/\Delta t_{\mathrm{min}}\).

Returns

omega: ndarray: The angular frequencies.

hash_array_along_axis(arr: ndarray, axis: int = 0) → List[int][source]: Return the hashes of arr along the first axis

mdot(arr: Sequence, axis: int = 0) → ndarray[source]: Multiple matrix products along axis

oper_equiv(psi: Union[ndarray, Qobj], phi: Union[ndarray, Qobj], eps: Optional[float] = None, normalized: bool = False) → Tuple[bool, float][source]

Checks whether psi and phi are equal up to a global phase, i.e.

\[|\psi\rangle = e^{i\chi}|\phi\rangle \Leftrightarrow \langle \phi|\psi\rangle = e^{i\chi},\]

and returns the phase. If the first return value is false, the second is meaningless in this context. psi and phi can also be operators.

Parameters

psi, phi: qutip.Qobj or array_like: Vectors or operators to be compared
eps: float: The tolerance below which the two objects are treated as equal, i.e., the function returns True if abs(1 - modulus) <= eps.
normalized: bool: Flag indicating if psi and phi are normalized with respect to the Hilbert-Schmidt inner product dot_HS().

Examples

>>> psi = paulis[1]
>>> phi = paulis[1]*np.exp(1j*1.2345)
>>> oper_equiv(psi, phi)
(True, 1.2345)

progressbar(iterable: Iterable, *args, **kwargs)[source]

Progress bar for loops. Uses tqdm.

Usage:

for i in progressbar(range(10)):
    do_something()

remove_float_errors(arr: ndarray, eps_scale: Optional[float] = None) → ndarray[source]

Clean up arr by removing floating point numbers smaller than the dtype’s precision multiplied by eps_scale. Treats real and imaginary parts separately.

Obviously only works for arrays with norm ~1.

tensor(*args, rank: int = 2, optimize: Union[bool, str] = False) → ndarray[source]

Fast, flexible tensor product using einsum. The product is taken over the last rank axes and broadcast over the remaining axes which thus need to follow numpy broadcasting rules. Note that vectors are treated as rank 2 tensors with shape (1, x) or (x, 1).

For example, the following shapes are compatible:

rank == 2 (e.g. matrices or vectors):

(a, b, c, d, d), (a, b, c, e, e) -> (a, b, c, d*e, d*e)
(a, b, c), (a, d, e) -> (a, b*d, c*e)
(a, b), (c, d, e) -> (c, a*d, b*e)
(1, a), (b, 1, c) -> (b, 1, a*c)

rank == 1:

(a, b), (a, c) -> (a, b*c)
(a, b, 1), (a, c) -> (a, b, c)

Parameters

args: array_like: The elements of the tensor product
rank: int, optional (default: 2): The rank of the tensors. E.g., for a Kronecker product between two matrices rank == 2. The remaining axes are broadcast over.
optimize: bool|str, optional (default: False): Optimize the tensor contraction order. Passed through to numpy.einsum().

See also

numpy.kron: NumPy tensor product.
tensor_insert: Insert array at given position in tensor product chain.
tensor_merge: Merge tensor product chains.
tensor_transpose: Transpose the order of a tensor product chain.

Examples

>>> Z = np.diag([1, -1])
>>> np.array_equal(tensor(Z, Z), np.kron(Z, Z))
True

>>> A, B = np.arange(2), np.arange(2, 5)
>>> tensor(A, B, rank=1)
array([[0, 0, 0, 2, 3, 4]])

>>> args = np.random.randn(4, 10, 3, 2)
>>> result = tensor(*args, rank=1)
>>> result.shape == (10, 3, 2**4)
True
>>> result = tensor(*args, rank=2)
>>> result.shape == (10, 3**4, 2**4)
True

>>> A, B = np.random.randn(1, 3), np.random.randn(3, 4)
>>> result = tensor(A, B)
>>> result.shape == (1*3, 3*4)
True

>>> A, B = np.random.randn(3, 1, 2), np.random.randn(2, 2, 2)
>>> try:
...     result = tensor(A, B, rank=2)
... except ValueError as err:  # cannot broadcast over axis 0
...     print(err)
Incompatible shapes (3, 1, 2) and (2, 2, 2) for tensor product of rank 2.
>>> result = tensor(A, B, rank=3)
>>> result.shape == (3*2, 1*2, 2*2)
True

tensor_insert(arr: ndarray, *args, pos: Union[int, Sequence[int]], arr_dims: Sequence[Sequence[int]], rank: int = 2, optimize: Union[bool, str] = False) → ndarray[source]

For a tensor product arr, insert args into the product chain at pos. E.g, if \(\verb|arr|\equiv A\otimes B\otimes C\) and \(\verb|pos|\equiv 2\), the result will be the tensor product

\[A\otimes B\otimes\left[\bigotimes_{X\in\verb|args|}X\right] \otimes C.\]

This function works in a similar way to numpy.insert() and the following would be functionally equivalent in the case that the constituent tensors of the product arr are known:

>>> tensor_insert(tensor(*arrs, rank=rank), *args, pos=pos, arr_dims=...,
...               rank=rank)

>>> tensor(*np.insert(arrs, pos, args, axis=0), rank=rank)

Parameters

arr: ndarray

The tensor product in whose chain the other args should be inserted

*args: ndarray

The tensors to be inserted in the product chain

pos: int|sequence of ints

The position(s) at which the args are inserted in the product chain. If an int and len(args) > 1, it is repeated so that all args are inserted in a row. If a sequence, it should indicate the indices in the original tensor product chain that led to arr before which args should be inserted.

arr_dims: array_like, shape (rank, n_const)

The last rank dimensions of the n_const constituent tensors of the tensor product arr as a list of lists with the list at position i containing the i-th relevant dimension of all args. Since the remaing axes are broadcast over, their shape is irrelevant.

For example, if arr = tensor(a, b, c, rank=2) and a,b,c have shapes (2, 3, 4), (5, 2, 2, 1), (2, 2), arr_dims = [[3, 2, 2], [4, 1, 2]].

rank: int, optional (default: 2)

The rank of the tensors. E.g., for a Kronecker product between two vectors, rank == 1, and between two matrices rank == 2. The remaining axes are broadcast over.

optimize: bool|str, optional (default: False)

Optimize the tensor contraction order. Passed through to numpy.einsum().

See also

numpy.insert: NumPy array insertion with similar syntax.
numpy.kron: NumPy tensor product.
tensor_insert: Insert array at given position in tensor product chain.
tensor_merge: Merge tensor product chains.
tensor_transpose: Transpose the order of a tensor product chain.

Examples

>>> I, X, Y, Z = paulis
>>> arr = tensor(X, I)
>>> r = tensor_insert(arr, Y, Z, arr_dims=[[2, 2], [2, 2]], pos=0)
>>> np.allclose(r, tensor(Y, Z, X, I))
True
>>> r = tensor_insert(arr, Y, Z, arr_dims=[[2, 2], [2, 2]], pos=1)
>>> np.allclose(r, tensor(X, Y, Z, I))
True
>>> r = tensor_insert(arr, Y, Z, arr_dims=[[2, 2], [2, 2]], pos=2)
>>> np.allclose(r, tensor(X, I, Y, Z))
True

Other ranks and different dimensions:

>>> from numpy.random import randn
>>> A, B, C = randn(2, 3, 1, 2), randn(2, 2, 2, 2), randn(3, 2, 1)
>>> arr = tensor(A, C, rank=3)
>>> r = tensor_insert(arr, B, pos=1, rank=3,
...                   arr_dims=[[3, 3], [1, 2], [2, 1]])
>>> np.allclose(r, tensor(A, B, C, rank=3))
True

>>> arrs, args = randn(2, 2, 2), randn(2, 2, 2)
>>> arr_dims = [[2, 2], [2, 2]]
>>> r = tensor_insert(tensor(*arrs), *args, pos=(0, 1), arr_dims=arr_dims)
>>> np.allclose(r, tensor(args[0], arrs[0], args[1], arrs[1]))
True
>>> r = tensor_insert(tensor(*arrs), *args, pos=(0, 0), arr_dims=arr_dims)
>>> np.allclose(r, tensor(*args, *arrs))
True
>>> r = tensor_insert(tensor(*arrs), *args, pos=(1, 2), arr_dims=arr_dims)
>>> np.allclose(r, tensor(*np.insert(arrs, (1, 2), args, axis=0)))
True

tensor_merge(arr: ndarray, ins: ndarray, pos: Sequence[int], arr_dims: Sequence[Sequence[int]], ins_dims: Sequence[Sequence[int]], rank: int = 2, optimize: Union[bool, str] = False) → ndarray[source]

For two tensor products arr and ins, merge ins into the product chain at indices pos. E.g, if \(\verb|arr|\equiv A\otimes B\otimes C\), \(\verb|ins|\equiv D\otimes E\), and \(\verb|pos|\equiv [1, 2]\), the result will be the tensor product

\[A\otimes D\otimes B\otimes E\otimes C.\]

This function works in a similar way to numpy.insert() and tensor_insert().

Parameters

arr: ndarray

The tensor product in whose chain the other args should be inserted

ins: ndarray

The tensor product to be inserted in the product chain

pos: sequence of ints

The positions at which the constituent tensors of ins are inserted in the product chain. Should indicate the indices in the original tensor product chain that led to arr before which the constituents of ins should be inserted.

arr_dims: array_like, shape (rank, n_const)

The last rank dimensions of the n_const constituent tensors of the tensor product arr as a list of lists with the list at position i containing the i-th relevant dimension of all args. Since the remaing axes are broadcast over, their shape is irrelevant.

For example, if arr = tensor(a, b, c, rank=2) and a,b,c have shapes (2, 3, 4), (5, 2, 2, 1), (2, 2), arr_dims = [[3, 2, 2], [4, 1, 2]].

ins_dims: array_like, shape (rank, n_const)

The last rank dimensions of the n_const constituent tensors of the tensor product ins as a list of lists with the list at position i containing the i-th relevant dimension of ins. Since the remaing axes are broadcast over, their shape is irrelevant.

rank: int, optional (default: 2)

The rank of the tensors. E.g., for a Kronecker product between two vectors, rank == 1, and between two matrices rank == 2. The remaining axes are broadcast over.

optimize: bool|str, optional (default: False)

Optimize the tensor contraction order. Passed through to numpy.einsum().

See also

numpy.insert: NumPy array insertion with similar syntax.
numpy.kron: NumPy tensor product.
tensor: Fast tensor product with broadcasting.
tensor_insert: Insert array at given position in tensor product chain.
tensor_transpose: Transpose the order of a tensor product chain.

Examples

>>> I, X, Y, Z = paulis
>>> arr = tensor(X, Y, Z)
>>> ins = tensor(I, I)
>>> r1 = tensor_merge(arr, ins, pos=[1, 2], arr_dims=[[2]*3, [2]*3],
...                   ins_dims=[[2]*2, [2]*2])
>>> np.allclose(r1, tensor(X, I, Y, I, Z))
True
>>> r2 = tensor_merge(ins, arr, pos=[0, 1, 2], arr_dims=[[2]*2, [2]*2],
...                   ins_dims=[[2]*3, [2]*3])
>>> np.allclose(r1, r2)
True

tensor_insert() can provide the same functionality in some cases:

>>> arr = tensor(Y, Z)
>>> ins = tensor(I, X)
>>> r1 = tensor_merge(arr, ins, pos=[0, 0], arr_dims=[[2]*2, [2]*2],
...                   ins_dims=[[2]*2, [2]*2])
>>> r2 = tensor_insert(arr, I, X, pos=[0, 0], arr_dims=[[2]*2, [2]*2])
>>> np.allclose(r1, r2)
True

Also tensors of rank other than 2 and numpy broadcasting are supported:

>>> arr = np.random.randn(2, 10, 3, 4)
>>> ins = np.random.randn(2, 10, 3, 2)
>>> r = tensor_merge(tensor(*arr, rank=1), tensor(*ins, rank=1), [0, 1],
...                  arr_dims=[[4, 4]], ins_dims=[[2, 2]], rank=1)
>>> np.allclose(r, tensor(ins[0], arr[0], ins[1], arr[1], rank=1))
True

tensor_transpose(arr: ndarray, order: Sequence[int], arr_dims: Sequence[Sequence[int]], rank: int = 2) → ndarray[source]

Transpose the order of a tensor product chain.

Parameters

arr: ndarray

The tensor product whose chain should be reordered.

order: sequence of ints

The transposition order. If arr == tensor(A, B) and order == (1, 0), the result will be tensor(B, A).

arr_dims: array_like, shape (rank, n_const)

The last rank dimensions of the n_const constituent tensors of the tensor product arr as a list of lists with the list at position i containing the i-th relevant dimension of all args. Since the remaing axes are broadcast over, their shape is irrelevant.

For example, if arr = tensor(a, b, c, rank=2) and a,b,c have shapes (2, 3, 4), (5, 2, 2, 1), (2, 2), arr_dims = [[3, 2, 2], [4, 1, 2]].

rank: int, optional (default: 2)

The rank of the tensors. E.g., for a Kronecker product between two vectors, rank == 1, and between two matrices rank == 2. The remaining axes are broadcast over.

Returns

transposed_arr: ndarray: The tensor product arr with its order transposed according to order

See also

numpy.insert: NumPy array insertion with similar syntax.
numpy.kron: NumPy tensor product.
tensor: Fast tensor product with broadcasting.
tensor_insert: Insert array at given position in tensor product chain.
tensor_merge: Merge tensor product chains.

Examples

>>> I, X, Y, Z = paulis
>>> arr = tensor(X, Y, Z)
>>> transposed = tensor_transpose(arr, [1, 2, 0], arr_dims=[[2, 2, 2]]*2)
>>> np.allclose(transposed, tensor(Y, Z, X))
True

2.11. Module contents

Package for efficient calculation of generalized filter functions

class Basis(basis_array: Sequence, traceless: Optional[bool] = None, btype: Optional[str] = None, labels: Optional[Sequence[str]] = None)[source]

Bases: ndarray

Class for operator bases. There are several ways to instantiate a Basis object:

by just calling this constructor with a (possibly incomplete) array of basis matrices. No checks regarding orthonormality or hermiticity are performed.

by calling one of the classes alternative constructors (classmethods):

pauli(): Pauli operator basis

ggm(): Generalized Gell-Mann basis

from_partial() Generate an complete basis from partial elements

These bases guarantee the following properties:

hermitian

orthonormal

[traceless] (can be controlled by a flag)

Since Basis is a subclass of NumPy’s ndarray, it inherits all of its attributes, e.g. shape. The following attributes behave slightly differently to a ndarray, however

A == B is True if all elements evaluate almost equal, i.e. equivalent to np.allclose(A, B).

basis.T transposes the last two axes of basis. For a full basis, this corresponds to transposing each element individually. For a basis element, it corresponds to normal transposition.

Parameters

basis_array: array_like, shape (n, d, d): An array or list of square matrices that are elements of an operator basis spanning \(\mathbb{C}^{d\times d}\). n should be smaller than or equal to d**2.
traceless: bool, optional (default: auto): Controls whether a traceless basis is forced. Here, traceless means that the first element of the basis is the identity and the remaining elements are matrices of trace zero. If an element of basis_array is neither traceless nor the identity and traceless == True, an exception will be raised. Defaults to True if basis_array is traceless and False if not.
btype: str, optional (default: ``’custom’``): A string describing the basis type. For example, a basis created by the factory method pauli() has btype ‘pauli’.
labels: sequence of str, optional: A list of labels for the individual basis elements. Defaults to ‘C_0’, ‘C_1’, …

Attributes

Other than the attributes inherited from ``ndarray``, a ``Basis``
instance has the following attributes:
btype: str: Basis type.
labels: sequence of str: The labels for the basis elements.
d: int: Dimension of the space spanned by the basis.
H: Basis: Hermitian conjugate.
isherm: bool: If the basis is hermitian.
isorthonorm: bool: If the basis is orthonormal.
istraceless: bool: If the basis is traceless except for an identity element
iscomplete: bool: If the basis is complete, ie spans the full space.
sparse: COO, shape (n, d, d): Representation in the COO format supplied by the sparse package.
four_element_traces: COO, shape (n, n, n, n): Traces over all possible combinations of four elements of self. This is required for the calculation of the error transfer matrix and thus cached in the Basis instance.
Most of the attributes above are properties which are lazily
evaluated and cached.

Methods

Other than the methods inherited from ``ndarray``, a ``Basis``
instance has the following methods:
normalize(b)	Normalizes the basis (used internally when creating a basis from elements)
tidyup(eps_scale=None)	Cleans up floating point errors in-place to make zeros actual zeros. `eps_scale` is an optional argument multiplied to the data type’s `eps` to get the absolute tolerance.

Constructor.

property H: Basis: Return the basis hermitian conjugated element-wise.

property T: Basis: Return the basis transposed element-wise.

property four_element_traces: COO: Return all traces of the form \(\mathrm{tr}(C_i C_j C_k C_l)\) as a sparse COO array for \(i,j,k,l > 0\) (i.e. excluding the identity).

classmethod from_partial(partial_basis_array: Sequence, traceless: Optional[bool] = None, btype: Optional[str] = None, labels: Optional[Sequence[str]] = None) → Basis[source]

Generate complete and orthonormal basis from a partial set.

The basis is completed using singular value decomposition to determine the null space of the expansion coefficients of the partial basis with respect to another complete basis.

Parameters

partial_basis_array: array_like: A sequence of basis elements.
traceless: bool, optional: If a traceless basis should be generated (i.e. the first element is the identity and all the others have trace zero).
btype: str, optional: A custom identifier.
labels: Sequence[str], optional: A list of custom labels for each element. If len(labels) == len(partial_basis_array), the newly created elements get labels ‘C_i’.

Returns

basis: Basis, shape (d**2, d, d): The orthonormal basis.

Raises

ValueError: If the given elements are not orthonormal.
ValueError: If the given elements are not traceless but traceless==True.
ValueError: If not len(partial_basis_array) or d**2 labels were given.

classmethod ggm(d: int) → Basis[source]

Returns a generalized Gell-Mann basis in \(d\) dimensions [Bert08] where the elements \(\Lambda_i\) are normalized with respect to the Hilbert-Schmidt inner product,

\[\begin{split}\langle\Lambda_i,\Lambda_j\rangle &= \mathrm{Tr}\,\Lambda_i^\dagger\Lambda_j \\ &= \delta_{ij}.\end{split}\]

Parameters

d: int: The dimensionality of the space spanned by the basis

Returns

basis: Basis: The Basis object representing the GGM.

References

Bert08: Bertlmann, R. A., & Krammer, P. (2008). Bloch vectors for qudits. Journal of Physics A: Mathematical and Theoretical, 41(23). https://doi.org/10.1088/1751-8113/41/23/235303

property iscomplete: bool: Returns True if basis is complete.

property isherm: bool: Returns True if all basis elements are hermitian.

property isorthonorm: bool: Returns True if basis is orthonormal.

property istraceless: bool: Returns True if basis is traceless except for possibly the identity.

normalize(copy: bool = False) → Union[None, Basis][source]: Normalize the basis.

classmethod pauli(n: int) → Basis[source]

Returns a Pauli basis for \(n\) qubits, i.e. the basis spans the space \(\mathbb{C}^{d\times d}\) with \(d = 2^n\):

\[\mathcal{P} = \{I, X, Y, Z\}^{\otimes n}.\]

The elements \(\sigma_i\) are normalized with respect to the Hilbert-Schmidt inner product,

\[\begin{split}\langle\sigma_i,\sigma_j\rangle &= \mathrm{Tr}\,\sigma_i^\dagger\sigma_j \\ &= \delta_{ij}.\end{split}\]

Parameters

n: int: The number of qubits.

Returns

basis: Basis: The Basis object representing the Pauli basis.

property sparse: COO: Return the basis as a sparse COO array

tidyup(eps_scale: Optional[float] = None) → None[source]: Wraps util.remove_float_errors.

class PulseSequence(*args, **kwargs)[source]

Bases: object

A class for pulse sequences and their filter functions.

The Hamiltonian is separated into a control and a noise part with

\[\begin{split}\mathcal{H}_c &= \sum_i a_i(t) A_i \\ \mathcal{H}_n &= \sum_j s_j(t) b_j(t) B_j\end{split}\]

where \(A_i\) and \(B_j\) are hermitian operators and \(b_j(t)\) are classically fluctuating noise variables captured in a power spectral density and not needed at instantiation of a PulseSequence.

Parameters

H_c: list of lists

A nested list of n_cops nested lists as taken by QuTiP functions (see for example qutip.propagator.propagator()) describing the control part of the Hamiltonian. The i-th entry of the list should be a list consisting of the i-th operator \(A_i\) making up the control Hamiltonian and a list or array \(a_i(t)\) describing the magnitude of that operator during the time intervals dt. Optionally, the list may also include operator identifiers. That is, H_c should look something like this:

H = [[c_oper1, c_coeff1, c_oper_identifier1],
     [c_oper2, c_coeff2, c_oper_identifier2], ...]

The operators may be given either as NumPy arrays or QuTiP Qobjs and each coefficient array should have the same number of elements as dt, and should be given in units of \(\hbar\). If not every sublist (read: operator) was given a identifier, they are automatically filled up with ‘A_i’ where i is the position of the operator.

H_n: list of lists

A nested list of n_nops nested lists as taken by QuTiP functions (see for example qutip.propagator.propagator()) describing the noise part of the Hamiltonian. The j-th entry of the list should be a list consisting of the j-th operator \(B_j\) making up the noise Hamiltonian and a list or array describing the sensitivity \(s_j(t)\) of the system to the noise operator during the time intervals dt. Optionally, the list may also include operator identifiers. That is, H_n should look something like this:

H = [[n_oper1, n_coeff1, n_oper_identifier1],
     [n_oper2, n_coeff2, n_oper_identifier2], ...]

The operators may be given either as NumPy arrays or QuTiP Qobjs and each coefficient array should have the same number of elements as dt, and should be given in units of \(\hbar\). If not every sublist (read: operator) was given a identifier, they are automatically filled up with ‘B_i’ where i is the position of the operator.

dt: array_like, shape (n_dt,)

The segment durations of the Hamiltonian (i.e. durations of constant control). Internally, the control operation is taken to start at \(t_0\equiv 0\), i.e. the edges of the constant control segments are at times t = [0, *np.cumsum(dt)].

basis: Basis, shape (d**2, d, d), optional

The operator basis in which to calculate. If a Generalized Gell-Mann basis (see ggm()) is chosen, some calculations will be faster for large dimensions due to a simpler basis expansion. However, when extending the pulse sequence to larger qubit registers, cached filter functions cannot be retained since the GGM basis does not factor into tensor products. In this case a Pauli basis is preferable.

Notes

Due to the heavy use of NumPy’s einsum() function, results have a floating point error of ~1e-13.

Examples

A rotation by \(\pi\) around the axis between x and y preceeded and followed by a period of free evolution with the system subject to dephasing noise.

>>> import qutip as qt; import numpy as np
>>> H_c = [[qt.sigmax(), [0, np.pi, 0]],
           [qt.sigmay(), [0, np.pi, 0]]]
>>> # Equivalent pulse:
>>> # H_c = [[qt.sigmax() + qt.sigmay(), [0, np.pi, 0]]]
>>> # The noise sensitivity is constant
>>> H_n = [[qt.sigmaz()/np.sqrt(2), [1, 1, 1], 'Z']]
>>> dt = [1, 1, 1]
>>> # Free evolution between t=0 and t=1, rotation between t=1 and t=2,
>>> # and free evolution again from t=2 to t=3.
>>> pulse = PulseSequence(H_c, H_n, dt)
>>> pulse.c_oper_identifiers
['A_0', 'A_1']
>>> pulse.n_oper_identifiers
['Z']
>>> omega = np.logspace(-1, 2, 500)
>>> F = pulse.get_filter_function(omega)    # shape (1, 500)
>>> # Plot the resulting filter function:
>>> from filter_functions import plotting
>>> fig, ax, leg = plotting.plot_filter_function(pulse)

Attributes

c_opers: ndarray, shape (n_cops, d, d): Control operators. Note that they are stored sorted by their corresponding identifiers.
n_opers: ndarray, shape (n_nops, d, d): Noise operators. Note that they are stored sorted by their corresponding identifiers.
c_oper_identifers: sequence of str: Identifiers for the control operators of the system. Stored sorted.
n_oper_identifers: sequence of str: Identifiers for the noise operators of the system. Stored sorted.
c_coeffs: ndarray, shape (n_cops, n_dt): Control parameters in units of \(\hbar\). Note that they are stored sorted by their corresponding identifiers.
n_coeffs: ndarray, shape (n_nops, n_dt): Noise sensitivities in units of \(\hbar\). Note that they are stored sorted by their corresponding identifiers.
dt: ndarray, shape (n_dt,): Time steps
t: ndarray, shape (n_dt + 1,): Absolute times taken to start at \(t_0\equiv 0\)
tau: float: Total duration. Equal to t[-1].
d: int: Dimension of the Hamiltonian
basis: Basis, shape (d**2, d, d): The operator basis used for calculation
nbytes: int: An estimate of the memory consumed by the PulseSequence instance and its attributes
If the Hamiltonian was diagonalized, the eigenvalues and -vectors as
well as the cumulative propagators are cached:
eigvals: ndarray, shape (n_dt, d): Eigenvalues \(D^{(g)}\)
eigvecs: ndarray, shape (n_dt, d, d): Eigenvectors \(V^{(g)}\)
propagators: ndarray, shape (n_dt+1, d, d): Cumulative propagators \(Q_g\)
total_propagator: ndarray, shape (d, d): The total propagator \(Q\) of the pulse alone. That is, \(|\psi(\tau)\rangle = Q|\psi(0)\rangle\).
total_propagator_liouville: array_like, shape (d**2, d**2): The transfer matrix for the total propagator of the pulse. Given by liouville_representation(pulse.total_propagator, pulse.basis).
Furthermore, when the filter function is calculated, the frequencies
are cached as well as other relevant quantities.

Methods

get_filter_function_derivative(omega[, ...])

Calculate the pulse sequence's filter function derivative.

cleanup(method=’conservative’)	Delete cached attributes
is_cached(attr)	Checks if a given attribute of the `PulseSequence` is cached
diagonalize()	Diagonalize the Hamiltonian of the pulse sequence, computing eigenvalues and -vectors as well as cumulative propagators
get_control_matrix(omega, show_progressbar=False)	Calculate the control matrix for frequencies omega
get_filter_function(omega, which=’fidelity’, show_progressbar=False)	Calculate the filter function for frequencies omega
get_pulse_correlation_filter_function(which=’fidelity’)	Get the pulse correlation filter function (only possible if computed during concatenation)
propagator_at_arb_t(t)	Calculate the propagator at arbitrary times

Initialize a PulseSequence instance.

cache_control_matrix(omega: Sequence[float], control_matrix: Optional[ndarray] = None, show_progressbar: bool = False, cache_intermediates: bool = False) → None[source]

Cache the control matrix and the frequencies it was calculated for.

Parameters

omega: array_like, shape (n_omega,): The frequencies for which to cache the filter function.
control_matrix: array_like, shape ([n_nops,] n_nops, d**2, n_omega), optional: The control matrix for the frequencies omega. If None, it is computed.
show_progressbar: bool: Show a progress bar for the calculation of the control matrix.
cache_intermediates: bool, optional: Keep intermediate terms of the calculation that are also required by other computations. Only applies if control_matrix is not supplied.

cache_filter_function(omega: Sequence[float], control_matrix: Optional[ndarray] = None, filter_function: Optional[ndarray] = None, which: str = 'fidelity', order: int = 1, show_progressbar: bool = False, cache_intermediates: bool = False) → None[source]

Cache the filter function. If control_matrix.ndim == 4, it is taken to be the ‘pulse correlation control matrix’ and summed along the first axis. In that case, also the pulse correlation filter function is calculated and cached. Total phase factors and transfer matrices of the the cumulative propagator are also cached so they can be reused during concatenation.

Parameters

omega: array_like, shape (n_omega,): The frequencies for which to cache the filter function.
control_matrix: array_like, shape ([n_nops,] n_nops, d**2, n_omega), optional: The control matrix for the frequencies omega. If None, it is computed and the filter function derived from it.
filter_function: array_like, shape (n_nops, n_nops, [d**2, d**2,] n_omega), optional: The filter function for the frequencies omega. If None, it is computed from R in the case order == 1 and from scratch else.
which: str, optional: Which filter function to cache. Either ‘fidelity’ (default) or ‘generalized’.
order: int, optional: First or second order filter function.
show_progressbar: bool, optional: Show a progress bar for the calculation of the control matrix.
cache_intermediates: bool, optional: Keep intermediate terms of the calculation that are also required by other computations.

See also

PulseSequence.get_filter_function: Getter method

cache_total_phases(omega: Sequence[float], total_phases: Optional[ndarray] = None) → None[source]

Cache the total phase factors for this pulse and omega.

Parameters

omega: array_like, shape (n_omega,): The frequencies for which to cache the phase factors.
total_phases: array_like, shape (n_omega,), optional: The total phase factors for the frequencies omega. If None, they are computed.

cleanup(method: str = 'conservative') → None[source]

Delete cached byproducts of the calculation of the filter function that are not necessarily needed anymore in order to free up memory.

Parameters

method: optional

If set to ‘conservative’ (the default), only the following attributes are deleted:

_eigvals

_eigvecs

_propagators

If set to ‘greedy’, all of the above as well as the following attributes are deleted:

_total_propagator

_total_propagator_liouville

_total_phases

_control_matrix

_control_matrix_pc

If set to ‘all’, all of the above as well as the following attributes are deleted:

omega

_filter_function

_filter_function_gen

_filter_function_pc

_filter_function_pc_gen

_filter_function_2

_intermediates[‘control_matrix_step’]

If set to ‘frequency dependent’ only attributes that are functions of frequency are initalized to None.

Note that if this PulseSequence is concatenated with another one, some of the attributes might need to be calculated again, resulting in slower execution of the concatenation.

diagonalize() → None[source]: Diagonalize the Hamiltonian defining the pulse sequence.

property duration: float: The duration of the pulse. Alias of tau.

property eigvals: ndarray: Get the eigenvalues of the pulse’s Hamiltonian.

property eigvecs: ndarray: Get the eigenvectors of the pulse’s Hamiltonian.

get_control_matrix(omega: Sequence[float], show_progressbar: bool = False, cache_intermediates: bool = False) → ndarray[source]

Get the control matrix for the frequencies omega. If it has been cached for the same frequencies, the cached version is returned, otherwise it is calculated from scratch.

Parameters

omega: array_like, shape (n_omega,): The frequencies at which to evaluate the control matrix.
show_progressbar: bool: Show a progress bar for the calculation of the control matrix.
cache_intermediates: bool, optional: Keep intermediate terms of the calculation that are also required by other computations.

Returns

control_matrix: ndarray, shape (n_nops, d**2, n_omega): The control matrix for the noise operators.

get_filter_function(omega: Sequence[float], which: str = 'fidelity', order: int = 1, show_progressbar: bool = False, cache_intermediates: bool = False) → ndarray[source]

Get the first or second order filter function.

The filter function is cached so it doesn’t need to be calculated twice for the same frequencies.

Parameters

omega: array_like, shape (n_omega,): The frequencies at which to evaluate the filter function.
which: str, optional: Which filter function to return. Either ‘fidelity’ (default) or ‘generalized’ (see Notes). Only if order == 1.
order: int, optional: First or second order filter function.
show_progressbar: bool, optional: Show a progress bar for the calculation of the control matrix.
cache_intermediates: bool, optional: Keep intermediate terms of the calculation that are also required by other computations.

Returns

filter_function: ndarray, shape (n_nops, n_nops, [d**2, d**2,] n_omega): The filter function for each combination of noise operators as a function of omega.

Notes

The first-order generalized filter function is given by

\[F_{\alpha\beta,kl}(\omega) = \tilde{\mathcal{B}}_{\alpha k}^\ast(\omega) \tilde{\mathcal{B}}_{\beta l}(\omega),\]

where \(\alpha,\beta\) are indices counting the noise operators \(B_\alpha\) and \(k,l\) indices counting the basis elements \(C_k\).

The fidelity filter function is obtained by tracing over the basis indices:

\[F_{\alpha\beta}(\omega) = \sum_{k} F_{\alpha\beta,kk}(\omega).\]

get_filter_function_derivative(omega: Sequence[float], control_identifiers: Optional[Sequence[str]] = None, n_oper_identifiers: Optional[Sequence[str]] = None, n_coeffs_deriv: Optional[Sequence[Sequence[float]]] = None) → ndarray[source]

Calculate the pulse sequence’s filter function derivative.

Parameters

omega: array_like, shape (n_omega,): Frequencies at which the pulse control matrix is to be evaluated.
control_identifiers: Sequence[str], shape (n_ctrl,): Sequence of strings with the control identifiers to distinguish between control and drift Hamiltonian. The default is None, in which case the derivative is computed for all known non-noise operators.
n_oper_identifiers: Sequence[str], shape (n_nops,): Sequence of strings with the noise identifiers for which to compute the derivative contribution. The default is None, in which case it is computed for all known noise operators.
n_coeffs_deriv: array_like, shape (n_nops, n_ctrl, n_dt): The derivatives of the noise susceptibilities by the control amplitudes. The rows and columns should be in the same order as the corresponding identifiers above. Defaults to None, in which case the coefficients are assumed to be constant and hence their derivative vanishing.

Warning

Internally, control and noise terms of the Hamiltonian are stored alphanumerically sorted by their identifiers. If the noise and/or control identifiers above are not explicitly given, the rows and/or columns of this parameter need to be sorted in the same fashion.

Returns

filter_function_deriv: ndarray, shape (n_nops, n_t, n_ctrl, n_omega): The regular filter functions’ derivatives for variation in each control contribution. Sorted in the same fashion as n_coeffs_deriv or, if not given, alphanumerically by the identifiers.

get_pulse_correlation_control_matrix() → ndarray[source]: Get the pulse correlation control matrix if it was cached.

get_pulse_correlation_filter_function(which: str = 'fidelity') → ndarray[source]

Get the pulse correlation filter function given by

\[F_{\alpha\beta}^{(gg')}(\omega) = e^{i\omega(t_{g-1} - t_{g'-1})} \tilde{\mathcal{B}}^{(g)}(\omega)\mathcal{Q}^{(g-1)} \mathcal{Q}^{(g'-1)\dagger} \tilde{\mathcal{B}}^{(g')\dagger}(\omega),\]

where \(g,g'\) index the pulse in the sequence and \(\alpha,\beta\) index the noise operators, if it was computed during concatenation. Since the calculation requires the individual pulse’s control matrices and phase factors, which are not retained after concatenation, the pulse correlation filter function cannot be computed afterwards.

Note that the frequencies for which the filter function was calculated are not stored.

Returns

filter_function_pc: ndarray, shape (n_pls, n_pls, n_nops, n_nops, n_omega)

The pulse correlation filter function for each noise operator as a function of omega. The first two axes correspond to the pulses in the sequence, i.e. if the concatenated pulse sequence is \(C\circ B\circ A\), the first two axes are arranged like

\[\begin{split}F_{\alpha\beta}^{(gg')} &= \begin{pmatrix} F_{\alpha\beta}^{(AA)} & F_{\alpha\beta}^{(AB)} & F_{\alpha\beta}^{(AC)} \\ F_{\alpha\beta}^{(BA)} & F_{\alpha\beta}^{(BB)} & F_{\alpha\beta}^{(BC)} \\ F_{\alpha\beta}^{(CA)} & F_{\alpha\beta}^{(CB)} & F_{\alpha\beta}^{(CC)} \end{pmatrix}\end{split}\]

for \(g,g'\in\{A, B, C\}\).

get_total_phases(omega: Sequence[float]) → ndarray[source]: Get the (cached) total phase factors for this pulse and omega.

is_cached(attr: str) → bool[source]: Returns True if the attribute is cached

property nbytes: int: Return an estimate of the amount of memory consumed by this object (or, more precisely, the array attributes of this object).

property omega: ndarray: Cached frequencies

propagator_at_arb_t(t: Sequence[float]) → ndarray[source]: Calculate the cumulative propagator Q(t) at times t by making use of the fact that we assume piecewise-constant control.

property propagators: ndarray: Get the eigenvectors of the pulse’s Hamiltonian.

property t: ndarray: The times of the pulse.

property tau: Union[float, int]: The duration of the pulse.

property total_propagator: ndarray: Get total propagator of the pulse.

property total_propagator_liouville: ndarray: Get the transfer matrix for the total propagator of the pulse.

concatenate(pulses: Iterable[PulseSequence], calc_pulse_correlation_FF: bool = False, calc_filter_function: Optional[bool] = None, which: str = 'fidelity', omega: Optional[Sequence[float]] = None, show_progressbar: bool = False) → PulseSequence[source]

Concatenate an arbitrary number of pulses. Note that pulses are concatenated left-to-right, that is,

\[\mathtt{concatenate((A, B))} \equiv B \circ A\]

so that \(A\) is executed before \(B\) when applying the concatenated pulse.

Parameters

pulses: sequence of PulseSequences: The PulseSequence instances to be concatenated. If any of the instances have a cached filter function, the filter function for the composite pulse will also be calculated in order to make use of the speedup gained from concatenating the filter functions. If omega is given, calculation of the composite filter function is forced.
calc_pulse_correlation_FF: bool, optional: Switch to control whether the pulse correlation filter function (see PulseSequence.get_pulse_correlation_filter_function()) is calculated. If omega is not given, the cached frequencies of all pulses need to be equal.
calc_filter_function: bool, optional: Switch to force the calculation of the filter function to be carried out or not. Overrides the automatic behavior of calculating it if at least one pulse has a cached control matrix. If True and no pulse has a cached control matrix, a list of frequencies must be supplied as omega.
which: str, optional: Which filter function to compute. Either ‘fidelity’ (default) or ‘generalized’ (see PulseSequence.get_filter_function() and PulseSequence.get_pulse_correlation_filter_function()).
omega: array_like, optional: Frequencies at which to evaluate the (pulse correlation) filter functions. If None, an attempt is made to use cached frequencies.
show_progressbar: bool: Show a progress bar for the calculation of the control matrix.

Returns

pulse: PulseSequence: The concatenated pulse.

concatenate_periodic(pulse: PulseSequence, repeats: int) → PulseSequence[source]

Concatenate a pulse sequence pulse whose Hamiltonian is periodic repeats times. Although performing the same task, this function is much faster for concatenating many identical pulses with filter functions than concatenate().

Note that for large dimensions, the calculation of the control matrix using this function might be very memory intensive.

Parameters

pulse: PulseSequence: The PulseSequence instance to be repeated. If it has a cached filter function, the filter function for the new pulse will also be computed.
repeats: int: The number of repetitions

Returns

newpulse: PulseSequence: The concatenated PulseSequence

See also

concatenate: Concatenate arbitrary PulseSequences.

Notes

The total control matrix is given by

\[\begin{split}\tilde{\mathcal{B}}(\omega) &= \tilde{\mathcal{B}}^{(1)}(\omega) \sum_{g=0}^{G-1} \left(e^{i\omega T}\right)^g \\ &= \tilde{\mathcal{B}}^{(1)}(\omega)\bigl( \mathbb{I} - e^{i\omega T} \mathcal{Q}^{(1)}\bigr)^{-1}\bigl( \mathbb{I} - \bigl(e^{i\omega T} \mathcal{Q}^{(1)}\bigr)^G\bigr).\end{split}\]

with \(T\) the period of the control Hamiltonian and \(G\) the number of periods. The last equality is valid only if \(\mathbb{I} - e^{i\omega T}\mathcal{Q}^{(1)}\) is invertible.

error_transfer_matrix(pulse: Optional[PulseSequence] = None, spectrum: Optional[ndarray] = None, omega: Optional[Sequence[float]] = None, n_oper_identifiers: Optional[Sequence[str]] = None, second_order: bool = False, cumulant_function: Optional[ndarray] = None, show_progressbar: bool = False, memory_parsimonious: bool = False, cache_intermediates: bool = False) → ndarray[source]

Compute the error transfer matrix up to unitary rotations.

Parameters

pulse: PulseSequence: The PulseSequence instance for which to compute the error transfer matrix.
spectrum: array_like, shape ([[n_nops,] n_nops,] n_omega): The two-sided noise power spectral density in units of inverse frequencies as an array of shape (n_omega,), (n_nops, n_omega), or (n_nops, n_nops, n_omega). In the first case, the same spectrum is taken for all noise operators, in the second, it is assumed that there are no correlations between different noise sources and thus there is one spectrum for each noise operator. In the third and most general case, there may be a spectrum for each pair of noise operators corresponding to the correlations between them. n_nops is the number of noise operators considered and should be equal to len(n_oper_identifiers).
omega: array_like, shape (n_omega,): The frequencies at which to calculate the filter functions.
n_oper_identifiers: array_like, optional: The identifiers of the noise operators for which to evaluate the error transfer matrix. The default is all. Note that, since in general contributions from different noise operators won’t commute, not selecting all noise operators results in neglecting terms of order \(\xi^4\).
second_order: bool, optional: Also take into account the frequency shifts \(\Delta\) that correspond to second order Magnus expansion and constitute unitary terms. Default False.
cumulant_function: ndarray, shape ([[n_pls, n_pls,] n_nops,] n_nops, d**2, d**2): A precomputed cumulant function. If given, pulse, spectrum, omega are not required.
show_progressbar: bool, optional: Show a progress bar for the calculation of the control matrix.
memory_parsimonious: bool, optional: Trade memory footprint for performance. See calculate_decay_amplitudes(). The default is False.
cache_intermediates: bool, optional: Keep and return intermediate terms of the calculation of the control matrix (if it is not already cached) that can be reused for second order or gradients. Can consume large amount of memory, but speed up the calculation.

Returns

error_transfer_matrix: ndarray, shape (d**2, d**2): The error transfer matrix. The individual noise operator contributions are summed up before exponentiating as they might not commute.

See also

calculate_cumulant_function: Calculate the cumulant function \(\mathcal{K}\)
calculate_decay_amplitudes: Calculate the \(\Gamma_{\alpha\beta,kl}\)
infidelity: Calculate only infidelity of a pulse.

Notes

The error transfer matrix is given by

\[\tilde{\mathcal{U}} = \exp\mathcal{K}(\tau)\]

with \(\mathcal{K}(\tau)\) the cumulant function (see calculate_cumulant_function()). For Gaussian noise this expression is exact when taking into account the decay amplitudes \(\Gamma\) and frequency shifts \(\Delta\). As the latter effects coherent errors it can be neglected if we assume that the experimenter has calibrated their pulse.

For non-Gaussian noise the expression above is perturbative and includes noise up to order \(\xi^2\) and hence \(\tilde{\mathcal{U}} = \mathbb{1} + \mathcal{K}(\tau) + \mathcal{O}(\xi^4)\) (although it is evaluated as a matrix exponential in any case).

Given the above expression of the error transfer matrix, the entanglement fidelity is given by

\[\mathcal{F}_\mathrm{e} = \frac{1}{d^2}\mathrm{tr}\,\tilde{\mathcal{U}}.\]

extend(pulse_to_qubit_mapping: Sequence[Sequence[Optional[Union[PulseSequence, Sequence[int], int, Mapping[str, str]]]]], N: Optional[int] = None, d_per_qubit: int = 2, additional_noise_Hamiltonian: Optional[Sequence[Sequence[Union[ndarray, Qobj, Sequence[float]]]]] = None, cache_diagonalization: Optional[bool] = None, cache_filter_function: Optional[bool] = None, omega: Optional[Sequence[float]] = None, show_progressbar: bool = False) → PulseSequence[source]

Map one or more pulse sequences to different qubits.

Parameters

pulse_to_qubit_mapping: sequence of mapping tuples

A sequence of tuples with the first entry a PulseSequence instance and the second an int or tuple of ints indicating the qubits that the PulseSequence should be mapped to. A mapping of operator identifiers may optionally be given as a third element of each tuple. By default, the index of the qubit the operator is mapped to is appended to its identifier.

Pulse sequences defined for multiple qubits may also be extended to non-neighboring qubits. Note that for multi-qubit pulses the order of the qubits is respected, i.e. mapping a pulse to (1, 0) is different from mapping it to (0, 1).

N: int

The total number of qubits the new PulseSequence should be defined for. By default, this is inferred from pulse_to_qubit_mapping.

d_per_qubit: int

The size of the Hilbert space a single qubit requires.

additional_noise_Hamiltonian: list of lists

Additional noise operators and corresponding sensitivities for the new pulse sequence.

cache_diagonalization: bool

Force diagonalizing the new pulse sequence. By default, diagonalization is cached if all pulses in pulse_to_qubit_mapping have been diagonalized since it is much cheaper to get the relevant quantities as tensor products from the mapped pulses instead of diagonalizing the new pulse.

cache_filter_function: bool

Force computing the filter functions for the new pulse sequence. Noise operators of individual pulses will be extended to the new Hilbert space. By default, this is done if all pulses in pulse_to_qubit_mapping have their filter functions cached.

Note that extending the filter functions is only possible if they the mapped pulses are using a separable basis like the Pauli basis.

omega: array_like

Frequencies for which to compute the filter functions if cache_filter_function == True. Defaults to None, in which case the cached frequencies of the individual pulses need to be the same.

show_progressbar: bool

Show a progress bar for the calculation of the control matrix.

Returns

newpulse: PulseSequence: The new pulse sequence on the larger qubit register. The noise operators (and possibly filter functions) are stored in the following order: first those of the multi-qubit pulses in the order they appeared in pulse_to_qubit_mapping, then those of the single-qubit pulses, and lastly any additional ones that may be given by additional_noise_Hamiltonian.

See also

remap: Map PulseSequence to a different qubit.
concatenate: Concatenate PulseSequences (in time).
concatenate_periodic: Periodically concatenate a PulseSequence.

Examples

>>> import filter_functions as ff
>>> I, X, Y, Z = ff.util.paulis
>>> X_pulse = ff.PulseSequence([[X, [np.pi/2], 'X']],
...                            [[X, [1], 'X'], [Z, [1], 'Z']],
...                            [1], basis=ff.Basis.pauli(1))
>>> XX_pulse = ff.extend([(X_pulse, 0), (X_pulse, 1)])
>>> XX_pulse.d
4
>>> XIX_pulse_1 = ff.extend([(X_pulse, 0), (X_pulse, 2)])
>>> XIX_pulse_1.d
8
>>> XXI_pulse = ff.extend([(X_pulse, 0), (X_pulse, 1)], N=3)
>>> XXI_pulse.d
8

Filter functions are automatically cached if they are for mapped pulses:

>>> omega = ff.util.get_sample_frequencies(X_pulse)
>>> X_pulse.cache_filter_function(omega)
>>> XX_pulse = ff.extend([(X_pulse, 0), (X_pulse, 1)])
>>> XX_pulse.is_cached('filter_function')
True

This behavior can also be overriden manually:

>>> XX_pulse = ff.extend([(X_pulse, 0), (X_pulse, 1)],
...                      cache_filter_function=False)
>>> XX_pulse.is_cached('filter_function')
False

Mapping pulses to non-neighboring qubits is also possible:

>>> Y_pulse = ff.PulseSequence([[Y, [np.pi/2], 'Y']],
...                            [[Y, [1], 'Y'], [Z, [1], 'Z']],
...                            [1], basis=ff.Basis.pauli(1))
>>> XXY_pulse = ff.extend([(XX_pulse, (0, 1)), (Y_pulse, 2)])
>>> XYX_pulse = ff.extend([(XX_pulse, (0, 2)), (Y_pulse, 1)])

Additionally, pulses can have the order of the qubits they are defined for permuted (see remap()):

>>> Z_pulse = ff.PulseSequence([[Z, [np.pi/2], 'Z']], [[Z, [1], 'Z']],
...                            [1], basis=ff.Basis.pauli(1))
>>> XY_pulse = ff.extend([(X_pulse, 0), (Y_pulse, 1)])
>>> YZX_pulse = ff.extend([(XY_pulse, (2, 0)), (Z_pulse, 1)])

Control and noise operator identifiers can be mapped according to a specified mapping:

>>> YX_pulse = ff.extend([(X_pulse, 1, {'X': 'IX', 'Z': 'IZ'}),
...                       (Y_pulse, 0, {'Y': 'YI', 'Z': 'ZI'})])
>>> YX_pulse.c_oper_identifiers
array(['IX', 'YI'], dtype='<U2')
>>> YX_pulse.n_oper_identifiers
array(['IX', 'IZ', 'YI', 'ZI'], dtype='<U2')

We can also add an additional noise Hamiltonian:

>>> H_n = [[ff.util.tensor(Z, Z, Z), [1], 'ZZZ']]
>>> XYX_pulse = ff.extend([(XX_pulse, (0, 2)), (Y_pulse, 1)],
...                       additional_noise_Hamiltonian=H_n)
>>> 'ZZZ' in XYX_pulse.n_oper_identifiers
True

infidelity(pulse: PulseSequence, spectrum: Union[Sequence[float], Callable], omega: Union[Sequence[float], Dict[str, Union[int, str]]], n_oper_identifiers: Optional[Sequence[str]] = None, which: str = 'total', show_progressbar: bool = False, cache_intermediates: bool = False, return_smallness: bool = False, test_convergence: bool = False) → Union[ndarray, Any][source]

Calculate the leading order entanglement infidelity.

This function calculates the infidelity approximately from the leading peturbation (see Notes). To compute it exactly for Gaussian noise and vanishing coherent errors (second order Magnus terms), use error_transfer_matrix() to obtain it from the full process matrix.

Parameters

pulse: PulseSequence

The PulseSequence instance for which to calculate the infidelity for.

spectrum: array_like, shape ([[n_nops,] n_nops,] omega) or callable

The two-sided noise power spectral density in units of inverse frequencies as an array of shape (n_omega,), (n_nops, n_omega), or (n_nops, n_nops, n_omega). In the first case, the same spectrum is taken for all noise operators, in the second, it is assumed that there are no correlations between different noise sources and thus there is one spectrum for each noise operator. In the third and most general case, there may be a spectrum for each pair of noise operators corresponding to the correlations between them. n_nops is the number of noise operators considered and should be equal to len(n_oper_identifiers).

If test_convergence is True, a function handle to compute the power spectral density from a sequence of frequencies is expected.

omega: array_like or dict

The frequencies at which the integration is to be carried out. If test_convergence is True, a dict with possible keys (‘omega_IR’, ‘omega_UV’, ‘spacing’, ‘n_min’, ‘n_max’, ‘n_points’), where all entries are integers except for spacing which should be a string, either ‘linear’ or ‘log’. ‘n_points’ controls how many steps are taken.

n_oper_identifiers: array_like, optional

The identifiers of the noise operators for which to calculate the infidelity contribution. If given, the infidelities for each noise operator will be returned. Otherwise, all noise operators will be taken into account.

which: str, optional

Which infidelities should be calculated, may be either ‘total’ (default) or ‘correlations’. In the former case, one value is returned for each noise operator, corresponding to the total infidelity of the pulse (or pulse sequence). In the latter, an array of infidelities is returned where element (i,j) corresponds to the infidelity contribution of the correlations between pulses i and j (see Notes). Note that this option is only available if the pulse correlation filter functions have been computed during concatenation (see calculate_pulse_correlation_filter_function() and concatenate()).

show_progressbar: bool, optional

Show a progressbar for the calculation of the control matrix.

cache_intermediates: bool, optional

Keep and return intermediate terms of the calculation of the control matrix (if it is not already cached) that can be reused for second order or gradients. The default is False.

return_smallness: bool, optional

Return the smallness parameter \(\xi\) for the given spectrum.

test_convergence: bool, optional

Test the convergence of the integral with respect to the number of frequency samples. Returns the number of frequency samples and the corresponding fidelities. See spectrum and omega for more information.

Returns

infid: ndarray, shape ([[n_pls, n_pls,], n_nops,] n_nops): Array with the infidelity contributions for each spectrum spectrum on the last axis or axes, depending on the shape of spectrum and which. If which is correlations, the first two axes are the individual pulse contributions. If spectrum is 2-d (3-d), the last axis (two axes) are the individual spectral contributions. Only if test_convergence is False.
n_samples: array_like: Array with number of frequency samples used for convergence test. Only if test_convergence is True.
convergence_infids: array_like: Array with infidelities calculated in convergence test. Only if test_convergence is True.

See also

calculate_decay_amplitudes: Calculate the full matrix of first order terms.
error_transfer_matrix: Calculate the full process matrix.
plotting.plot_infidelity_convergence: Convenience function to plot results.
pulse_sequence.concatenate: Concatenate PulseSequence objects.
calculate_pulse_correlation_filter_function.

Notes

The infidelity is given by

\[\begin{split}\mathcal{I}_{\alpha\beta} &= 1 - \frac{1}{d^2}\mathrm{tr}\:\tilde{\mathcal{U}} \\ &= \frac{1}{d}\int_{-\infty}^{\infty} \frac{\mathrm{d}\omega}{2\pi}S_{\alpha\beta}(\omega) F_{\alpha\beta}(\omega) + \mathcal{O}\big(\xi^4\big) \\ &= \sum_{g,g'=1}^G \mathcal{I}_{\alpha\beta}^{(gg')}\end{split}\]

with \(S_{\alpha\beta}(\omega)\) the two-sided noise spectral density and \(F_{\alpha\beta}(\omega)\) the first-order filter function for noise sources \(\alpha,\beta\). The noise spectrum may include correlated noise sources, that is, its entry at \((\alpha,\beta)\) corresponds to the correlations between sources \(\alpha\) and \(\beta\). \(\mathcal{I}_{\alpha\beta}^{(gg')}\) are the correlation infidelities that can be computed by setting which='correlations'.

To convert to the average gate infidelity, use the following relation given by Horodecki et al. [Hor99] and Nielsen [Nie02]:

\[\mathcal{I}_\mathrm{avg} = \frac{d}{d+1}\mathcal{I}.\]

The smallness parameter is given by

\[\xi^2 = \sum_\alpha\left[ \lvert\lvert B_\alpha\rvert\rvert^2 \int_{-\infty}^\infty\frac{\mathrm{d}\omega}{2\pi} S_\alpha(\omega)\left(\sum_gs_\alpha^{(g)}\Delta t_g \right)^2 \right].\]

Note that in practice, the integral is only evaluated on the interval \(\omega\in[\omega_\mathrm{min},\omega_\mathrm{max}]\).

References

Hor99: Horodecki, M., Horodecki, P., & Horodecki, R. (1999). General teleportation channel, singlet fraction, and quasidistillation. Physical Review A - Atomic, Molecular, and Optical Physics, 60(3), 1888–1898. https://doi.org/10.1103/PhysRevA.60.1888
Nie02: Nielsen, M. A. (2002). A simple formula for the average gate fidelity of a quantum dynamical operation. Physics Letters, Section A: General, Atomic and Solid State Physics, 303(4), 249–252. https://doi.org/10.1016/S0375-9601(02)01272-0

infidelity_derivative(pulse: PulseSequence, spectrum: Sequence[float], omega: Sequence[float], control_identifiers: Optional[Sequence[str]] = None, n_oper_identifiers: Optional[Sequence[str]] = None, n_coeffs_deriv: Optional[Sequence[Sequence[float]]] = None) → ndarray[source]

Calculate the entanglement infidelity derivative of the PulseSequence pulse.

Parameters

pulse: PulseSequence: The PulseSequence instance for which to calculate the infidelity.
spectrum: array_like, shape ([[n_nops,] n_nops,] omega): The two-sided noise power spectral density in units of inverse frequencies as an array of shape (n_omega,), (n_nops, n_omega), or (n_nops, n_nops, n_omega). In the first case, the same spectrum is taken for all noise operators, in the second, it is assumed that there are no correlations between different noise sources and thus there is one spectrum for each noise operator. In the third and most general case, there may be a spectrum for each pair of noise operators corresponding to the correlations between them. n_nops is the number of noise operators considered and should be equal to len(n_oper_identifiers).
omega: array_like, shape (n_omega,): The frequencies at which the integration is to be carried out.
control_identifiers: Sequence[str], shape (n_ctrl,): Sequence of strings with the control identifiers to distinguish between control and drift Hamiltonian. The default is None, in which case the derivative is computed for all known non-noise operators.
n_oper_identifiers: Sequence[str], shape (n_nops,): Sequence of strings with the noise identifiers for which to compute the derivative contribution. The default is None, in which case it is computed for all known noise operators.
n_coeffs_deriv: array_like, shape (n_nops, n_ctrl, n_dt): The derivatives of the noise susceptibilities by the control amplitudes. The rows and columns should be in the same order as the corresponding identifiers above. Defaults to None, in which case the coefficients are assumed to be constant and hence their derivative vanishing.

Warning

Internally, control and noise terms of the Hamiltonian are stored alphanumerically sorted by their identifiers. If the noise and/or control identifiers above are not explicitly given, the rows and/or columns of this parameter need to be sorted in the same fashion.

Returns

infid_deriv: ndarray, shape (n_nops, n_dt, n_ctrl): Array with the derivative of the infidelity for each noise source taken for each control direction at each time step \(\frac{\partial I_e}{\partial u_h(t_{g'})}\). Sorted in the same fashion as n_coeffs_deriv or, if not given, alphanumerically by the identifiers.

Raises

ValueError: If the provided noise spectral density does not fit expected shape.

Notes

The infidelity’s derivative is given by

\[\frac{\partial I_e}{\partial u_h(t_{g'})} = \frac{1}{d} \int_{-\infty}^\infty \frac{d\omega}{2\pi} S_\alpha(\omega) \frac{\partial F_\alpha(\omega)} {\partial u_h(t_{g'})}\]

with \(S_{\alpha}(\omega)\) the noise spectral density and \(F_{\alpha}(\omega)\) the canonical filter function for noise source \(\alpha\).

To convert to the average gate infidelity, use the following relation given by Horodecki et al. [Hor99] and Nielsen [Nie02]:

\[\big\langle\mathcal{I}_\mathrm{avg}\big\rangle = \frac{d}{d+1} \big\langle\mathcal{I}_\mathrm{e}\big\rangle.\]

References

Hor99: Horodecki, M., Horodecki, P., & Horodecki, R. (1999). General teleportation channel, singlet fraction, and quasidistillation. Physical Review A - Atomic, Molecular, and Optical Physics, 60(3), 1888–1898. https://doi.org/10.1103/PhysRevA.60.1888
Nie02: Nielsen, M. A. (2002). A simple formula for the average gate fidelity of a quantum dynamical operation. Physics Letters, Section A: General, Atomic and Solid State Physics, 303(4), 249–252. https://doi.org/10.1016/S0375-9601(02)01272-0

liouville_representation(U: ndarray, basis: Basis) → ndarray[source]

Get the Liouville representaion of the unitary U with respect to the basis.

Parameters

U: ndarray, shape (…, d, d): The unitary.
basis: Basis, shape (d**2, d, d): The basis used for the representation, e.g. a Pauli basis.

Returns

R: ndarray, shape (…, d**2, d**2): The Liouville representation of U.

Notes

The Liouville representation of a unitary quantum operation \(\mathcal{U}:\rho\rightarrow U\rho U^\dagger\) is given by

\[\mathcal{U}_{ij} = \mathrm{tr}(C_i U C_j U^\dagger)\]

with \(C_i\) elements of the basis spanning \(\mathbb{C}^{d\times d}\) with \(d\) the dimension of the Hilbert space.

remap(pulse: PulseSequence, order: Sequence[int], d_per_qubit: int = 2, oper_identifier_mapping: Optional[Mapping[str, str]] = None) → PulseSequence[source]

Remap a PulseSequence by changing the order of qubits in the register. Cached attributes are automatically attempted to be retained.

Caution

This function simply permutes the order of the tensor product elements of control and noise operators. Thus, the resultant pulse will have its filter functions defined for different noise operators than the original one.

Parameters

pulse: PulseSequence: The pulse whose qubit order should be permuted.
order: sequence of ints: A list of permutation indices. E.g., if pulse is defined for two qubits, order == [1, 0] will reverse the order of qubits.
d_per_qubit: int (default: 2): The size of the Hilbert space a single qubit inhabitates.
oper_identifier_mapping: dict_like: A mapping that maps operator identifiers from the old pulse to the remapped pulse. The default is the identity mapping.

Returns

remapped_pulse: PulseSequence: A new PulseSequence instance with the order of the qubits permuted according to order.

See also

extend: Map PulseSequences to composite Hilbert spaces.
util.tensor_transpose: Transpose the order of a tensor product.

Examples

>>> X, Y = util.paulis[1:3]
>>> XY, YX = util.tensor(X, Y), util.tensor(Y, X)
>>> pulse = PulseSequence([[XY, [np.pi/2], 'XY']], [[YX, [1], 'YX']], [1],
...                       Basis.pauli(2))
>>> mapping = {'XY': 'YX', 'YX': 'XY'}
>>> remapped_pulse = remap(pulse, (1, 0), oper_identifier_mapping=mapping)
>>> target_pulse = PulseSequence([[YX, [np.pi/2], 'YX']],
...                              [[XY, [1], 'XY']], [1], Basis.pauli(2))
>>> remapped_pulse == target_pulse
True

Caching of attributes is automatically handled >>> remapped_pulse.is_cached(‘filter_function’) False >>> pulse.cache_filter_function(util.get_sample_frequencies(pulse)) >>> remapped_pulse = remap(pulse, (1, 0)) >>> remapped_pulse.is_cached(‘filter_function’) True