diff --git a/QGL/BasicSequences/RB.py b/QGL/BasicSequences/RB.py index 9a2c1e58..c7c87466 100644 --- a/QGL/BasicSequences/RB.py +++ b/QGL/BasicSequences/RB.py @@ -2,6 +2,7 @@ from ..Compiler import compile_to_hardware from ..PulseSequencePlotter import plot_pulse_files from ..Cliffords import clifford_seq, clifford_mat, inverse_clifford +from ..Euler import XYXClifford from .helpers import create_cal_seqs, cal_descriptor import os @@ -481,6 +482,67 @@ def SingleQubitRB_DiAC(qubit, plot_pulse_files(metafile) return metafile +def SingleQubitRB_XYX(qubit, seqs, purity=False, showPlot=False, add_cals=True): + """ + Single qubit randomized benchmarking using XYX Euler pulses. + + Parameters + ---------- + qubit : Channels.LogicalChannel + Logical channel to implement sequence + seqs : int iterable + list of lists of Clifford group integers produced by create_RB_seqs + purity : boolean, optional + If True, this create sequences for purity RB + showPlot : boolean, optional + Whether to plot + add_cals : boolean, optional + Whether to append calibration pulses to the end of the sequence + + Returns + ------- + metafile : string + Path to a json metafile with details about the sequences and paths + to compiled machine files + + Examples + -------- + >>> seqs = create_RB_seqs(1, [2,4,8], repeats=2, interleaveGate=1); + >>> mf = SingleQubitRB(q1, seqs); + Compiled 10 sequences. + >>> mf + '/path/to/exp/exp-meta.json' + """ + + seqsBis = [] + op = [Id(qubit, length=0), Y90m(qubit), X90(qubit)] + for ct in range(3 if purity else 1): + for seq in seqs: + seqsBis.append([XYXClifford(qubit, c) for c in seq]) + #append tomography pulse to measure purity + seqsBis[-1].append(op[ct]) + #append measurement + seqsBis[-1].append(MEAS(qubit)) + + axis_descriptor = [{ + 'name': 'length', + 'unit': None, + 'points': list(map(len, seqs)), + 'partition': 1 + }] + + #Tack on the calibration sequences + if add_cals: + seqsBis += create_cal_seqs((qubit, ), 2) + axis_descriptor.append(cal_descriptor((qubit,), 2)) + + metafile = compile_to_hardware(seqsBis, 'RB_XYX/RB_XYX', axis_descriptor = axis_descriptor, extra_meta = {'sequences':seqs}) + + if showPlot: + plot_pulse_files(metafile) + return metafile + + def SingleQubitIRB_AC(qubit, seqFile, showPlot=False): """ Single qubit interleaved randomized benchmarking using atomic Clifford diff --git a/QGL/Euler.py b/QGL/Euler.py new file mode 100644 index 00000000..1ebbd359 --- /dev/null +++ b/QGL/Euler.py @@ -0,0 +1,195 @@ +""" +Tools for creating Euler-angle based gates. + +Original Author: Guilhem Ribeill, Luke Govia + +Copyright 2020 Raytheon BBN Technologies + +Licensed under the Apache License, Version 2.0 (the "License"); +you may not use this file except in compliance with the License. +You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + +Unless required by applicable law or agreed to in writing, software +distributed under the License is distributed on an "AS IS" BASIS, +WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +See the License for the specific language governing permissions and +limitations under the License. +""" + +import numpy as np +from scipy.linalg import expm +from .Cliffords import C1 +from .PulsePrimitives import * + +#Pauli matrices +pX = np.array([[0, 1], [1, 0]], dtype=np.complex128) +pZ = np.array([[1, 0], [0, -1]], dtype=np.complex128) +pY = 1j * pX @ pZ +pI = np.eye(2, dtype=np.complex128) + +#Machine precision +_eps = np.finfo(np.complex128).eps + +#### FUNCTIONS COPIED FROM PYGSTI +#### See: https://github.com/pyGSTio/pyGSTi + +#### PYGSTI NOTICE +# Python GST Implementation (PyGSTi) v. 0.9 +# Copyright 2015, 2019 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +# Under the terms of Contract DE-NA0003525 with NTESS, the U.S. Government retains certain rights in this software. +#### END PYGSTI NOTICE + +#### PYGSTI COPYRRIGHT +# Copyright 2015, 2019 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +# Under the terms of Contract DE-NA0003525 with NTESS, the U.S. Government retains certain rights +# in this software. +# Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except +# in compliance with the License. You may obtain a copy of the License at +# http://www.apache.org/licenses/LICENSE-2.0 or in the LICENSE file in the root pyGSTi directory. +#### END PYGSTI COPYRIGHT + +def tracenorm(A, tol=np.sqrt(_eps)): + """Compute the trace norm of a matrix A given by: + Tr(sqrt{A^dagger * A}) + From: https://github.com/pyGSTio/pyGSTi/blob/master/pygsti/tools/optools.py + """ + if np.linalg.norm(A - np.conjugate(A.T)) < tol: + #Hermitian, so just sum eigenvalue magnitudes + return np.sum(np.abs(np.linalg.eigvals(A))) + else: + #Sum of singular values (positive by construction) + return np.sum(np.linalg.svd(A, compute_uv=False)) + +def tracedist(A, B, tol=np.sqrt(_eps)): + """Compute the trace distance between matrices A and B given by: + 0.5 * Tr(sqrt{(A-B)^dagger * (A-B)}) + From: https://github.com/pyGSTio/pyGSTi/blob/master/pygsti/tools/optools.py + """ + return 0.5 * tracenorm(A - B) + +#### END FUNCTIONS COPIED FROM PYGSTI + +def is_close(A, B, tol=np.sqrt(_eps)): + """Check if two matrices are close in the sense of trace distance. + """ + if tracedist(A, B) < tol: + return True + else: + A /= np.exp(1j*np.angle(A[0,0])) + B /= np.exp(1j*np.angle(B[0,0])) + return tracedist(A, B) < tol + +def haar_unitary(d): + """Generate a Haar-random unitary matrix of dimension d. + Algorithm from: + F. Medrazzi. "How to generate random matrices from the classical compact groups" + arXiv: math-ph/0609050 + """ + assert d > 1, 'Dimension must be > 1!' + re_Z = np.random.randn(d*d).reshape((d,d)) + im_Z = np.random.randn(d*d).reshape((d,d)) + Z = (re_Z + 1j*im_Z)/np.sqrt(2.0) + Q, R = np.linalg.qr(Z) + L = np.diag(np.diag(R) / np.abs(np.diag(R))) + return Q @ L @ Q + +def xyx_unitary(α, β, γ): + """Unitary decomposed as Rx, Ry, Rx rotations. + Angles are in matrix order, not in circuit order! + """ + return expm(-0.5j*α*pX)@expm(-0.5j*β*pY)@expm(-0.5j*γ*pX) + +def zyz_unitary(ϕ, θ, λ): + """Unitary decomposed as Rz, Ry, Rz rotations. + Angles are in matrix order, not in circuit order! + """ + return expm(-0.5j*ϕ*pZ)@expm(-0.5j*θ*pY)@expm(-0.5j*λ*pZ) + +def diatomic_unitary(a, b, c): + """Unitary decomposed as a diatomic gate of the form + Ztheta + X90 + Ztheta + X90 + Ztheta + """ + X90 = expm(-0.25j*np.pi*pX) + return expm(-0.5j*a*pZ)@X90@expm(-0.5j*b*pZ)@X90@expm(-0.5j*c*pZ) + +def zyz_angles(U): + """Euler angles for a unitary matrix U in the sequence Z-Y-Z. + Note that angles are returned in matrix multiplication, not circuit order. + """ + assert U.shape == (2,2), "Must use a 2x2 matrix!" + k = 1.0/np.sqrt(np.linalg.det(U)) + SU = k*U + θ = 2 * np.arctan2(np.abs(SU[1,0]), np.abs(SU[0,0])) + a = 2 * np.angle(SU[1,1]) + b = 2 * np.angle(SU[1,0]) + ϕ = (a + b) * 0.5 + λ = (a - b) * 0.5 + return (ϕ, θ, λ) + +def xyx_angles(U): + """Euler angles for a unitary matrix U in the sequence X-Y-X. + Note that angles are returned in matrix multiplication, not circuit order. + We make use of the identity: + Rx(a)Ry(b)Rx(c) = H Rz(a) Ry(-b) Rz(c) H + """ + H = np.array([[1., 1.], [1., -1.]], dtype=np.complex128)/np.sqrt(2) + ϕ, θ, λ = zyz_angles(H@U@H) + return (ϕ, -1.0*θ, λ) + +def diatomic_angles(U): + ϕ, θ, λ = zyz_angles(U) + a = ϕ + b = np.pi - θ + c = λ - np.pi + return (a, b, c) + +# C1 = {} +# C1[0] = pI +# C1[1] = expm(-1j * (pi / 4) * pX) +# C1[2] = expm(-2j * (pi / 4) * pX) +# C1[3] = expm(-3j * (pi / 4) * pX) +# C1[4] = expm(-1j * (pi / 4) * pY) +# C1[5] = expm(-2j * (pi / 4) * pY) +# C1[6] = expm(-3j * (pi / 4) * pY) +# C1[7] = expm(-1j * (pi / 4) * pZ) +# C1[8] = expm(-2j * (pi / 4) * pZ) +# C1[9] = expm(-3j * (pi / 4) * pZ) +# C1[10] = expm(-1j * (pi / 2) * (1 / np.sqrt(2)) * (pX + pY)) +# C1[11] = expm(-1j * (pi / 2) * (1 / np.sqrt(2)) * (pX - pY)) +# C1[12] = expm(-1j * (pi / 2) * (1 / np.sqrt(2)) * (pX + pZ)) +# C1[13] = expm(-1j * (pi / 2) * (1 / np.sqrt(2)) * (pX - pZ)) +# C1[14] = expm(-1j * (pi / 2) * (1 / np.sqrt(2)) * (pY + pZ)) +# C1[15] = expm(-1j * (pi / 2) * (1 / np.sqrt(2)) * (pY - pZ)) +# C1[16] = expm(-1j * (pi / 3) * (1 / np.sqrt(3)) * (pX + pY + pZ)) +# C1[17] = expm(-2j * (pi / 3) * (1 / np.sqrt(3)) * (pX + pY + pZ)) +# C1[18] = expm(-1j * (pi / 3) * (1 / np.sqrt(3)) * (pX - pY + pZ)) +# C1[19] = expm(-2j * (pi / 3) * (1 / np.sqrt(3)) * (pX - pY + pZ)) +# C1[20] = expm(-1j * (pi / 3) * (1 / np.sqrt(3)) * (pX + pY - pZ)) +# C1[21] = expm(-2j * (pi / 3) * (1 / np.sqrt(3)) * (pX + pY - pZ)) +# C1[22] = expm(-1j * (pi / 3) * (1 / np.sqrt(3)) * (-pX + pY + pZ)) +# C1[23] = expm(-2j * (pi / 3) * (1 / np.sqrt(3)) * (-pX + pY + pZ)) + +def XYXClifford(qubit, cliff_num): + """ + The set of 24 Diatomic Clifford single qubit pulses. Each pulse is decomposed + as Rx(α)Ry(β)Rx(γ). + + Parameters + ---------- + qubit : logical channel to implement sequence (LogicalChannel) + cliffNum : the zero-indexed Clifford number + + Returns + ------- + pulse object + """ + α, β, γ = xyx_angles(C1[cliff_num]) + + p1 = Id(qubit) if np.isclose(γ, 0.0) else Xtheta(qubit, angle=γ) + p2 = Id(qubit) if np.isclose(β, 0.0) else Ytheta(qubit, angle=β) + p3 = Id(qubit) if np.isclose(α, 0.0) else Xtheta(qubit, angle=α) + + return p1 + p2 + p3 + diff --git a/tests/test_euler.py b/tests/test_euler.py new file mode 100644 index 00000000..b51757b9 --- /dev/null +++ b/tests/test_euler.py @@ -0,0 +1,50 @@ +import unittest + +from QGL import * +from QGL.Euler import * +from QGL.Cliffords import C1 +import QGL.config +try: + from helpers import setup_test_lib +except: + from .helpers import setup_test_lib + +class EulerDecompositions(unittest.TestCase): + + N_test = 1000 + + def setUp(self): + pass + #setup_test_lib() + #self.q1 = QubitFactory('q1') + + def test_zyz_decomp(self): + for j in range(self.N_test): + Uh = haar_unitary(2) + Ux = zyz_unitary(*zyz_angles(Uh)) + assert is_close(Uh, Ux) + + def test_xyx_decomp(self): + for j in range(self.N_test): + Uh = haar_unitary(2) + Ux = xyx_unitary(*xyx_angles(Uh)) + assert is_close(Uh, Ux) + + def test_diatomic_decomp(self): + for j in range(self.N_test): + Uh = haar_unitary(2) + Ux = diatomic_unitary(*diatomic_angles(Uh)) + assert is_close(Uh, Ux) + + def test_xyx_cliffords(self): + for j in range(24): + Uxyx = xyx_unitary(*xyx_angles(C1[j])) + assert is_close(Uxyx, C1[j]), f"{j}" + + def test_diatomic_cliffords(self): + for j in range(24): + Ud = diatomic_unitary(*diatomic_angles(C1[j])) + assert is_close(Ud, C1[j]), f"{j}" + +if __name__ == "__main__": + unittest.main()