# Source code for weylchamber.local_invariants

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.qip.operations.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.qip.operations.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()