From 01d39f6ade58a3ecf84fd564d35a5efb4bbe1cff Mon Sep 17 00:00:00 2001 From: Guilhem Ribeill Date: Mon, 13 Apr 2020 17:31:32 -0400 Subject: [PATCH 01/13] Add euler angle conversions --- QGL/Euler.py | 178 +++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 178 insertions(+) create mode 100644 QGL/Euler.py diff --git a/QGL/Euler.py b/QGL/Euler.py new file mode 100644 index 00000000..41e801a3 --- /dev/null +++ b/QGL/Euler.py @@ -0,0 +1,178 @@ +""" +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 .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(ϕ, θ, λ): + return expm(-0.5j*ϕ*pZ)@expm(-0.5j*θ*pY)@expm(-0.5j*λ*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*θ, λ) + +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 + From ae7f22b5fc606a8415e286ffabd43d3fbb71977c Mon Sep 17 00:00:00 2001 From: Guilhem Ribeill Date: Mon, 13 Apr 2020 17:31:54 -0400 Subject: [PATCH 02/13] Add Euler angle conversion tests. --- tests/test_euler.py | 39 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 39 insertions(+) create mode 100644 tests/test_euler.py diff --git a/tests/test_euler.py b/tests/test_euler.py new file mode 100644 index 00000000..cb937a68 --- /dev/null +++ b/tests/test_euler.py @@ -0,0 +1,39 @@ +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_xyx_cliffords(self): + for j in range(24): + Uxyx = xyx_unitary(*xyx_angles(C1[j])) + assert is_close(Uxyx, C1[j]), f"{j}" + +if __name__ == "__main__": + unittest.main() From c0bd5056eed0657ff8c85bc95a2547aa27786369 Mon Sep 17 00:00:00 2001 From: Guilhem Ribeill Date: Mon, 13 Apr 2020 17:37:08 -0400 Subject: [PATCH 03/13] Add sequence helper function --- QGL/BasicSequences/RB.py | 62 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 62 insertions(+) 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 From 83fba4e8d6933e51a4333fc02dee5f3ce4142730 Mon Sep 17 00:00:00 2001 From: Guilhem Ribeill Date: Tue, 14 Apr 2020 12:09:40 -0400 Subject: [PATCH 04/13] Add diatomic Euler angle decomposition --- QGL/Euler.py | 69 ++++++++++++++++++++++++++++----------------- tests/test_euler.py | 11 ++++++++ 2 files changed, 54 insertions(+), 26 deletions(-) diff --git a/QGL/Euler.py b/QGL/Euler.py index 41e801a3..1ebbd359 100644 --- a/QGL/Euler.py +++ b/QGL/Euler.py @@ -20,7 +20,7 @@ import numpy as np from scipy.linalg import expm - +from .Cliffords import C1 from .PulsePrimitives import * #Pauli matrices @@ -102,8 +102,18 @@ def xyx_unitary(α, β, γ): 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. @@ -128,31 +138,38 @@ def xyx_angles(U): ϕ, θ, λ = zyz_angles(H@U@H) return (ϕ, -1.0*θ, λ) -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 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): """ diff --git a/tests/test_euler.py b/tests/test_euler.py index cb937a68..b51757b9 100644 --- a/tests/test_euler.py +++ b/tests/test_euler.py @@ -30,10 +30,21 @@ def test_xyx_decomp(self): 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() From 1d072678152b15d23d2aca6c65ee5241d265e4a9 Mon Sep 17 00:00:00 2001 From: Guilhem Ribeill Date: Tue, 14 Apr 2020 15:57:25 -0400 Subject: [PATCH 05/13] Make sure XYX angles are in (-pi, pi). --- QGL/Euler.py | 17 +++++++++++++---- 1 file changed, 13 insertions(+), 4 deletions(-) diff --git a/QGL/Euler.py b/QGL/Euler.py index 1ebbd359..6476345b 100644 --- a/QGL/Euler.py +++ b/QGL/Euler.py @@ -72,14 +72,16 @@ def tracedist(A, B, tol=np.sqrt(_eps)): #### 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. + """Check if two matrices are close in the sense of trace distance. """ if tracedist(A, B) < tol: - return True + return True else: + A[np.abs(A) < tol] = 0.0 + B[np.abs(B) < tol] = 0.0 A /= np.exp(1j*np.angle(A[0,0])) B /= np.exp(1j*np.angle(B[0,0])) - return tracedist(A, B) < tol + return ((tracedist(A, B) < tol) or (tracedist(A, -1.0*B) < tol)) def haar_unitary(d): """Generate a Haar-random unitary matrix of dimension d. @@ -128,6 +130,13 @@ def zyz_angles(U): λ = (a - b) * 0.5 return (ϕ, θ, λ) +def _mod_2pi(angle): + if angle > np.pi: + angle -= 2*np.pi + if angle < -np.pi: + angle += 2*np.pi + return angle + 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. @@ -136,7 +145,7 @@ def xyx_angles(U): """ H = np.array([[1., 1.], [1., -1.]], dtype=np.complex128)/np.sqrt(2) ϕ, θ, λ = zyz_angles(H@U@H) - return (ϕ, -1.0*θ, λ) + return (_mod_2pi(ϕ), _mod_2pi(-1.0*θ), _mod_2pi(λ)) def diatomic_angles(U): ϕ, θ, λ = zyz_angles(U) From a58a74e0152f3a8c0c648598834ffb4953cf993f Mon Sep 17 00:00:00 2001 From: Guilhem Ribeill Date: Wed, 15 Apr 2020 17:05:08 -0400 Subject: [PATCH 06/13] Move matrix tools to separate file. --- QGL/Euler.py | 101 +------------------------------------------- tests/test_euler.py | 1 + 2 files changed, 2 insertions(+), 100 deletions(-) diff --git a/QGL/Euler.py b/QGL/Euler.py index 6476345b..e878a9c9 100644 --- a/QGL/Euler.py +++ b/QGL/Euler.py @@ -22,80 +22,7 @@ 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.abs(A) < tol] = 0.0 - B[np.abs(B) < tol] = 0.0 - A /= np.exp(1j*np.angle(A[0,0])) - B /= np.exp(1j*np.angle(B[0,0])) - return ((tracedist(A, B) < tol) or (tracedist(A, -1.0*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 +from .tools.matrix_tools import * def xyx_unitary(α, β, γ): """Unitary decomposed as Rx, Ry, Rx rotations. @@ -154,32 +81,6 @@ def diatomic_angles(U): 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 diff --git a/tests/test_euler.py b/tests/test_euler.py index b51757b9..a19eb4c0 100644 --- a/tests/test_euler.py +++ b/tests/test_euler.py @@ -2,6 +2,7 @@ from QGL import * from QGL.Euler import * +from QGL.tools.matrix_tools import * from QGL.Cliffords import C1 import QGL.config try: From 04b1d6669a2bf4c2ffa01ddcae444b47c77ec5da Mon Sep 17 00:00:00 2001 From: Guilhem Ribeill Date: Wed, 15 Apr 2020 17:18:22 -0400 Subject: [PATCH 07/13] Refactor Clifford operations to separate out those that deal with pulses from pure matrix manipulation --- QGL/Cliffords.py | 176 +------------------------------ QGL/tools/clifford_tools.py | 200 ++++++++++++++++++++++++++++++++++++ QGL/tools/matrix_tools.py | 95 +++++++++++++++++ 3 files changed, 296 insertions(+), 175 deletions(-) create mode 100644 QGL/tools/clifford_tools.py create mode 100644 QGL/tools/matrix_tools.py diff --git a/QGL/Cliffords.py b/QGL/Cliffords.py index 6a050fde..932c7c77 100644 --- a/QGL/Cliffords.py +++ b/QGL/Cliffords.py @@ -9,84 +9,9 @@ import operator from functools import reduce +from .tools.clifford_tools import * from .PulsePrimitives import * -#Single qubit paulis -pX = np.array([[0, 1], [1, 0]], dtype=np.complex128) -pY = np.array([[0, -1j], [1j, 0]], dtype=np.complex128) -pZ = np.array([[1, 0], [0, -1]], dtype=np.complex128) -pI = np.eye(2, dtype=np.complex128) - - -def pauli_mats(n): - """ - Return a list of n-qubit Paulis as numpy array. - """ - assert n > 0, "You need at least 1 qubit!" - if n == 1: - return [pI, pX, pY, pZ] - else: - paulis = pauli_mats(n - 1) - return [np.kron(p1, p2) - for p1, p2 in product([pI, pX, pY, pZ], paulis)] - -#Basis single-qubit Cliffords with an arbitrary enumeration -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)) - - -#A little memoize decorator -def memoize(function): - cache = {} - - def decorated(*args): - if args not in cache: - cache[args] = function(*args) - return cache[args] - - return decorated - - -@memoize -def clifford_multiply(c1, c2): - ''' - Multiplication table for single qubit cliffords. Note this assumes c1 is applied first. - i.e. clifford_multiply(c1, c2) calculates c2*c1 - ''' - tmpMult = np.dot(C1[c2], C1[c1]) - checkArray = np.array( - [np.abs(np.trace(np.dot(tmpMult.transpose().conj(), C1[x]))) - for x in range(24)]) - return checkArray.argmax() - -#We can usually (without atomic Cliffords) only apply a subset of the single-qubit Cliffords -#i.e. the pulses that we can apply: Id, X90, X90m, Y90, Y90m, X, Y -generatorPulses = [0, 1, 3, 4, 6, 2, 5] - - def generator_pulse(G): """ A function that returns the pulse corresponding to a generator @@ -101,66 +26,6 @@ def generator_pulse(G): 5: (Y, Ym)} return choice(generatorPulses[G]) -#Get all combinations of generator sequences up to length three -generatorSeqs = [x for x in product(generatorPulses,repeat=1)] + \ - [x for x in product(generatorPulses,repeat=2)] + \ - [x for x in product(generatorPulses,repeat=3)] - -#Find the effective unitary for each generator sequence -reducedSeqs = np.array([reduce(clifford_multiply, x) for x in generatorSeqs]) - -#Pick first generator sequence (and thus shortest) that gives each Clifford and then -#also add all those that have the same length - -#First for each of the 24 single-qubit Cliffords find which sequences create them -allC1Seqs = [np.nonzero(reducedSeqs == x)[0] for x in range(24)] -#And the length of the first one for all 24 -minSeqLengths = [len(generatorSeqs[seqs[0]]) for seqs in allC1Seqs] -#Now pull out all those that are the same length as the first one -C1Seqs = [] -for minLength, seqs in zip(minSeqLengths, allC1Seqs): - C1Seqs.append([s for s in seqs if len(generatorSeqs[s]) == minLength]) - -C2Seqs = [] -""" -The IBM paper has the Sgroup (rotation n*(pi/3) rotations about the X+Y+Z axis) -Sgroup = [C[0], C[16], C[17]] - -The two qubit Cliffords can be written down as the product of -1. A choice of one of 24^2 C \otimes C single-qubit Cliffords -2. Optionally an entangling gate from CNOT, iSWAP and SWAP -3. Optional one of 9 S \otimes S gate - -Therefore, we'll enumerate the two-qubit Clifford as a three tuple ((c1,c2), Entangling, (s1,s2)) - -""" - -#1. All pairs of single-qubit Cliffords -for c1, c2 in product(range(24), repeat=2): - C2Seqs.append(((c1, c2), None, None)) - -#2. The CNOT-like class, replacing the CNOT with a echoCR -#TODO: sort out whether we need to explicitly encorporate the single qubit rotations into the trailing S gates -# The leading single-qubit Cliffords are fully sampled so they should be fine - -for (c1, c2), (s1, s2) in product( - product( - range(24), repeat=2), - product([0, 16, 17], repeat=2)): - C2Seqs.append(((c1, c2), "CNOT", (s1, s2))) - -#3. iSWAP like class - replacing iSWAP with (echoCR - (Y90m*Y90m) - echoCR) -for (c1, c2), (s1, s2) in product( - product( - range(24), repeat=2), - product([0, 16, 17], repeat=2)): - C2Seqs.append(((c1, c2), "iSWAP", (s1, s2))) - -#4. SWAP like class -for c1, c2 in product(range(24), repeat=2): - C2Seqs.append(((c1, c2), "SWAP", None)) - - def Cx2(c1, c2, q1, q2): """ Helper function to create pulse block for a pair of single-qubit Cliffords @@ -221,42 +86,3 @@ def clifford_seq(c, q1, q2=None): if c[2]: seq += [Cx2(c[2][0], c[2][1], q1, q2)] return seq - - -@memoize -def clifford_mat(c, numQubits): - """ - Return the matrix unitary the implements the qubit clifford C - """ - assert numQubits <= 2, "Oops! I only handle one or two qubits" - if numQubits == 1: - return C1[c] - else: - c = C2Seqs[c] - mat = np.kron(clifford_mat(c[0][0], 1), clifford_mat(c[0][1], 1)) - if c[1]: - mat = np.dot(entangling_mat(c[1]), mat) - if c[2]: - mat = np.dot( - np.kron( - clifford_mat(c[2][0], 1), clifford_mat(c[2][1], 1)), mat) - return mat - - -def inverse_clifford(cMat): - dim = cMat.shape[0] - if dim == 2: - for ct in range(24): - if np.isclose( - np.abs(np.dot(cMat, clifford_mat(ct, 1)).trace()), dim): - return ct - elif dim == 4: - for ct in range(len(C2Seqs)): - if np.isclose( - np.abs(np.dot(cMat, clifford_mat(ct, 2)).trace()), dim): - return ct - else: - raise Exception("Expected 2 or 4 qubit dimensional matrix.") - - #If we got here something is wrong - raise Exception("Couldn't find inverse clifford") diff --git a/QGL/tools/clifford_tools.py b/QGL/tools/clifford_tools.py new file mode 100644 index 00000000..566372de --- /dev/null +++ b/QGL/tools/clifford_tools.py @@ -0,0 +1,200 @@ +""" +Tools for manipulating 1 and 2 qubit cliffords. + +Original Author: Blake Johnson, Colm Ryan, Guilhem Ribeill + +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 numpy import pi +from itertools import product +from random import choice +import operator +from functools import reduce + +#Single qubit paulis +pX = np.array([[0, 1], [1, 0]], dtype=np.complex128) +pY = np.array([[0, -1j], [1j, 0]], dtype=np.complex128) +pZ = np.array([[1, 0], [0, -1]], dtype=np.complex128) +pI = np.eye(2, dtype=np.complex128) + +def pauli_mats(n): + """ + Return a list of n-qubit Paulis as numpy array. + """ + assert n > 0, "You need at least 1 qubit!" + if n == 1: + return [pI, pX, pY, pZ] + else: + paulis = pauli_mats(n - 1) + return [np.kron(p1, p2) + for p1, p2 in product([pI, pX, pY, pZ], paulis)] + +#Basis single-qubit Cliffords with an arbitrary enumeration +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)) + + +#A little memoize decorator +def memoize(function): + cache = {} + + def decorated(*args): + if args not in cache: + cache[args] = function(*args) + return cache[args] + + return decorated + + +@memoize +def clifford_multiply(c1, c2): + ''' + Multiplication table for single qubit cliffords. Note this assumes c1 is applied first. + i.e. clifford_multiply(c1, c2) calculates c2*c1 + ''' + tmpMult = np.dot(C1[c2], C1[c1]) + checkArray = np.array( + [np.abs(np.trace(np.dot(tmpMult.transpose().conj(), C1[x]))) + for x in range(24)]) + return checkArray.argmax() + +#We can usually (without atomic Cliffords) only apply a subset of the single-qubit Cliffords +#i.e. the pulses that we can apply: Id, X90, X90m, Y90, Y90m, X, Y +generatorPulses = [0, 1, 3, 4, 6, 2, 5] + +#Get all combinations of generator sequences up to length three +generatorSeqs = [x for x in product(generatorPulses,repeat=1)] + \ + [x for x in product(generatorPulses,repeat=2)] + \ + [x for x in product(generatorPulses,repeat=3)] + +#Find the effective unitary for each generator sequence +reducedSeqs = np.array([reduce(clifford_multiply, x) for x in generatorSeqs]) + +#Pick first generator sequence (and thus shortest) that gives each Clifford and then +#also add all those that have the same length + +#First for each of the 24 single-qubit Cliffords find which sequences create them +allC1Seqs = [np.nonzero(reducedSeqs == x)[0] for x in range(24)] +#And the length of the first one for all 24 +minSeqLengths = [len(generatorSeqs[seqs[0]]) for seqs in allC1Seqs] +#Now pull out all those that are the same length as the first one +C1Seqs = [] +for minLength, seqs in zip(minSeqLengths, allC1Seqs): + C1Seqs.append([s for s in seqs if len(generatorSeqs[s]) == minLength]) + +C2Seqs = [] +""" +The IBM paper has the Sgroup (rotation n*(pi/3) rotations about the X+Y+Z axis) +Sgroup = [C[0], C[16], C[17]] + +The two qubit Cliffords can be written down as the product of +1. A choice of one of 24^2 C \otimes C single-qubit Cliffords +2. Optionally an entangling gate from CNOT, iSWAP and SWAP +3. Optional one of 9 S \otimes S gate + +Therefore, we'll enumerate the two-qubit Clifford as a three tuple ((c1,c2), Entangling, (s1,s2)) + +""" + +#1. All pairs of single-qubit Cliffords +for c1, c2 in product(range(24), repeat=2): + C2Seqs.append(((c1, c2), None, None)) + +#2. The CNOT-like class, replacing the CNOT with a echoCR +#TODO: sort out whether we need to explicitly encorporate the single qubit rotations into the trailing S gates +# The leading single-qubit Cliffords are fully sampled so they should be fine + +for (c1, c2), (s1, s2) in product( + product( + range(24), repeat=2), + product([0, 16, 17], repeat=2)): + C2Seqs.append(((c1, c2), "CNOT", (s1, s2))) + +#3. iSWAP like class - replacing iSWAP with (echoCR - (Y90m*Y90m) - echoCR) +for (c1, c2), (s1, s2) in product( + product( + range(24), repeat=2), + product([0, 16, 17], repeat=2)): + C2Seqs.append(((c1, c2), "iSWAP", (s1, s2))) + +#4. SWAP like class +for c1, c2 in product(range(24), repeat=2): + C2Seqs.append(((c1, c2), "SWAP", None)) + + +@memoize +def clifford_mat(c, numQubits): + """ + Return the matrix unitary the implements the qubit clifford C + """ + assert numQubits <= 2, "Oops! I only handle one or two qubits" + if numQubits == 1: + return C1[c] + else: + c = C2Seqs[c] + mat = np.kron(clifford_mat(c[0][0], 1), clifford_mat(c[0][1], 1)) + if c[1]: + mat = np.dot(entangling_mat(c[1]), mat) + if c[2]: + mat = np.dot( + np.kron( + clifford_mat(c[2][0], 1), clifford_mat(c[2][1], 1)), mat) + return mat + + +def inverse_clifford(cMat): + """Return the inverse clifford index.""" + dim = cMat.shape[0] + if dim == 2: + for ct in range(24): + if np.isclose( + np.abs(np.dot(cMat, clifford_mat(ct, 1)).trace()), dim): + return ct + elif dim == 4: + for ct in range(len(C2Seqs)): + if np.isclose( + np.abs(np.dot(cMat, clifford_mat(ct, 2)).trace()), dim): + return ct + else: + raise Exception("Expected 2 or 4 qubit dimensional matrix.") + + #If we got here something is wrong + raise Exception("Couldn't find inverse clifford") diff --git a/QGL/tools/matrix_tools.py b/QGL/tools/matrix_tools.py new file mode 100644 index 00000000..a994231d --- /dev/null +++ b/QGL/tools/matrix_tools.py @@ -0,0 +1,95 @@ +""" +Tools for manipulating matrices + +Original Author: Guilhem Ribeill + +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 + +#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.abs(A) < tol] = 0.0 + B[np.abs(B) < tol] = 0.0 + A /= np.exp(1j*np.angle(A[0,0])) + B /= np.exp(1j*np.angle(B[0,0])) + return ((tracedist(A, B) < tol) or (tracedist(A, -1.0*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 \ No newline at end of file From f42319ad27e6cb70df7ba8a28c503d30acfafbbc Mon Sep 17 00:00:00 2001 From: Guilhem Ribeill Date: Wed, 15 Apr 2020 17:24:52 -0400 Subject: [PATCH 08/13] Move Euler angle stuff into tools directory. --- QGL/{Euler.py => tools/euler_angles.py} | 25 +------------------------ 1 file changed, 1 insertion(+), 24 deletions(-) rename QGL/{Euler.py => tools/euler_angles.py} (79%) diff --git a/QGL/Euler.py b/QGL/tools/euler_angles.py similarity index 79% rename from QGL/Euler.py rename to QGL/tools/euler_angles.py index e878a9c9..c9b197db 100644 --- a/QGL/Euler.py +++ b/QGL/tools/euler_angles.py @@ -79,27 +79,4 @@ def diatomic_angles(U): a = ϕ b = np.pi - θ c = λ - np.pi - return (a, b, c) - -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 - + return (a, b, c) \ No newline at end of file From 9b066bd0c506aa3052cc08443c297647c1df19c5 Mon Sep 17 00:00:00 2001 From: Guilhem Ribeill Date: Wed, 15 Apr 2020 17:53:47 -0400 Subject: [PATCH 09/13] WIP Refactor cliffords sequences --- QGL/BasicSequences/RB.py | 4 +- QGL/Cliffords.py | 406 ++++++++++++++++++++++++++++++++---- QGL/tools/clifford_tools.py | 18 ++ QGL/tools/euler_angles.py | 5 +- 4 files changed, 384 insertions(+), 49 deletions(-) diff --git a/QGL/BasicSequences/RB.py b/QGL/BasicSequences/RB.py index c7c87466..f3393f3f 100644 --- a/QGL/BasicSequences/RB.py +++ b/QGL/BasicSequences/RB.py @@ -1,8 +1,8 @@ from ..PulsePrimitives import * 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 ..tools.clifford_tools import clifford_seq, clifford_mat, inverse_clifford +#from ..Euler import XYXClifford from .helpers import create_cal_seqs, cal_descriptor import os diff --git a/QGL/Cliffords.py b/QGL/Cliffords.py index 932c7c77..5cf85e9c 100644 --- a/QGL/Cliffords.py +++ b/QGL/Cliffords.py @@ -10,8 +10,13 @@ from functools import reduce from .tools.clifford_tools import * +from .tools.euler_angles import xyx_angles from .PulsePrimitives import * +### +### Single Qubit Cliffords +### + def generator_pulse(G): """ A function that returns the pulse corresponding to a generator @@ -26,22 +31,359 @@ def generator_pulse(G): 5: (Y, Ym)} return choice(generatorPulses[G]) -def Cx2(c1, c2, q1, q2): +def StdClifford(qubit, cliffNum): """ - Helper function to create pulse block for a pair of single-qubit Cliffords - """ - #Create list of pulse objects on the qubits - seq1 = clifford_seq(c1, q1) - seq2 = clifford_seq(c2, q2) - #Create the pulse block - return reduce(operator.add, seq1) * reduce(operator.add, seq2) + The set of 24 clifford pulses from the Pauli group generators. + Note that this will randomly select between +/- X and +- Y pulses. + + Parameters + ---------- + qubit : logical channel to implement sequence (LogicalChannel) + cliffNum : the zero-indexed Clifford number + + Returns + ------- + pulse object + """ + genSeq = generatorSeqs[choice(C1Seqs[cliffNum])] + return reduce(operator.add, [generator_pulse(g)(qubit) for g in genSeq]) + +def AC(qubit, cliffNum): + """ + + The set of 24 Atomic Clifford single qubit pulses. + + Parameters + ---------- + qubit : logical channel to implement sequence (LogicalChannel) + cliffNum : the zero-indexed Clifford number + + Returns + ------- + pulse object + """ + + #Figure out the approximate nutation frequency calibration from the X180 and the sampling_rate + Xp = X(qubit) + xpulse = Xp.amp * Xp.shape + nutFreq = 0.5 / (sum(xpulse) / qubit.phys_chan.sampling_rate) + #Now a big else if chain for to get the specific Clifford + if cliffNum == 0: + #Identity gate + return Id(qubit, length=0) + elif cliffNum == 1: + #X90 + return X90(qubit) + elif cliffNum == 2: + #X180 + return X(qubit) + elif cliffNum == 3: + #X90m + return X90m(qubit) + elif cliffNum == 4: + #Y90 + return Y90(qubit) + elif cliffNum == 5: + #Y180 + return Y(qubit) + elif cliffNum == 6: + #Y90m + return Y90m(qubit) + elif cliffNum == 7: + #Z90 + return Z90(qubit) + elif cliffNum == 8: + #Z180 + return Z(qubit) + elif cliffNum == 9: + #Z90m + return Z90m(qubit) + elif cliffNum == 10: + #X+Y 180 + return U(qubit, phase=pi / 4, label="AC_10") + elif cliffNum == 11: + #X-Y 180 + return U(qubit, phase=-pi / 4, label="AC_11") + elif cliffNum == 12: + #X+Z 180(Hadamard) + return arb_axis_drag(qubit, + nutFreq, + rotAngle=pi, + polarAngle=pi / 4, + label="AC_12") + elif cliffNum == 13: + #X-Z 180 + return arb_axis_drag(qubit, + nutFreq, + rotAngle=pi, + polarAngle=pi / 4, + aziAngle=pi, + label="AC_13") + elif cliffNum == 14: + #Y+Z 180 + return arb_axis_drag(qubit, + nutFreq, + rotAngle=pi, + polarAngle=pi / 4, + aziAngle=pi / 2, + label="AC_14") + elif cliffNum == 15: + #Y-Z 180 + return arb_axis_drag(qubit, + nutFreq, + rotAngle=pi, + polarAngle=pi / 4, + aziAngle=-pi / 2, + label="AC_15") + elif cliffNum == 16: + #X+Y+Z 120 + return arb_axis_drag(qubit, + nutFreq, + rotAngle=2 * pi / 3, + polarAngle=acos(1 / sqrt(3)), + aziAngle=pi / 4, + label="AC_16") + elif cliffNum == 17: + #X+Y+Z -120 (equivalent to -X-Y-Z 120) + return arb_axis_drag(qubit, + nutFreq, + rotAngle=2 * pi / 3, + polarAngle=pi - acos(1 / sqrt(3)), + aziAngle=5 * pi / 4, + label="AC_17") + elif cliffNum == 18: + #X-Y+Z 120 + return arb_axis_drag(qubit, + nutFreq, + rotAngle=2 * pi / 3, + polarAngle=acos(1 / sqrt(3)), + aziAngle=-pi / 4, + label="AC_18") + elif cliffNum == 19: + #X-Y+Z 120 (equivalent to -X+Y-Z -120) + return arb_axis_drag(qubit, + nutFreq, + rotAngle=2 * pi / 3, + polarAngle=pi - acos(1 / sqrt(3)), + aziAngle=3 * pi / 4, + label="AC_19") + elif cliffNum == 20: + #X+Y-Z 120 + return arb_axis_drag(qubit, + nutFreq, + rotAngle=2 * pi / 3, + polarAngle=pi - acos(1 / sqrt(3)), + aziAngle=pi / 4, + label="AC_20") + elif cliffNum == 21: + #X+Y-Z -120 (equivalent to -X-Y+Z 120) + return arb_axis_drag(qubit, + nutFreq, + rotAngle=2 * pi / 3, + polarAngle=acos(1 / sqrt(3)), + aziAngle=5 * pi / 4, + label="AC_21") + elif cliffNum == 22: + #-X+Y+Z 120 + return arb_axis_drag(qubit, + nutFreq, + rotAngle=2 * pi / 3, + polarAngle=acos(1 / sqrt(3)), + aziAngle=3 * pi / 4, + label="AC_22") + elif cliffNum == 23: + #-X+Y+Z -120 (equivalent to X-Y-Z 120) + return arb_axis_drag(qubit, + nutFreq, + rotAngle=2 * pi / 3, + polarAngle=pi - acos(1 / sqrt(3)), + aziAngle=-pi / 4, + label="AC_23") + else: + raise ValueError('Clifford number must be between 0 and 23') + +def get_DiAC_phases(cliffNum): + """ + + Returns the phases (in multiples of pi) of the three Z gates dressing the two X90 + pulses comprising the DiAC pulse correspoding to cliffNum + e.g., get_DiAC_phases(1) returns a=0, b=1, c=1, in + Ztheta(a) + X90 + Ztheta(b) + X90 + Ztheta(c) = Id + """ + DiAC_table = [ + [0, 1, 1], + [0.5, -0.5, 0.5], + [0, 0, 0], + [0.5, 0.5, 0.5], + [0, -0.5, 1], + [0, 0, 1], + [0, 0.5, 1], + [0, 1, -0.5], + [0, 1, 0], + [0, 1, 0.5], + [0, 0, 0.5], + [0, 0, -0.5], + [1, -0.5, 1], + [1, 0.5, 1], + [0.5, -0.5, -0.5], + [0.5, 0.5, -0.5], + [0.5, -0.5, 1], + [1, -0.5, -0.5], + [0, 0.5, -0.5], + [-0.5, -0.5, 1], + [1, 0.5, -0.5], + [0.5, 0.5, 1], + [0, -0.5, -0.5], + [-0.5, 0.5, 1]] + return DiAC_table[cliffNum] + +def DiAC(qubit, cliffNum, compiled = True): + """ + + The set of 24 Diatomic Clifford single qubit pulses. Each pulse is decomposed + as Ztheta(a) + X90 + Ztheta(b) + X90 + Ztheta(c) if compiled = False, + uses also Y90, Y90m and shorter sequences if compiled = True + + Parameters + ---------- + qubit : logical channel to implement sequence (LogicalChannel) + cliffNum : the zero-indexed Clifford number + + Returns + ------- + pulse object + """ + #Now a big else if chain for to get the specific Clifford + if not compiled: + DiAC_phases = get_DiAC_phases(cliffNum) + return Ztheta(qubit, angle = DiAC_phases[0]*np.pi) + X90(qubit) + Ztheta(qubit, angle = DiAC_phases[1]*np.pi) + \ + X90(qubit) + Ztheta(qubit, angle = DiAC_phases[2]*np.pi) + else: + if cliffNum == 0: + #Identity gate + return Id(qubit, length=0) + elif cliffNum == 1: + #X90 + return X90(qubit) + elif cliffNum == 2: + #X180 + return X90(qubit)+X90(qubit) + elif cliffNum == 3: + #X90m + return X90m(qubit) + elif cliffNum == 4: + #Y90 + return Y90(qubit) + elif cliffNum == 5: + #Y180 + return Y90(qubit) + Y90(qubit) + elif cliffNum == 6: + #Y90m + return Y90m(qubit) + elif cliffNum == 7: + #Z90 + return Z90(qubit) + elif cliffNum == 8: + #Z180 + return Z(qubit) + elif cliffNum == 9: + #Z90m + return Z90m(qubit) + elif cliffNum == 10: + #X+Y 180 + return Y90(qubit) + Y90(qubit) + Z90m(qubit) + elif cliffNum == 11: + #X-Y 180 + return Y90(qubit) + Y90(qubit) + Z90(qubit) + elif cliffNum == 12: + #X+Z 180(Hadamard) + return Z(qubit) + Y90(qubit) + elif cliffNum == 13: + #X-Z 180 + return Z(qubit) + Y90m(qubit) + elif cliffNum == 14: + #Y+Z 180 + return Z90(qubit) + Y90(qubit) + Z90(qubit) + elif cliffNum == 15: + #Y-Z 180 + return Z90(qubit) + Y90m(qubit) + Z90(qubit) + elif cliffNum == 16: + #X+Y+Z -120 (equivalent to -X-Y-Z 120) + return Z90(qubit) + Y90(qubit) + elif cliffNum == 17: + #X+Y+Z 120 + return Z(qubit) + Y90(qubit) + Z90(qubit) + elif cliffNum == 18: + #X-Y+Z 120 (equivalent to -X+Y-Z 120) + return Y90m(qubit) + Z90(qubit) + elif cliffNum == 19: + #X-Y+Z -120 + return Z90m(qubit) + Y90(qubit) + elif cliffNum == 20: + #X+Y-Z -120 (equivalent to -X-Y+Z 120) + return Z(qubit) + Y90m(qubit) + Z90(qubit) + elif cliffNum == 21: + #X+Y-Z 120 + return Z90(qubit) + Y90m(qubit) + elif cliffNum == 22: + #-X+Y+Z -120 (equivalent to X-Y-Z 120) + return Y90(qubit) + Z90(qubit) + elif cliffNum == 23: + #-X+Y+Z 120 + return Z90m(qubit) + Y90m(qubit) + else: + raise ValueError('Clifford number must be between 0 and 23') + +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 + +### +### Two qubit Cliffords +### + +clifford_map = {} +clifford_map['STD'] = StdClifford +clifford_map['DIAC'] = DiAC +clifford_map['AC'] = AC +clifford_map['XYX'] = XYXClifford + +def Cx2(c1, c2, q1, q2, kind='std'): + """ + Helper function to create pulse block for a pair of single-qubit Cliffords + """ + + clifford_fun = clifford_map[kind.upper()] + seq1 = clifford_fun(q1, c1) + seq2 = clifford_fun(q2, c2) + + #Create the pulse block + return seq1 * seq2 def entangling_seq(gate, q1, q2): """ - Helper function to create the entangling gate sequence - """ + Helper function to create the entangling gate sequence + """ if gate == "CNOT": return ZX90_CR(q2, q1) elif gate == "iSWAP": @@ -50,39 +392,15 @@ def entangling_seq(gate, q1, q2): return [ZX90_CR(q2, q1), Y90m(q1) * Y90m(q2), ZX90_CR( q2, q1), (X90(q1) + Y90m(q1)) * X90(q2), ZX90_CR(q2, q1)] +def TwoQubitClifford(q1, q2, cliffNum, kind='std'): -def entangling_mat(gate): - """ - Helper function to create the entangling gate matrix - """ - echoCR = expm(1j * pi / 4 * np.kron(pX, pZ)) - if gate == "CNOT": - return echoCR - elif gate == "iSWAP": - return reduce(lambda x, y: np.dot(y, x), - [echoCR, np.kron(C1[6], C1[6]), echoCR]) - elif gate == "SWAP": - return reduce(lambda x, y: np.dot(y, x), - [echoCR, np.kron(C1[6], C1[6]), echoCR, np.kron( - np.dot(C1[6], C1[1]), C1[1]), echoCR]) - else: - raise ValueError("Entangling gate must be one of: CNOT, iSWAP, SWAP.") - + if kind.upper() not in clifford_map.keys(): + raise ValueError(f"Unknown clifford type: must be one of {clifford.map.keys()}.") -def clifford_seq(c, q1, q2=None): - """ - Return a sequence of pulses that implements a clifford C - """ - #If qubit2 not defined assume 1 qubit - if not q2: - genSeq = generatorSeqs[choice(C1Seqs[c])] - return [generator_pulse(g)(q1) for g in genSeq] - else: - #Look up the sequence for the integer - c = C2Seqs[c] - seq = [Cx2(c[0][0], c[0][1], q1, q2)] - if c[1]: - seq += entangling_seq(c[1], q1, q2) - if c[2]: - seq += [Cx2(c[2][0], c[2][1], q1, q2)] - return seq + c = C2Seqs[cliffNum] + seq = [Cx2(c[0][0], c[0][1], q1, q2, kind=kind)] + if c[1]: + seq += entangling_seq(c[1], q1, q2) + if c[2]: + seq += [Cx2(c[2][0], c[2][1], q1, q2, kind=kind)] + return seq \ No newline at end of file diff --git a/QGL/tools/clifford_tools.py b/QGL/tools/clifford_tools.py index 566372de..244513ba 100644 --- a/QGL/tools/clifford_tools.py +++ b/QGL/tools/clifford_tools.py @@ -180,6 +180,24 @@ def clifford_mat(c, numQubits): return mat +def entangling_mat(gate): + """ + Helper function to create the entangling gate matrix + """ + echoCR = expm(1j * pi / 4 * np.kron(pX, pZ)) + if gate == "CNOT": + return echoCR + elif gate == "iSWAP": + return reduce(lambda x, y: np.dot(y, x), + [echoCR, np.kron(C1[6], C1[6]), echoCR]) + elif gate == "SWAP": + return reduce(lambda x, y: np.dot(y, x), + [echoCR, np.kron(C1[6], C1[6]), echoCR, np.kron( + np.dot(C1[6], C1[1]), C1[1]), echoCR]) + else: + raise ValueError("Entangling gate must be one of: CNOT, iSWAP, SWAP.") + + def inverse_clifford(cMat): """Return the inverse clifford index.""" dim = cMat.shape[0] diff --git a/QGL/tools/euler_angles.py b/QGL/tools/euler_angles.py index c9b197db..4ee03f12 100644 --- a/QGL/tools/euler_angles.py +++ b/QGL/tools/euler_angles.py @@ -20,9 +20,8 @@ import numpy as np from scipy.linalg import expm -from .Cliffords import C1 -from .PulsePrimitives import * -from .tools.matrix_tools import * +from .clifford_tools import C1 +from .matrix_tools import * def xyx_unitary(α, β, γ): """Unitary decomposed as Rx, Ry, Rx rotations. From 561ffe2298642090a236253dec9c40d520935e80 Mon Sep 17 00:00:00 2001 From: Guilhem Ribeill Date: Wed, 15 Apr 2020 18:05:19 -0400 Subject: [PATCH 10/13] WIP: Update RB to use the new functions --- QGL/BasicSequences/RB.py | 102 ++++---------- QGL/PulsePrimitives.py | 289 --------------------------------------- QGL/__init__.py | 1 + 3 files changed, 31 insertions(+), 361 deletions(-) diff --git a/QGL/BasicSequences/RB.py b/QGL/BasicSequences/RB.py index f3393f3f..1c5fd707 100644 --- a/QGL/BasicSequences/RB.py +++ b/QGL/BasicSequences/RB.py @@ -1,9 +1,10 @@ from ..PulsePrimitives import * +from ..Cliffords import * from ..Compiler import compile_to_hardware from ..PulseSequencePlotter import plot_pulse_files -#from ..tools.clifford_tools import clifford_seq, clifford_mat, inverse_clifford -#from ..Euler import XYXClifford +from ..tools.clifford_tools import clifford_mat, inverse_clifford from .helpers import create_cal_seqs, cal_descriptor +from ..config import logger import os from csv import reader @@ -88,7 +89,7 @@ def create_RB_seqs(numQubits, return seqs -def SingleQubitRB(qubit, seqs, purity=False, showPlot=False, add_cals=True): +def SingleQubitRB(qubit, seqs, cliff_type='std', purity=False, showPlot=False, add_cals=True): """ Single qubit randomized benchmarking using 90 and 180 generators. @@ -120,12 +121,16 @@ def SingleQubitRB(qubit, seqs, purity=False, showPlot=False, add_cals=True): '/path/to/exp/exp-meta.json' """ + if cliff_type.upper() not in clifford_map.keys(): + raise ValueError(f"Unknown clifford type: must be one of {clifford.map.keys()}.") + + clifford = clifford_map[cliff_type] + 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(reduce(operator.add, [clifford_seq(c, qubit) - for c in seq])) + seqsBis.append([clifford(qubit,c) for c in seq]) #append tomography pulse to measure purity seqsBis[-1].append(op[ct]) #append measurement @@ -149,7 +154,7 @@ def SingleQubitRB(qubit, seqs, purity=False, showPlot=False, add_cals=True): plot_pulse_files(metafile) return metafile -def SingleQubitLeakageRB(qubit, seqs, pi2args, showPlot=False): +def SingleQubitLeakageRB(qubit, seqs, pi2args, cliff_type='std', showPlot=False): """ Single qubit randomized benchmarking using 90 and 180 generators to measure leakage outside the qubit subspace. @@ -183,9 +188,14 @@ def SingleQubitLeakageRB(qubit, seqs, pi2args, showPlot=False): '/path/to/exp/exp-meta.json' """ + if cliff_type.upper() not in clifford_map.keys(): + raise ValueError(f"Unknown clifford type: must be one of {clifford.map.keys()}.") + + clifford = clifford_map[cliff_type] + seqsBis = [] for seq in seqs: - combined_seq = reduce(operator.add, [clifford_seq(c, qubit) for c in seq]) + combined_seq = [clifford(qubit, c) for c in seq] # Append sequence with tomography ids and measurement seqsBis.append(combined_seq + [Id(qubit), Id(qubit), MEAS(qubit)]) @@ -220,7 +230,7 @@ def SingleQubitLeakageRB(qubit, seqs, pi2args, showPlot=False): -def TwoQubitRB(q1, q2, seqs, showPlot=False, suffix="", add_cals=True): +def TwoQubitRB(q1, q2, seqs, cliff_type='std', showPlot=False, suffix="", add_cals=True): """ Two qubit randomized benchmarking using 90 and 180 single qubit generators and ZX90. @@ -254,9 +264,10 @@ def TwoQubitRB(q1, q2, seqs, showPlot=False, suffix="", add_cals=True): >>> mf '/path/to/exp/exp-meta.json' """ + seqsBis = [] for seq in seqs: - seqsBis.append(reduce(operator.add, [clifford_seq(c, q2, q1) + seqsBis.append(reduce(operator.add, [TwoQubitClifford(q2, q1, kind=cliff_type) for c in seq])) #Add the measurement to all sequences @@ -281,7 +292,7 @@ def TwoQubitRB(q1, q2, seqs, showPlot=False, suffix="", add_cals=True): plot_pulse_files(metafile) return metafile -def TwoQubitLeakageRB(q1, q2, meas_qubit, seqs, pi2args, showPlot=False): +def TwoQubitLeakageRB(q1, q2, meas_qubit, seqs, pi2args, cliff_type='std', showPlot=False): """ Two qubit randomized benchmarking using 90 and 180 single qubit generators and ZX90 to measure leakage outside the qubit subspace. See https:// @@ -320,7 +331,7 @@ def TwoQubitLeakageRB(q1, q2, meas_qubit, seqs, pi2args, showPlot=False): """ seqsBis = [] for seq in seqs: - combined_seq = reduce(operator.add, [clifford_seq(c, q2, q1) for c in seq]) + combined_seq = reduce(operator.add, [TwoQubitClifford(q2, q1, c, kind=cliff_type) for c in seq]) # Append sequence with tomography ids and measurement seqsBis.append(combined_seq + [Id(meas_qubit), Id(meas_qubit), MEAS(meas_qubit)]) @@ -386,6 +397,10 @@ def SingleQubitRB_AC(qubit, seqs, purity=False, showPlot=False, add_cals=True): >>> mf '/path/to/exp/exp-meta.json' """ + + logger.warning("This function is deprecated and may be removed in a future release of QGL! " + + "Use `SingleQubitRB` with the `cliff_type` keyword argument instead.") + seqsBis = [] op = [Id(qubit, length=0), Y90m(qubit), X90(qubit)] for ct in range(3 if purity else 1): @@ -454,6 +469,10 @@ def SingleQubitRB_DiAC(qubit, >>> mf '/path/to/exp/exp-meta.json' """ + + logger.warning("This function is deprecated and may be removed in a future release of QGL! " + + "Use `SingleQubitRB` with the `cliff_type` keyword argument instead.") + seqsBis = [] op = [Id(qubit, length=0), Y90m(qubit), X90(qubit)] for ct in range(3 if purity else 1): @@ -482,67 +501,6 @@ 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/PulsePrimitives.py b/QGL/PulsePrimitives.py index ba908919..22adee16 100644 --- a/QGL/PulsePrimitives.py +++ b/QGL/PulsePrimitives.py @@ -337,295 +337,6 @@ def arb_axis_drag(qubit, return Pulse(kwargs["label"] if "label" in kwargs else "ArbAxis", qubit, params, 1.0, aziAngle, frameChange) - -def AC(qubit, cliffNum): - """ - - The set of 24 Atomic Clifford single qubit pulses. - - Parameters - ---------- - qubit : logical channel to implement sequence (LogicalChannel) - cliffNum : the zero-indexed Clifford number - - Returns - ------- - pulse object - """ - - #Figure out the approximate nutation frequency calibration from the X180 and the sampling_rate - Xp = X(qubit) - xpulse = Xp.amp * Xp.shape - nutFreq = 0.5 / (sum(xpulse) / qubit.phys_chan.sampling_rate) - - #Now a big else if chain for to get the specific Clifford - if cliffNum == 0: - #Identity gate - return Id(qubit, length=0) - elif cliffNum == 1: - #X90 - return X90(qubit) - elif cliffNum == 2: - #X180 - return X(qubit) - elif cliffNum == 3: - #X90m - return X90m(qubit) - elif cliffNum == 4: - #Y90 - return Y90(qubit) - elif cliffNum == 5: - #Y180 - return Y(qubit) - elif cliffNum == 6: - #Y90m - return Y90m(qubit) - elif cliffNum == 7: - #Z90 - return Z90(qubit) - elif cliffNum == 8: - #Z180 - return Z(qubit) - elif cliffNum == 9: - #Z90m - return Z90m(qubit) - elif cliffNum == 10: - #X+Y 180 - return U(qubit, phase=pi / 4, label="AC_10") - elif cliffNum == 11: - #X-Y 180 - return U(qubit, phase=-pi / 4, label="AC_11") - elif cliffNum == 12: - #X+Z 180(Hadamard) - return arb_axis_drag(qubit, - nutFreq, - rotAngle=pi, - polarAngle=pi / 4, - label="AC_12") - elif cliffNum == 13: - #X-Z 180 - return arb_axis_drag(qubit, - nutFreq, - rotAngle=pi, - polarAngle=pi / 4, - aziAngle=pi, - label="AC_13") - elif cliffNum == 14: - #Y+Z 180 - return arb_axis_drag(qubit, - nutFreq, - rotAngle=pi, - polarAngle=pi / 4, - aziAngle=pi / 2, - label="AC_14") - elif cliffNum == 15: - #Y-Z 180 - return arb_axis_drag(qubit, - nutFreq, - rotAngle=pi, - polarAngle=pi / 4, - aziAngle=-pi / 2, - label="AC_15") - elif cliffNum == 16: - #X+Y+Z 120 - return arb_axis_drag(qubit, - nutFreq, - rotAngle=2 * pi / 3, - polarAngle=acos(1 / sqrt(3)), - aziAngle=pi / 4, - label="AC_16") - elif cliffNum == 17: - #X+Y+Z -120 (equivalent to -X-Y-Z 120) - return arb_axis_drag(qubit, - nutFreq, - rotAngle=2 * pi / 3, - polarAngle=pi - acos(1 / sqrt(3)), - aziAngle=5 * pi / 4, - label="AC_17") - elif cliffNum == 18: - #X-Y+Z 120 - return arb_axis_drag(qubit, - nutFreq, - rotAngle=2 * pi / 3, - polarAngle=acos(1 / sqrt(3)), - aziAngle=-pi / 4, - label="AC_18") - elif cliffNum == 19: - #X-Y+Z 120 (equivalent to -X+Y-Z -120) - return arb_axis_drag(qubit, - nutFreq, - rotAngle=2 * pi / 3, - polarAngle=pi - acos(1 / sqrt(3)), - aziAngle=3 * pi / 4, - label="AC_19") - elif cliffNum == 20: - #X+Y-Z 120 - return arb_axis_drag(qubit, - nutFreq, - rotAngle=2 * pi / 3, - polarAngle=pi - acos(1 / sqrt(3)), - aziAngle=pi / 4, - label="AC_20") - elif cliffNum == 21: - #X+Y-Z -120 (equivalent to -X-Y+Z 120) - return arb_axis_drag(qubit, - nutFreq, - rotAngle=2 * pi / 3, - polarAngle=acos(1 / sqrt(3)), - aziAngle=5 * pi / 4, - label="AC_21") - elif cliffNum == 22: - #-X+Y+Z 120 - return arb_axis_drag(qubit, - nutFreq, - rotAngle=2 * pi / 3, - polarAngle=acos(1 / sqrt(3)), - aziAngle=3 * pi / 4, - label="AC_22") - elif cliffNum == 23: - #-X+Y+Z -120 (equivalent to X-Y-Z 120) - return arb_axis_drag(qubit, - nutFreq, - rotAngle=2 * pi / 3, - polarAngle=pi - acos(1 / sqrt(3)), - aziAngle=-pi / 4, - label="AC_23") - else: - raise ValueError('Clifford number must be between 0 and 23') - - -def get_DiAC_phases(cliffNum): - """ - - Returns the phases (in multiples of pi) of the three Z gates dressing the two X90 - pulses comprising the DiAC pulse correspoding to cliffNum - e.g., get_DiAC_phases(1) returns a=0, b=1, c=1, in - Ztheta(a) + X90 + Ztheta(b) + X90 + Ztheta(c) = Id - """ - DiAC_table = [ - [0, 1, 1], - [0.5, -0.5, 0.5], - [0, 0, 0], - [0.5, 0.5, 0.5], - [0, -0.5, 1], - [0, 0, 1], - [0, 0.5, 1], - [0, 1, -0.5], - [0, 1, 0], - [0, 1, 0.5], - [0, 0, 0.5], - [0, 0, -0.5], - [1, -0.5, 1], - [1, 0.5, 1], - [0.5, -0.5, -0.5], - [0.5, 0.5, -0.5], - [0.5, -0.5, 1], - [1, -0.5, -0.5], - [0, 0.5, -0.5], - [-0.5, -0.5, 1], - [1, 0.5, -0.5], - [0.5, 0.5, 1], - [0, -0.5, -0.5], - [-0.5, 0.5, 1]] - return DiAC_table[cliffNum] - -def DiAC(qubit, cliffNum, compiled = True): - """ - - The set of 24 Diatomic Clifford single qubit pulses. Each pulse is decomposed - as Ztheta(a) + X90 + Ztheta(b) + X90 + Ztheta(c) if compiled = False, - uses also Y90, Y90m and shorter sequences if compiled = True - - Parameters - ---------- - qubit : logical channel to implement sequence (LogicalChannel) - cliffNum : the zero-indexed Clifford number - - Returns - ------- - pulse object - """ - #Now a big else if chain for to get the specific Clifford - if not compiled: - DiAC_phases = get_DiAC_phases(cliffNum) - return Ztheta(qubit, angle = DiAC_phases[0]*np.pi) + X90(qubit) + Ztheta(qubit, angle = DiAC_phases[1]*np.pi) + \ - X90(qubit) + Ztheta(qubit, angle = DiAC_phases[2]*np.pi) - else: - if cliffNum == 0: - #Identity gate - return Id(qubit, length=0) - elif cliffNum == 1: - #X90 - return X90(qubit) - elif cliffNum == 2: - #X180 - return X90(qubit)+X90(qubit) - elif cliffNum == 3: - #X90m - return X90m(qubit) - elif cliffNum == 4: - #Y90 - return Y90(qubit) - elif cliffNum == 5: - #Y180 - return Y90(qubit) + Y90(qubit) - elif cliffNum == 6: - #Y90m - return Y90m(qubit) - elif cliffNum == 7: - #Z90 - return Z90(qubit) - elif cliffNum == 8: - #Z180 - return Z(qubit) - elif cliffNum == 9: - #Z90m - return Z90m(qubit) - elif cliffNum == 10: - #X+Y 180 - return Y90(qubit) + Y90(qubit) + Z90m(qubit) - elif cliffNum == 11: - #X-Y 180 - return Y90(qubit) + Y90(qubit) + Z90(qubit) - elif cliffNum == 12: - #X+Z 180(Hadamard) - return Z(qubit) + Y90(qubit) - elif cliffNum == 13: - #X-Z 180 - return Z(qubit) + Y90m(qubit) - elif cliffNum == 14: - #Y+Z 180 - return Z90(qubit) + Y90(qubit) + Z90(qubit) - elif cliffNum == 15: - #Y-Z 180 - return Z90(qubit) + Y90m(qubit) + Z90(qubit) - elif cliffNum == 16: - #X+Y+Z -120 (equivalent to -X-Y-Z 120) - return Z90(qubit) + Y90(qubit) - elif cliffNum == 17: - #X+Y+Z 120 - return Z(qubit) + Y90(qubit) + Z90(qubit) - elif cliffNum == 18: - #X-Y+Z 120 (equivalent to -X+Y-Z 120) - return Y90m(qubit) + Z90(qubit) - elif cliffNum == 19: - #X-Y+Z -120 - return Z90m(qubit) + Y90(qubit) - elif cliffNum == 20: - #X+Y-Z -120 (equivalent to -X-Y+Z 120) - return Z(qubit) + Y90m(qubit) + Z90(qubit) - elif cliffNum == 21: - #X+Y-Z 120 - return Z90(qubit) + Y90m(qubit) - elif cliffNum == 22: - #-X+Y+Z -120 (equivalent to X-Y-Z 120) - return Y90(qubit) + Z90(qubit) - elif cliffNum == 23: - #-X+Y+Z 120 - return Z90m(qubit) + Y90m(qubit) - else: - raise ValueError('Clifford number must be between 0 and 23') - ## two-qubit primitivies # helper used by echoCR diff --git a/QGL/__init__.py b/QGL/__init__.py index 82defa84..2df5915b 100644 --- a/QGL/__init__.py +++ b/QGL/__init__.py @@ -9,6 +9,7 @@ from .Plotting import build_waveforms #, show, , plot_waveforms from .PulseSequencePlotter import plot_pulse_files from .Tomography import state_tomo, process_tomo +from .Cliffords import StdClifford, DiAC, AC, XYXClifford, TwoQubitClifford from .Scheduler import schedule from .TdmInstructions import MajorityMask, MajorityVote, WriteAddr, Invalidate, Decode, DecodeSetRounds From 5d77040e16f57b0299463d2c255fec30850aa1a0 Mon Sep 17 00:00:00 2001 From: Guilhem Ribeill Date: Wed, 15 Apr 2020 18:29:49 -0400 Subject: [PATCH 11/13] Update tests. --- tests/test_clifford.py | 39 +++++++++++++++++++++++++++++++++++++++ tests/test_euler.py | 4 ++-- 2 files changed, 41 insertions(+), 2 deletions(-) create mode 100644 tests/test_clifford.py diff --git a/tests/test_clifford.py b/tests/test_clifford.py new file mode 100644 index 00000000..ecd01649 --- /dev/null +++ b/tests/test_clifford.py @@ -0,0 +1,39 @@ +import unittest + +from QGL import * +from QGL.tools.matrix_tools import * +from QGL.tools.clifford_tools import * +import QGL.config +try: + from helpers import setup_test_lib +except: + from .helpers import setup_test_lib + +class EulerDecompositions(unittest.TestCase): + + def setUp(self): + self.N1 = 24 #number of single qubit Cliffords + self.N2 = 11520 #number of two qubit Cliffords + self.N_test_2 = 30 #number of two qubit Cliffords to test + + def test_n_cliffords(self): + assert len(C1Seqs) == 24 + assert len(C2Seqs) == 11520 + + def test_multiply(self): + for j, k in product(range(self.N1), range(self.N1)): + m = C1[clifford_multiply(j, k)] + mtemp = (C1[k]@C1[j]).transpose().conj() + assert np.isclose(np.abs((mtemp@m).trace()), 2.0) + + def test_inverse(self): + for j in range(self.N1): + inv = C1[inverse_clifford(C1[j])] + assert is_close(inv@C1[j], pI) + for j in np.random.choice(range(11520), self.N_test_2): + C = clifford_mat(j, 2) + Ci = clifford_mat(inverse_clifford(C), 2) + assert np.isclose(np.abs((Ci@C).trace()), 4.0) + +if __name__ == "__main__": + unittest.main() diff --git a/tests/test_euler.py b/tests/test_euler.py index a19eb4c0..071f01eb 100644 --- a/tests/test_euler.py +++ b/tests/test_euler.py @@ -1,9 +1,9 @@ import unittest from QGL import * -from QGL.Euler import * +from QGL.tools.euler_angles import * from QGL.tools.matrix_tools import * -from QGL.Cliffords import C1 +from QGL.tools.clifford_tools import C1 import QGL.config try: from helpers import setup_test_lib From 14f127a33255d9ffc4fa1d124fe8b2c392221b47 Mon Sep 17 00:00:00 2001 From: Guilhem Ribeill Date: Wed, 15 Apr 2020 18:40:10 -0400 Subject: [PATCH 12/13] Fix keyerror in clifford types --- QGL/BasicSequences/RB.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/QGL/BasicSequences/RB.py b/QGL/BasicSequences/RB.py index 1c5fd707..af8fcea6 100644 --- a/QGL/BasicSequences/RB.py +++ b/QGL/BasicSequences/RB.py @@ -124,7 +124,7 @@ def SingleQubitRB(qubit, seqs, cliff_type='std', purity=False, showPlot=False, a if cliff_type.upper() not in clifford_map.keys(): raise ValueError(f"Unknown clifford type: must be one of {clifford.map.keys()}.") - clifford = clifford_map[cliff_type] + clifford = clifford_map[cliff_type.upper()] seqsBis = [] op = [Id(qubit, length=0), Y90m(qubit), X90(qubit)] @@ -191,7 +191,7 @@ def SingleQubitLeakageRB(qubit, seqs, pi2args, cliff_type='std', showPlot=False) if cliff_type.upper() not in clifford_map.keys(): raise ValueError(f"Unknown clifford type: must be one of {clifford.map.keys()}.") - clifford = clifford_map[cliff_type] + clifford = clifford_map[cliff_type.upper()] seqsBis = [] for seq in seqs: From 750a31c1357e38b74895f5f978aa31732335178c Mon Sep 17 00:00:00 2001 From: Guilhem Ribeill Date: Thu, 16 Apr 2020 15:05:22 -0400 Subject: [PATCH 13/13] Fixes for two-qubit RB --- QGL/BasicSequences/RB.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/QGL/BasicSequences/RB.py b/QGL/BasicSequences/RB.py index af8fcea6..5c82a537 100644 --- a/QGL/BasicSequences/RB.py +++ b/QGL/BasicSequences/RB.py @@ -230,7 +230,7 @@ def SingleQubitLeakageRB(qubit, seqs, pi2args, cliff_type='std', showPlot=False) -def TwoQubitRB(q1, q2, seqs, cliff_type='std', showPlot=False, suffix="", add_cals=True): +def TwoQubitRB(q1, q2, seqs, meas_qubits='all', cliff_type='std', showPlot=False, suffix="", add_cals=True): """ Two qubit randomized benchmarking using 90 and 180 single qubit generators and ZX90. @@ -267,12 +267,15 @@ def TwoQubitRB(q1, q2, seqs, cliff_type='std', showPlot=False, suffix="", add_ca seqsBis = [] for seq in seqs: - seqsBis.append(reduce(operator.add, [TwoQubitClifford(q2, q1, kind=cliff_type) + seqsBis.append(reduce(operator.add, [TwoQubitClifford(q2, q1, c, kind=cliff_type) for c in seq])) #Add the measurement to all sequences for seq in seqsBis: - seq.append(MEAS(q1) * MEAS(q2)) + if meas_qubits.upper() == "ALL": + seq.append(MEAS(q1) * MEAS(q2)) + else: + seq.append(reduce(operator.mul, [MEAS(q) for q in meas_qubits])) axis_descriptor = [{ 'name': 'length',