import numpy as np
from numpy import sin, cos, pi
import scipy
from scipy.optimize import leastsq
import logging
import qutip
from .prec import DEFAULT_WEYL_PRECISSION
from ._types import Gate, GTuple
from .coordinates import to_magic, c1c2c3, _SQ_unitary
from .cartan_decomposition import canonical_gate
__all__ = [
'g1g2g3',
'g1g2g3_from_c1c2c3',
'J_T_LI',
'closest_LI',
'make_LI_krotov_chi_constructor',
]
[docs]def g1g2g3(U: Gate, ndigits=DEFAULT_WEYL_PRECISSION) -> GTuple:
"""Calculate local invariants $(g_1, g_3, g_3)$
Given a two-qubit gate, calculate local invariants $(g_1, g_2, g_3)$.
U must be in the canonical basis. For numerical stability, the resulting
values are rounded to the given precision, cf. the `ndigits` parameter of
the built-in :func:`round` function.
>>> print("%.2f %.2f %.2f" % g1g2g3(qutip.gates.cnot()))
0.00 0.00 1.00
"""
# mathematically, the determinant of U and UB is the same, but
# we seem to get better numerical accuracy if we calculate detU with
# the rotated U
UB = to_magic(U).full() # instance of np.ndarray
detU = np.linalg.det(UB)
m = UB.T @ UB
g1_2 = (np.trace(m)) ** 2 / (16.0 * detU)
g3 = (np.trace(m) ** 2 - np.trace(m @ m)) / (4.0 * detU)
g1 = round(g1_2.real + 0.0, ndigits) # adding 0.0 turns -0.0 into +0.0
g2 = round(g1_2.imag + 0.0, ndigits)
g3 = round(g3.real + 0.0, ndigits)
return (g1, g2, g3)
[docs]def g1g2g3_from_c1c2c3(
c1: float, c2: float, c3: float, ndigits=DEFAULT_WEYL_PRECISSION
) -> GTuple:
"""Calculate local invariants from the Weyl chamber coordinates
Calculate the local invariants $(g_1, g_2, g_3)$ from the Weyl chamber
coordinates $(c_1, c_2, c_3)$, in units of π. The result is rounded to the
given precision, in order to enhance numerical stability (cf. `ndigits`
parameter of the built-in :func:`round` function)
Example:
>>> CNOT = qutip.gates.cnot()
>>> print("%.2f %.2f %.2f" % g1g2g3_from_c1c2c3(*c1c2c3(CNOT)))
0.00 0.00 1.00
"""
c1 *= pi
c2 *= pi
c3 *= pi
g1 = round(
cos(c1) ** 2 * cos(c2) ** 2 * cos(c3) ** 2
- sin(c1) ** 2 * sin(c2) ** 2 * sin(c3) ** 2
+ 0.0,
ndigits,
)
g2 = round(0.25 * sin(2 * c1) * sin(2 * c2) * sin(2 * c3) + 0.0, ndigits)
g3 = round(4 * g1 - cos(2 * c1) * cos(2 * c2) * cos(2 * c3) + 0.0, ndigits)
return g1, g2, g3
[docs]def J_T_LI(O: Gate, U: Gate, form='g'):
"""Calculate value of the local-invariants functional
Args:
O: The optimal gate
U: The achieved gate
form (str): form of the functional to use, 'g' or 'c'
"""
if form == 'g':
return np.sum(np.abs(np.array(g1g2g3(O)) - np.array(g1g2g3(U))) ** 2)
elif form == 'c':
delta_c = np.array(c1c2c3(O)) - np.array(c1c2c3(U))
return np.prod(cos(np.pi * (delta_c) / 2.0))
else:
raise ValueError("Illegal value for 'form'")
[docs]def closest_LI(
U: Gate, c1: float, c2: float, c3: float, method='leastsq', limit=1.0e-6
):
"""Find the closest gate that has the given Weyl chamber coordinates
The `c1`, `c2`, `c3` are given in units of π
"""
A = canonical_gate(c1, c2, c3)
def f_U(p):
return _SQ_unitary(*p[:8]) * A * _SQ_unitary(*p[8:])
return _closest_gate(U, f_U, n=16, method=method, limit=limit)
def _closest_gate(U, f_U, n, x_max=2 * pi, method='leastsq', limit=1.0e-6):
"""Find the closest gate to U that fulfills the parametrization implied by
the function f_U
Args:
U (Gate): Target gate
f_U (callable): function that takes an array of n values and returns an
gate.
n (integer): Number of parameters (size of the argument of f_U)
x_max (float): Maximum value for each element of the array passed as an
argument to f_U. There is no way to have different a different
range for the different elements
method (str): Name of mimimization method, either 'leastsq' or any of
the gradient-free methods implemented by scipy.optimize.mimize
limit (float): absolute error of the distance between the target gate
and the optimized gate for convergence. The limit is automatically
increased by an order of magnitude every 100 iterations
"""
logger = logging.getLogger(__name__)
logger.debug("_closests_gate with method %s", method)
from scipy.optimize import minimize
if method == 'leastsq':
def f_minimize(p):
d = _vectorize(f_U(p) - U)
return np.concatenate([d.real, d.imag])
else:
def f_minimize(p):
return _norm(U - f_U(p))
dist_min = None
iter = 0
while True:
iter += 1
if iter > 100:
iter = 0
limit *= 10
logger.debug("_closests_gate limit -> %.2e", limit)
p0 = x_max * np.random.random(n)
success = False
if method == 'leastsq':
p, info = leastsq(f_minimize, p0)
U_min = f_U(p)
if info in [1, 2, 3, 4]:
success = True
else:
res = minimize(f_minimize, p0, method=method)
U_min = f_U(res.x)
success = res.success
if success:
dist = _norm(U_min - U)
logger.debug("_closests_gate dist = %.5e", dist)
if dist_min is None:
dist_min = dist
logger.debug("_closests_gate dist_min -> %.5e", dist_min)
else:
logger.debug(
"_closests_gate delta_dist = %.5e", abs(dist - dist_min)
)
if abs(dist - dist_min) < limit:
return U_min
else:
if dist < dist_min:
dist_min = dist
logger.debug(
"_closests_gate dist_min -> %.5e", dist_min
)
def _vectorize(a, order='F'):
"""Return vectorization of multi-dimensional numpy array or matrix `a`
Examples:
>>> a = np.array([1,2,3,4])
>>> _vectorize(a)
array([1, 2, 3, 4])
>>> a = np.array([[1,2],[3,4]])
>>> _vectorize(a)
array([1, 3, 2, 4])
>>> _vectorize(a, order='C')
array([1, 2, 3, 4])
"""
if isinstance(a, qutip.Qobj):
a = a.full()
N = a.size
return np.squeeze(np.asarray(a).reshape((1, N), order=order))
def _norm(v):
"""Calculate the norm of a vector or matrix `v`, matching the inner product
defined in the `inner` routine. An algorithm like
Gram-Schmidt-Orthonormalization will only work if the choice of norm and
inner product are compatible.
If `v` is a vector, the norm is the 2-norm (i.e. the standard Euclidian
vector norm).
If `v` is a matrix, the norm is the Hilbert-Schmidt (aka Frobenius) norm.
Note that the HS norm of a matrix is identical to the 2-norm of any
vectorization of that matrix (e.g. writing the columns of the matrix
underneat each other). Also, the HS norm of the m x 1 matrix is the same as
the 2-norm of the equivalent m-dimensional vector.
"""
if isinstance(v, qutip.Qobj):
v = v.data
if isinstance(v, scipy.sparse.spmatrix):
return scipy.sparse.linalg.norm(v)
else:
return scipy.linalg.norm(v)
[docs]def make_LI_krotov_chi_constructor(gate, canonical_basis, unitarity_weight=0):
r"""Return a constructor for the χ's in an LI optimization.
Return a `chi_constructor` that determines the boundary condition of the
backwards propagation in an optimization towards the local equivalence
class of `gate` in Krotov's method, based on the foward-propagtion of the
Bell states.
Args:
gate (qutip.Qobj): A 4×4 quantum gate, in the `canonical_basis`.
canonical_basis (list[qutip.Qobj]): A list of four basis states that
define the canonical basis $\ket{00}$, $\ket{01}$, $\ket{10}$, and
$\ket{11}$ of the logical subspace.
unitarity_weight (float): A weight in [0, 1] that determines how much
emphasis is placed on maintaining population in the logical
subspace.
Returns:
callable: a function ``chi_constructor(fw_states_T, *args)`` that
receive the result of a foward propagation of the Bell states (obtained
from `canonical_basis` via :func:`weylchamber.gates.bell_basis`), and
returns a list of statex $\ket{\chi}$ that are the boundary condition
for the backward propagation in Krotov's method. Positional arguments
beyond `fw_states_T` are ignored.
"""
# see make_PE_krotov_chi_constructor
raise NotImplementedError()
def _get_a_kl_PE(UB):
"""Return the 4×4 `A_kl` coefficient matrix (:class:`qutip.Qobj`)
for the perfect-entanglers functional, for a given gate `UB` in the Bell
basis.
"""
raise NotImplementedError()