Source code for weylchamber.coordinates

"""Calculation of Weyl chamber coordinates"""

import numpy as np
from numpy import less, greater
from numpy import logical_and as And
from numpy import logical_or as Or
from numpy import less_equal as less_eq
from numpy import greater_equal as greater_eq
import qutip

from .gates import Qmagic, SxSx, SySy, SzSz
from .prec import DEFAULT_WEYL_PRECISSION
from ._types import Gate, CTuple

__all__ = [
'c1c2c3', 'canonical_gate', 'from_magic', 'point_in_region',
'point_in_weyl_chamber', 'random_gate', 'random_weyl_point', 'to_magic',
'weyl_region']

[docs]def c1c2c3(U: Gate, ndigits=DEFAULT_WEYL_PRECISSION) -> CTuple: """Calculate Weyl chamber coordinates $(c_1, c_2, c_3)$ Given U (in canonical basis), calculate the Weyl Chamber coordinates $(c_1, c_2, c_3)$, in units of π. In order to facility numerical stability, the resulting coordinates are rounded to the given precision (cf. ndigits parameter of the built-in round function). Otherwise, rounding errors would likely to result in points that are not in the Weyl chamber, e.g. (0.1, 0.0, 1.0e-13) Algorithm from Childs et al., PRA 68, 052311 (2003). Example: >>> print("%.2f %.2f %.2f" % c1c2c3(qutip.gates.cnot())) 0.50 0.00 0.00 """ U = qutip.Qobj(U, dims=[[2, 2], [2, 2]]) if U.shape != (4, 4): raise ValueError("Gates must have a 4×4 shape") U_tilde = SySy * U.trans() * SySy ev = np.linalg.eigvals( ((U * U_tilde).full())/np.sqrt(complex(np.linalg.det(U.full())))) two_S = np.angle(ev) / np.pi for i in range(len(two_S)): if two_S[i] <= -0.5: two_S[i] += 2.0 S = np.sort(two_S / 2.0)[::-1] # sort decreasing n = int(round(sum(S))) S -= np.r_[np.ones(n), np.zeros(4-n)] S = np.roll(S, -n) M = np.array([[1, 1, 0], [1, 0, 1], [0, 1, 1]]) c1, c2, c3 = np.dot(M, S[:3]) if c3 < 0: c1 = 1 - c1 c3 = -c3 return (round(c1+0.0, ndigits), # adding 0.0 turns -0.0 into +0.0 round(c2+0.0, ndigits), round(c3+0.0, ndigits))
[docs]def point_in_weyl_chamber( c1: float, c2: float, c3: float, raise_exception=False) -> bool: """Check if the coordinates $(c_1, c_2, c_3)$ are inside the Weyl chamber Examples: >>> BGATE = qutip.gates.berkeley() >>> point_in_weyl_chamber(*c1c2c3(BGATE)) True >>> point_in_weyl_chamber(*c1c2c3(qutip.gates.identity([2, 2]))) True The coordinates may also be array-like, in which case a boolean numpy array is returned. >>> res = point_in_weyl_chamber( ... [0.0,0.5,0.8], [1.0,0.25,0.0], [1.0,0.25,0.0]) >>> assert np.all(res == np.array([False, True, True])) If raise_exception is True, raise an :exc:ValueError if any values are outside the Weyl chamber. >>> try: ... point_in_weyl_chamber(1.0, 0.5, 0, raise_exception=True) ... except ValueError as e: ... print(e) (1, 0.5, 0) is not in the Weyl chamber """ result = Or( And(And( less(c1, 0.5), less_eq(c2, c1)), less_eq(c3, c2)), And(And( greater_eq(c1, 0.5), less_eq(c2, 1.0-np.array(c1))), less_eq(c3, c2)) ) if raise_exception: if not np.all(result): if np.isscalar(c1): # assume c2, c3 are scalar, too raise ValueError( "(%g, %g, %g) is not in the Weyl chamber" % (c1, c2, c3)) else: raise ValueError( "Not all values (c1, c2, c3) are in the Weyl chamber") return result
[docs]def point_in_region( region: str, c1: float, c2: float, c3: float, check_weyl=False) -> bool: """Check if $(c_1, c_2, c_3)$ are in the given region of the Weyl chamber The regions are 'W0' (between origin O and perfect entanglers polyhedron), 'W0*' (between point A1 and perfect entangler polyhedron), 'W1' (between A3 point and perfect entanglers polyhedron), and 'PE' (inside perfect entanglers polyhedron) If the check_weyl parameter is given a True, raise a ValueError for any points outside of the Weyl chamber Examples: >>> from weylchamber.visualize import WeylChamber >>> point_in_region('W0', *WeylChamber.O) True >>> point_in_region('W0', 0.2, 0.05, 0.0) True >>> point_in_region('W0', *WeylChamber.L) False >>> point_in_region('W0', *WeylChamber.Q) False >>> point_in_region('PE', *WeylChamber.Q) True >>> point_in_region('W0*', *WeylChamber.A1) True >>> point_in_region('W0*', 0.8, 0.1, 0.1) True >>> point_in_region('W1', *WeylChamber.A3) True >>> point_in_region('W1', 0.5, 0.4, 0.25) True >>> point_in_region('W1', 0.5, 0.25, 0) False >>> point_in_region('PE', 0.5, 0.25, 0) True The function may be also applied against arrays:: >>> res = point_in_region('W1', [0.5,0.5], [0.4,0.25], [0.25,0.0]) >>> assert np.all(res == np.array([ True, False])) """ regions = ['W0', 'W0*', 'W1', 'PE'] if region == 'PE': return _point_in_PE(c1, c2, c3, check_weyl=check_weyl) else: in_weyl = point_in_weyl_chamber(c1, c2, c3, raise_exception=check_weyl) c1 = np.array(c1) c2 = np.array(c2) c3 = np.array(c3) if region == 'W0': return And(in_weyl, less(c1+c2, 0.5)) elif region == 'W0*': return And(in_weyl, greater(c1-c2, 0.5)) elif region == 'W1': return And(in_weyl, greater(c2+c3, 0.5)) else: raise ValueError("region %s is not in %s" % (region, regions))
def _point_in_PE(c1, c2, c3, check_weyl=False): """Return True if the coordinates c1, c2, c3 are inside the perfect-entangler polyhedron >>> BGATE = qutip.gates.berkeley() >>> _point_in_PE(*c1c2c3(BGATE)) True >>> _point_in_PE(*c1c2c3(qutip.gates.identity(4))) False >>> res = _point_in_PE([0.0, 0.5, 0.8], [1.0, 0.25, 0.0], [1.0, 0.25, 0.0]) >>> assert np.all(res == np.array([False, True, False])) """ in_weyl = point_in_weyl_chamber(c1, c2, c3, raise_exception=check_weyl) c1 = np.array(c1); c2 = np.array(c2); c3 = np.array(c3) return And(in_weyl, And(And( greater_eq(c1+c2, 0.5), less_eq(c1-c2, 0.5)), less_eq(c2+c3, 0.5) ) )
[docs]def weyl_region(c1: float, c2: float, c3: float, check_weyl=True) -> str: """Return the region of the Weyl chamber the given point is in. Returns: One of 'W0', 'W0*', 'W1', 'PE' Examples: >>> from weylchamber.visualize import WeylChamber >>> print(weyl_region(*WeylChamber.O)) W0 >>> print(weyl_region(*WeylChamber.A1)) W0* >>> print(weyl_region(*WeylChamber.A3)) W1 >>> print(weyl_region(*WeylChamber.L)) PE >>> print(weyl_region(0.2, 0.05, 0.0)) W0 >>> print(weyl_region(0.8, 0.1, 0.1)) W0* >>> print(weyl_region(0.5, 0.25, 0)) PE >>> print(weyl_region(0.5, 0.4, 0.25)) W1 >>> try: ... weyl_region(1.0, 0.5, 0) ... except ValueError as e: ... print(e) (1, 0.5, 0) is not in the Weyl chamber >>> print(weyl_region(1.0, 0.1, 0, check_weyl=False)) W0* Only scalar values are accepted for c1, c2, c3 """ point_in_weyl_chamber(c1, c2, c3, raise_exception=check_weyl) if c1+c2 < 0.5: return 'W0' elif c1-c2 > 0.5: return 'W0*' elif c2+c3 > 0.5: return 'W1' else: return 'PE'
[docs]def to_magic(A: Gate) -> qutip.Qobj: """Convert A from the canonical basis to the the "magic" Bell basis""" if A.shape != (4, 4): raise ValueError("Gates must have a 4×4 shape") return Qmagic.dag() * qutip.Qobj(A, dims=[[2, 2], [2, 2]]) * Qmagic
[docs]def from_magic(A: Gate) -> qutip.Qobj: """ The inverse of :func:.to_magic""" if A.shape != (4, 4): raise ValueError("Gates must have a 4×4 shape") return Qmagic * qutip.Qobj(A, dims=[[2, 2], [2, 2]]) * Qmagic.dag()
[docs]def random_weyl_point(region=None) -> CTuple: """Return a random point $(c_1, c_2, c_3)$ in the Weyl chamber (units of π) If region is given in ['W0', 'W0*', 'W1', 'PE'], the point will be in the specified region of the Weyl chamber Example: >>> c1, c2, c3 = random_weyl_point() >>> point_in_weyl_chamber(c1, c2, c3) True >>> c1, c2, c3 = random_weyl_point(region='PE') >>> point_in_region('PE', c1, c2, c3) True >>> c1, c2, c3 = random_weyl_point(region='W0') >>> point_in_region('W0', c1, c2, c3) True >>> c1, c2, c3 = random_weyl_point(region='W0*') >>> point_in_region('W0*', c1, c2, c3) True >>> c1, c2, c3 = random_weyl_point(region='W1') >>> point_in_region('W1', c1, c2, c3) True """ while True: c1 = np.random.rand() c2 = 0.5*np.random.rand() c3 = 0.5*np.random.rand() if point_in_weyl_chamber(c1, c2, c3): if region is None: return c1, c2, c3 else: if point_in_region(region, c1, c2, c3): return c1, c2, c3
def _U2(phi, theta, phi1, phi2): r"""Return a unitary gate as numpy array using the parametrization .. math:: U = e^{i \phi} \begin{bmatrix} \cos\theta e^{ i\phi_1} & \sin\theta e^{ i\phi_2}\\ -\sin\theta e^{-i\phi_2} & \cos\theta e^{-i\phi_1}\\ \end{bmatrix} """ from numpy import exp, sin, cos return exp(1j*phi) * np.array([ [ cos(theta) * exp( 1j*phi1), sin(theta) * exp( 1j*phi2) ], [ -sin(theta) * exp(-1j*phi2), cos(theta) * exp(-1j*phi1) ]]) def _SQ_unitary( phi_left, theta_left, phi1_left, phi2_left, phi_right, theta_right, phi1_right, phi2_right): """Return a non-entangling two-qubit gate (a two-qubit gate locally equivalent to the identity)""" Id = np.identity(2) return qutip.Qobj( np.kron(_U2(phi_left, theta_left, phi1_left, phi2_left), Id).dot( np.kron(Id, _U2(phi_right, theta_right, phi1_right, phi2_right))), dims=[[2, 2], [2, 2]])
[docs]def canonical_gate(c1: float, c2: float, c3: float) -> qutip.Qobj: r"""Return the canonical gate for the given $(c_1, c_2, c_3)$ The canonical gate is defined as .. math:: \Op{A} = e^{i \frac{\pi}{2} \left( c_1 \Op{\sigma}_x \Op{\sigma}_x + c_2 \Op{\sigma}_y \Op{\sigma}_y + c_3 \Op{\sigma}_z \Op{\sigma}_z \right)} Example: >>> gate = qutip.Qobj( ... [[0.707107+0.000000j, 0.000000+0.000000j, 0.000000+0.000000j, 0.000000+0.707107j], ... [0.000000+0.000000j, 0.707107+0.000000j, 0.000000+0.707107j, 0.000000+0.000000j], ... [0.000000+0.000000j, 0.000000+0.707107j, 0.707107+0.000000j, 0.000000+0.000000j], ... [0.000000+0.707107j, 0.000000+0.000000j, 0.000000+0.000000j, 0.707107+0.000000j]], ... dims=[[2,2], [2,2]]) >>> U = canonical_gate(0.5,0,0) >>> assert (U - gate).norm() < 1e-6 >>> assert np.max(np.abs( ... np.array(c1c2c3(U)) - np.array([0.5, 0, 0]))) < 1e-15 """ return (np.pi*0.5j * (c1 * SxSx + c2 * SySy + c3 * SzSz)).expm()
[docs]def random_gate(region=None) -> qutip.Qobj: """Return a random two-qubit gate If region is not None, the gate will be in the specified region of the Weyl chamber. The following region names are allowed: * 'SQ': single qubit gates, consisting of the Weyl chamber points O, A2 * 'PE': polyhedron of perfect entanglers * 'W0': region between point O and the PE polyhedron * 'W0*': region between point A2 and the PE polyhedron * 'W1': region between point A3 ([SWAP]) and the PE polyhedron """ if region == 'SQ': p = 2.0 * np.pi * np.random.random(8) return _SQ_unitary(*p) else: A = canonical_gate(*random_weyl_point(region=region)) # Qobj p = 2.0 * np.pi * np.random.random(16) return _SQ_unitary(*p[:8]) * A * _SQ_unitary(*p[8:])