Source code for filter_functions.gradient

# -*- coding: utf-8 -*-
# =============================================================================
#     filter_functions
#     Copyright (C) 2019 Quantum Technology Group, RWTH Aachen University
#
#     This program is free software: you can redistribute it and/or modify
#     it under the terms of the GNU General Public License as published by
#     the Free Software Foundation, either version 3 of the License, or
#     (at your option) any later version.
#
#     This program is distributed in the hope that it will be useful,
#     but WITHOUT ANY WARRANTY; without even the implied warranty of
#     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
#     GNU General Public License for more details.
#
#     You should have received a copy of the GNU General Public License
#     along with this program. If not, see <http://www.gnu.org/licenses/>.
#
#     This module is an extension of the filter_functions package written by
#     Tobias Hangleiter. Its implementation was center of the Bachelor thesis
#     'Filter Function Derivative for Quantum Optimal Control' by Isabel Le.
#
#     Contact email: isabel.le@rwth-aachen.de
#
#     The code has been extended by Julian Teske such that the noise
#     susceptibilities or noise coefficients can depend on the control
#     amplitudes as well.
#
#     Contact email: j.teske@fz-juelich.de
# =============================================================================
"""
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.

Functions
---------
:func:`calculate_derivative_of_control_matrix_from_scratch`
    Calculate the derivative of the control matrix from scratch.
:func:`calculate_canonical_filter_function_derivative`
    Compute the filter function derivative from the control matrix.
:func:`infidelity_derivative`
    Calculate the infidelity derivative.
"""
from typing import Dict, Optional, Sequence, Tuple

import numpy as np
import opt_einsum as oe
from opt_einsum.contract import ContractExpression
from numpy import ndarray

from . import numeric, superoperator, util
from .basis import Basis
from .types import Coefficients, Operator

__all__ = ['calculate_derivative_of_control_matrix_from_scratch',
           'calculate_filter_function_derivative', 'infidelity_derivative']


def _derivative_integral(E: Coefficients, eigvals: Coefficients, dt: float,
                         out: ndarray) -> ndarray:
    """
    Compute the integral appearing in the derivative of the control
    matrix. Result (out) has shape (len(E), d, d, d, d).
    """
    # Precompute masks and energy differences
    dE = np.subtract.outer(eigvals, eigvals)
    mask_dE = np.abs(dE) < 1e-7
    EdE = np.add.outer(E, dE)
    mask_EdE = np.abs(EdE) < 1e-7
    EdEdE = np.add.outer(EdE, dE[~mask_dE])
    mask_EdEdE = np.abs(EdEdE) < 1e-7

    # Case Omega_pq == 0
    tmp1 = np.divide(util.cexp(EdE*dt), EdE, where=~mask_EdE)
    tmp2 = tmp1 - np.divide(1, EdE, where=~mask_EdE)
    tmp2[mask_EdE] = 1j * dt

    tmp1 *= -1j * dt
    tmp1 += np.divide(tmp2, EdE, where=~mask_EdE)
    tmp1[mask_EdE] = dt**2 / 2

    out[:, mask_dE] = tmp1[:, None]

    # Case Omega_pq != 0
    tmp1 = np.divide(1 - util.cexp(EdEdE*dt, where=~mask_EdEdE), EdEdE, where=~mask_EdEdE)
    tmp1[mask_EdEdE] = -1j * dt
    tmp1 += tmp2[..., None]

    out[:, ~mask_dE] = (tmp1 / dE[~mask_dE]).transpose(0, 3, 1, 2)
    return out


def _liouville_derivative(dt: Coefficients, propagators: ndarray, basis: Basis, eigvecs: ndarray,
                          eigvals: ndarray, c_opers_transformed: ndarray) -> ndarray:
    r"""
    Calculate the derivatives of the control propagators in Liouville
    representation.

    Parameters
    ----------
    dt: array_like, shape (n_dt)
        Sequence duration, i.e. for the :math:`g`-th pulse
        :math:`t_g - t_{g-1}`.
    propagators: array_like, shape (n_dt+1, d, d)
        The propagators :math:`Q_g = P_g P_{g-1}\cdots P_0` as a (d, d)
        array with *d* the dimension of the Hilbert space.
    basis: Basis, shape (d**2, d, d)
        The basis elements, in which the pulse control matrix will be
        expanded.
    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, ...])``.
    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, ...])``.
    c_opers_transformed: array_like, shape (n_dt, c_ctrl, d, d)
        The control operators transformed into the eigenspace of the
        control Hamiltonian. The drift operators are ignored, if
        identifiers for accessible control operators are provided.

    Returns
    -------
    liouville_deriv: array_like, shape (n_dt, n_ctrl, n_dt, d**2, d**2)
        The derivative of the control propagators in Liouville
        representation :math:`\frac{\partial \mathcal{Q}_{jk}^{(g)}}
        {\partial u_h(t_{g^\prime})}`.

    Notes
    -----
    The derivatives of the control propagators in Liouville
    representation are calculated according to

    .. math::

        \frac{\partial\mathcal{Q}_{jk}^{(g-1)}}{\partial u_h(t_{g^\prime})} &=
        \Theta_{g-1}(g^\prime) \mathrm{tr}\Big(
            \frac{\partial U_c^\dagger(t_{g-1},0)}{\partial u_h(t_{g^\prime})}
            C_j U_c(t_{g-1},0) C_k\Big)\\
            &+\Theta_{g-1}(g^\prime)\mathrm{tr}\Big(U_c^\dagger(t_{g-1},0)C_j
                                                \frac{\partial U_c(t_{g-1},0)}
                                                {\partial u_h(t_{g^\prime})}C_k
                                                \Big),

    where :math:`\Theta_{g-1}(g^\prime)` being 1 if :math:`g^\prime < g-1` and
    zero otherwise.

    """
    n, d = eigvecs.shape[:2]
    # omega_i - omega_j
    omega_diff = eigvals[:, :, None] - eigvals[:, None, :]
    dt_broadcast = np.broadcast_to(dt[:, None, None], omega_diff.shape)
    # mask = omega_diff == 0
    mask = np.broadcast_to(np.eye(d, dtype=bool), omega_diff.shape)
    A_mat = np.empty(omega_diff.shape, dtype=complex)
    A_mat[mask] = dt_broadcast[mask]
    A_mat[~mask] = 1j*(1 - util.cexp(omega_diff[~mask]*dt_broadcast[~mask])) / omega_diff[~mask]
    U_deriv = -1j * (propagators[1:]
                     @ propagators[:-1].conj().swapaxes(-1, -2)
                     @ eigvecs
                     @ (A_mat * c_opers_transformed.swapaxes(0, 1))
                     @ eigvecs.conj().swapaxes(-1, -2))

    # Calculate the whole propagator derivative
    propagators_deriv = np.zeros((c_opers_transformed.shape[1], n-1, n, d, d), dtype=complex)
    U_deriv_transformed = np.zeros((c_opers_transformed.shape[1], n-1, d, d), dtype=complex)
    for g in range(n-1):
        U_deriv_transformed[:, g] = (propagators[g+1].conj().swapaxes(-1, -2)
                                     @ U_deriv[:, g]
                                     @ propagators[g])
        propagators_deriv[:, g, :g+1] = propagators[g+1] @ U_deriv_transformed[:, :g+1]

    # Equivalent but usually slower contraction: 'htsba,jbc,tcd,kda->thsjk'
    # Can just take 2*Re(·) when calculating x + x*
    liouville_deriv = np.einsum('htsba,tjkba->thsjk', propagators_deriv.conj(),
                                (basis @ propagators[1:-1, None])[:, :, None] @ basis).real
    liouville_deriv *= 2
    return liouville_deriv


def _control_matrix_at_timestep_derivative(
        omega: Coefficients,
        dt: Coefficients,
        eigvals: ndarray,
        eigvecs: ndarray,
        basis_transformed,
        c_opers_transformed: ndarray,
        n_opers_transformed: ndarray,
        n_coeffs: Sequence[Coefficients],
        n_coeffs_deriv: Sequence[Coefficients],
        phase_factor: ndarray,
        integral: ndarray,
        deriv_integral: ndarray,
        ctrlmat_step: ndarray,
        ctrlmat_expr: ContractExpression,
) -> Tuple[ndarray, ndarray]:
    r"""Calculate the control matrices and corresponding derivatives.

    Calculate control matrices at each time step and the corresponding
    partial derivatives of those with respect to control strength at
    each time step.

    Parameters
    ----------
    omega: array_like, shape (n_omega)
        Frequencies, at which the pulse control matrix is to be
        evaluated.
    dt: array_like, shape (n_dt)
        Sequence duration, i.e. for the :math:`g`-th pulse
        :math:`t_g - t_{g-1}`.
    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.
        ``D == 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.
        ``V == array([V_0, V_1, ...])``.
    basis_transformed: array_like, shape (d**2, d, d)
        The basis elements in which the pulse control matrix will be
        expanded transformed by the eigenvectors.
    c_opers_transformed: array_like, shape (n_ctrl, d, d)
        The control operators transformed into the eigenspace of the
        control Hamiltonian.
    n_opers_transformed: array_like, shape (n_nops, d, d)
        The noise operators transformed into the eigenspace of the
        control Hamiltonian.
    n_coeffs: array_like, shape (n_nops,)
        The sensitivities of the system to the noise operators given by
        *n_opers_transformed* at the given time step.
    n_coeffs_deriv: array_like, shape (n_nops, n_ctrl,)
        The derivatives of the noise susceptibilities by the control
        amplitudes. Defaults to None.
    phase_factor: array_like, shape (n_omega,)
        The phase factor :math:`e^{i\omega t_{g-1}}`.
    integral: array_like, shape (n_omega, d, d)
        The integral during the time step appearing in the regular
        control matrix.
    deriv_integral: array_like, shape (n_omega, d, d, d, d)
        An array to write the integral during the time step appearing in
        the derivative into.
    ctrlmat_step: array_like, shape (n_nops, d**2, n_omega)
        An array to write the control matrix during the time step into.
    ctrlmat_expr: ContractExpression
        An :class:`opt_einsum.contract.ContractExpression` to compute
        the control matrix during the time step.

    Returns
    -------
    ctrlmat_g: ndarray, shape (n_dt, n_nops, d**2, n_omega)
        The individual control matrices of all time steps
    ctrlmat_g_deriv: ndarray, shape (n_dt, n_nops, d**2, n_ctrl, n_omega)
        The corresponding derivative with respect to the control
        strength :math:`\frac{\partial\mathcal{B}_{\alpha j}^{(g)}(\omega)}

    Notes
    -----
    The control matrix at each time step is evaluated according to

    .. math::

            \mathcal{B}_{\alpha j}^{(g)}(\omega) = s_\alpha^{(g)}\mathrm{tr}
            \left([\bar{B}_\alpha \circ I_1^{(g)}(\omega)] \bar{C}_j \right),

    where

    .. math::

        I_{1,nm}^{(g)}(\omega) = \frac{\exp(\mathrm{i}(\omega + \omega_n^{(g)}
                                            - \omega_m^{(g)}) \Delta t_g) - 1}
        {\mathrm{i}(\omega + \omega_n^{(g)} - \omega_m^{(g)})}

    The derivative of the control matrix with respect to the control
    strength at different time steps is calculated according to

    .. math::

        \frac{\partial \mathcal{B}_{\alpha j}^{(g)}(\omega)}
        {\partial u_h(t_{g^\prime})} =
        \mathrm{i}\delta_{gg^\prime} s_\alpha^{(g)} \mathrm{tr}
        \left( \bar{B}_{\alpha} \cdot \mathbb{M} \right)
        + \frac{\partial s_\alpha^{(g)}}{u_h (t_{g^\prime})} \text{tr}
        \left( (\overline{B}_{\alpha} \circ I_1^{(g)}{}(\omega)) \cdot
        \overline{C}_{j}) \right).

    We assume that the noise susceptibility :math:`s` only depends
    locally on the time i.e. :math:`\partial_{u(t_g)} s(t_{g^\prime})
    = \delta_{gg^\prime} \partial_{u(t_g)} s(t_g)`
    If denoting :math:`\Delta\omega_{ij} = \omega_i^{(g)} -
    \omega_j^{(g)}` the integral part is encapsulated in

    .. math::

        \mathbb{M}_{mn} = \left[ \bar{C}_j, \mathbb{I}^{(mn)}
                                \circ \bar{H}_h \right]_{mn},

    with

    .. math::

        \mathbb{I}_{pq}^{(mn)} &= \delta_{pq} \left(
            \frac{\Delta t_g \cdot
                  \exp(\mathrm{i}(\omega + \Delta\omega_{nm})\Delta t_g)}
            {\mathrm{i}(\omega + \Delta\omega_{nm})}
            + \frac{\exp(\mathrm{i}(\omega + \Delta\omega_{nm})\Delta t_g) - 1}
            {(\omega + \Delta\omega_{nm})^2}\right)\\
            &+  \frac{1-\delta_{pq}}{\mathrm{i}\Delta\omega_{pq}} \cdot
            \frac{\exp(\mathrm{i}(\omega + \Delta\omega_{nm}
                                  + \Delta\omega_{pq})\Delta t_g) - 1}
            {\mathrm{i}(\omega + \Delta\omega_{nm} + \Delta\omega_{pq})}\\
            &- \frac{1-\delta_{pq}}{\mathrm{i}\Delta\omega_{pq}} \cdot
            \frac{\exp(\mathrm{i}(\omega + \Delta\omega_{nm})\Delta t_g) - 1}
            {\mathrm{i}(\omega + \Delta\omega_{nm})}
    """
    d = len(eigvecs)
    d2 = d**2
    n_ctrl = len(c_opers_transformed)
    n_nops = len(n_opers_transformed)
    n_omega = len(omega)

    deriv_integral = _derivative_integral(omega, eigvals, dt, out=deriv_integral)
    ctrlmat_step = ctrlmat_expr(phase_factor, basis_transformed, n_opers_transformed, integral,
                                out=ctrlmat_step)

    # Basically a tensor product
    # K = np.einsum('taij,thkl->tahikjl', VBV, VHV)
    # L = K.transpose(0, 1, 2, 4, 3, 6, 5)
    K = util.tensor(n_opers_transformed[:, None], c_opers_transformed[None])
    L = util.tensor_transpose(K, (1, 0), [[d, d], [d, d]])
    k = np.diagonal(K.reshape(n_nops, n_ctrl, d, d, d, d), 0, -2, -3)
    l = np.diagonal(L.reshape(n_nops, n_ctrl, d, d, d, d), 0, -2, -3)
    i1 = np.diagonal(deriv_integral, 0, -2, -3)
    i2 = np.diagonal(deriv_integral, 0, -1, -4)
    # reshaping in F-major is significantly faster than C-major (~ factor 2-4)
    M = np.einsum(
        'ahpm,opm->ahop',
        l.reshape(n_nops, n_ctrl, d2, d, order='F'),
        i1.reshape(n_omega, d2, d, order='F')
    ).reshape(n_nops, n_ctrl, n_omega, d, d, order='F')
    if d == 2:
        # Life a bit simpler
        mask = np.eye(d, dtype=bool)
        M[..., mask] -= M[..., mask][..., ::-1]
        M[..., ~mask] *= 2
    else:
        M -= np.einsum(
            'ahpn,opn->ahop',
            k.swapaxes(-2, -3).reshape(n_nops, n_ctrl, d2, d, order='F'),
            i2.reshape(n_omega, d2, d, order='F')
        ).reshape(n_nops, n_ctrl, n_omega, d, d, order='F').swapaxes(-1, -2)

    # Expand in basis transformed to eigenspace. Include phase factor and
    # factor 1j here to make use of optimized contraction order
    ctrlmat_step_deriv = oe.contract('o,jnk,ahokn->ajho', phase_factor, 1j*basis_transformed, M,
                                     optimize=[(1, 2), (0, 1)])

    if n_coeffs_deriv is not None:
        # equivalent contraction: 'ah,a,ako->akho', but this faster
        ctrlmat_step_deriv += ((n_coeffs_deriv / n_coeffs[:, None])[:, None, :, None]
                               * ctrlmat_step[:, :, None])

    return ctrlmat_step, ctrlmat_step_deriv


[docs]def calculate_derivative_of_control_matrix_from_scratch( omega: Coefficients, propagators: ndarray, eigvals: ndarray, eigvecs: ndarray, basis: Basis, t: Coefficients, dt: Coefficients, n_opers: Sequence[Operator], n_coeffs: Sequence[Coefficients], c_opers: Sequence[Operator], all_identifiers: Sequence[str], control_identifiers: Optional[Sequence[str]] = None, n_coeffs_deriv: Optional[Sequence[Coefficients]] = None, intermediates: Optional[Dict[str, ndarray]] = None ) -> ndarray: r""" 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 :math:`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 :math:`g`-th pulse :math:`t_g - t_{g-1}`. n_opers: array_like, shape (n_nops, d, d) Noise operators :math:`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_cops, d, d) Control operators :math:`H_k`. all_identifiers: array_like, shape (n_cops) Identifiers of all control operators. control_identifiers: Sequence[str], shape (n_ctrl), Optional Sequence of strings with the control identifiers to distinguish between accessible control and drift Hamiltonian. The default is None. 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. Raises ------ ValueError If the given identifiers *control_identifier* are not subset of the total identifiers *all_identifiers* of all control operators. 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 :math:`\frac{\partial\mathcal{B}_{\alpha k}(\omega)}{\partial u_h(t_{g'})}`. Notes ----- The derivative of the control matrix is calculated according to .. math :: \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) See Also -------- _liouville_derivative _control_matrix_at_timestep_derivative """ # Distinction between control and drift operators and only # calculate the derivatives in control direction try: idx = util.get_indices_from_identifiers(all_identifiers, control_identifiers) except ValueError as err: raise ValueError('Given control identifiers have to be a subset of (drift+control) ' + 'Hamiltonian!') from err d = eigvecs.shape[-1] n_dt = len(dt) n_ctrl = len(idx) n_nops = len(n_opers) n_omega = len(omega) # Precompute some transformations or grab from cache if possible basis_transformed = numeric._transform_by_unitary(eigvecs[:, None], basis[None], out=np.empty((n_dt, d**2, d, d), complex)) c_opers_transformed = numeric._transform_hamiltonian(eigvecs, c_opers[idx]).swapaxes(0, 1) if not intermediates: # None or empty n_opers_transformed = numeric._transform_hamiltonian(eigvecs, n_opers, n_coeffs).swapaxes(0, 1) exp_buf, integral = np.empty((2, n_omega, d, d), dtype=complex) else: n_opers_transformed = intermediates['n_opers_transformed'].swapaxes(0, 1) propagators_liouville = superoperator.liouville_representation(propagators[:-1], basis) propagators_liouville_deriv = _liouville_derivative(dt, propagators, basis, eigvecs, eigvals, c_opers_transformed) deriv_integral = np.empty((n_omega, d, d, d, d), dtype=complex) ctrlmat_deriv = np.empty((n_ctrl, n_omega, n_dt, n_nops, d**2), dtype=complex) ctrlmat_step = np.empty((n_dt, n_nops, d**2, n_omega), dtype=complex) # Optimized expression that is passed to control_matrix_at_timestep_derivative # in each iteration ctrlmat_expr = oe.contract_expression('o,icd,adc,odc->aio', (len(omega),), basis.shape, n_opers.shape, (len(omega), d, d), optimize=[(0, 3), (0, 1), (0, 1)]) for g in range(n_dt): if not intermediates: integral = numeric._first_order_integral(omega, eigvals[g], dt[g], exp_buf, integral) else: integral = intermediates['first_order_integral'][g] n_coeff_deriv = n_coeffs_deriv if n_coeffs_deriv is None else n_coeffs_deriv[:, :, g] # ctrlmat_step is computed from scratch because the quantity in # the cache (from numeric.calculate_control_matrix_from_scratch) # contains the Liouville propagators already and it is expensive # to remove them by multiplying with the transpose. ctrlmat_step[g], ctrlmat_step_deriv = _control_matrix_at_timestep_derivative( omega, dt[g], eigvals[g], eigvecs[g], basis_transformed[g], c_opers_transformed[g], n_opers_transformed[g], n_coeffs[:, g], n_coeff_deriv, util.cexp(omega*t[g]), integral, deriv_integral, ctrlmat_step[g], ctrlmat_expr ) # Phase factor already part of ctrlmat_step_deriv ctrlmat_deriv[:, :, g] = (ctrlmat_step_deriv.transpose(2, 3, 0, 1)
[docs] @ propagators_liouville[g]) # opt_einsum a lot faster here # Phase factor again already part of ctrlmat_step ctrlmat_deriv += oe.contract('tajo,thsjk->hosak', ctrlmat_step[1:], propagators_liouville_deriv) return ctrlmat_deriv
def calculate_filter_function_derivative(ctrlmat: ndarray, ctrlmat_deriv: ndarray) -> ndarray: r""" 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 :math:`\frac{\partial F_\alpha(\omega)}{\partial u_h(t_{g'})}`. Notes ----- The filter function derivative is calculated according to .. math :: \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) """ return 2*np.einsum('ako,hotak->atho', ctrlmat.conj(), ctrlmat_deriv).real
[docs]def infidelity_derivative( pulse: 'PulseSequence', spectrum: Coefficients, omega: Coefficients, control_identifiers: Optional[Sequence[str]] = None, n_coeffs_deriv: Optional[Sequence[Coefficients]] = None ) -> ndarray: r"""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 identifiern to distinguish between accessible control and drift Hamiltonian. 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. Raises ------ ValueError If the provided noise spectral density does not fit expected shape. 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 :math:`\frac{\partial I_e}{\partial u_h(t_{g'})}`. Notes ----- The infidelity's derivative is given by .. math:: \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 :math:`S_{\alpha}(\omega)` the noise spectral density and :math:`F_{\alpha}(\omega)` the canonical filter function for noise source :math:`\alpha`. To convert to the average gate infidelity, use the following relation given by Horodecki et al. [Hor99]_ and Nielsen [Nie02]_: .. math:: \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 """ spectrum = numeric._parse_spectrum(spectrum, omega, range(len(pulse.n_opers))) filter_function_deriv = pulse.get_filter_function_derivative(omega, control_identifiers, n_coeffs_deriv) integrand = np.einsum('ao,atho->atho', spectrum, filter_function_deriv) infid_deriv = util.integrate(integrand, omega) / (2*np.pi*pulse.d) return infid_deriv