From 04f1a42a2eacca6db7cafe5385b43238895d3017 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Thu, 15 Feb 2024 17:11:32 -0500 Subject: [PATCH 01/51] bugfix for my attemped efficiency improvement --- pygsti/tools/matrixtools.py | 40 +++++++++++++++++++++++++++++++++++++ 1 file changed, 40 insertions(+) diff --git a/pygsti/tools/matrixtools.py b/pygsti/tools/matrixtools.py index 94940c45c..bffa3fb89 100644 --- a/pygsti/tools/matrixtools.py +++ b/pygsti/tools/matrixtools.py @@ -138,6 +138,46 @@ def is_valid_density_mx(mx, tol=1e-9): return abs(_np.trace(mx) - 1.0) < tol and is_pos_def(mx, tol) +def frobeniusnorm(ar): + """ + Compute the frobenius norm of an array (or matrix), + + sqrt( sum( each_element_of_a^2 ) ) + + Parameters + ---------- + ar : numpy array + What to compute the frobenius norm of. Note that ar can be any shape + or number of dimenions. + + Returns + ------- + float or complex + depending on the element type of ar. + """ + return _np.linalg.norm(ar.ravel()) + + +def frobeniusnorm_squared(ar): + """ + Compute the squared frobenius norm of an array (or matrix), + + sum( each_element_of_a^2 ) ) + + Parameters + ---------- + ar : numpy array + What to compute the squared frobenius norm of. Note that ar can be any + shape or number of dimenions. + + Returns + ------- + float or complex + depending on the element type of ar. + """ + return _np.linalg.norm(ar.ravel())**2 + + def nullspace(m, tol=1e-7): """ Compute the nullspace of a matrix. From 5252e41269d1682b1e1aebafcff1b7c214b738db Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Thu, 15 Feb 2024 17:55:48 -0500 Subject: [PATCH 02/51] extend basistools to gracefully accomodate objects that derive from LinearOperator (not just plain numpy ndarrays) --- pygsti/modelmembers/operations/linearop.py | 10 +++++++++- pygsti/tools/basistools.py | 10 ++++++++++ 2 files changed, 19 insertions(+), 1 deletion(-) diff --git a/pygsti/modelmembers/operations/linearop.py b/pygsti/modelmembers/operations/linearop.py index c86352a9f..81a871c6f 100644 --- a/pygsti/modelmembers/operations/linearop.py +++ b/pygsti/modelmembers/operations/linearop.py @@ -44,7 +44,7 @@ class LinearOperator(_modelmember.ModelMember): """ - Base class for all operation representations + Base class for all *square* operation representations Parameters ---------- @@ -92,6 +92,14 @@ def hilbert_schmidt_size(self): """ return (self.dim)**2 + @property + def shape(self): + # Provide this function to mimic numpy array semantics. + # + # We can't rely on self._rep.shape since superclasses + # are given broad freedom to define semantics of self._rep. + return (self.dim, self.dim) + def set_dense(self, m): """ Set the dense-matrix value of this operation. diff --git a/pygsti/tools/basistools.py b/pygsti/tools/basistools.py index 2ec896cd5..a380f2dfd 100644 --- a/pygsti/tools/basistools.py +++ b/pygsti/tools/basistools.py @@ -141,6 +141,13 @@ def change_basis(mx, from_basis, to_basis): The given operation matrix converted to the `to_basis` basis. Array size is the same as `mx`. """ + if isinstance(mx, _np.ndarray) and mx.ndim == 0: + # mx is probably a wrapper around a pyGSTi object. + mx = mx.item() + if hasattr(mx, 'to_dense'): + # Question: how do we know the returned representation is + # in fact in the "from_basis"? + mx = mx.to_dense() if len(mx.shape) not in (1, 2): raise ValueError("Invalid dimension of object - must be 1 or 2, i.e. a vector or matrix") @@ -272,6 +279,9 @@ def create_basis_for_matrix(mx, basis): ------- Basis """ + if isinstance(mx, _np.ndarray) and mx.ndim == 0: + # ^ mx is probably just holding a pyGSTi object + mx = mx.item() dim = mx.shape[0] if isinstance(basis, _basis.Basis): assert(basis.dim == dim), "Supplied Basis has wrong dimension!" From f531c8bbfb85e5c5e180db46a3e8fd0a97cb0c2d Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Wed, 21 Feb 2024 11:08:47 -0500 Subject: [PATCH 03/51] readability tweak and comments --- pygsti/tools/optools.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pygsti/tools/optools.py b/pygsti/tools/optools.py index 8d68a73ff..e64401b21 100644 --- a/pygsti/tools/optools.py +++ b/pygsti/tools/optools.py @@ -491,8 +491,8 @@ def entanglement_fidelity(a, b, mx_basis='pp', is_tp=None, is_unitary=None): #if the tp flag isn't set we'll calculate whether it is true here if is_tp is None: - def is_tp_fn(x): return _np.isclose(x[0, 0], 1.0) and all( - [_np.isclose(x[0, i], 0) for i in range(1,d2)]) + def is_tp_fn(x): + return _np.isclose(x[0, 0], 1.0) and _np.allclose(x[0, 1:d2], 0) is_tp= (is_tp_fn(a) and is_tp_fn(b)) From 2b34cae5244a6511406036459827936e93cc45ed Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Wed, 21 Feb 2024 19:19:22 -0500 Subject: [PATCH 04/51] update function in reportables to indicate API limitation. Bugfix in frobenius norm computation. Bugfix in test_optools.py. --- pygsti/report/reportables.py | 7 +++++-- test/unit/tools/test_optools.py | 15 ++++++++++++++- 2 files changed, 19 insertions(+), 3 deletions(-) diff --git a/pygsti/report/reportables.py b/pygsti/report/reportables.py index e4fa68146..1a0cf4152 100644 --- a/pygsti/report/reportables.py +++ b/pygsti/report/reportables.py @@ -893,13 +893,16 @@ def upper_bound_fidelity(gate, mx_basis): gate : numpy.ndarray the transfer-matrix specifying a gate's action. - mx_basis : Basis or {'pp', 'gm', 'std'} - the basis that `gate` is in. + mx_basis : Basis or string + Currently restricted to Pauli-product Returns ------- float """ + basis_str = mx_basis if isinstance(mx_basis, str) else mx_basis.name + if basis_str != 'pp': + raise NotImplementedError(f'Basis must be Pauli-Product, got {mx_basis}.') return _tools.fidelity_upper_bound(gate)[0] diff --git a/test/unit/tools/test_optools.py b/test/unit/tools/test_optools.py index a93c9bc5d..23c8c36fa 100644 --- a/test/unit/tools/test_optools.py +++ b/test/unit/tools/test_optools.py @@ -50,7 +50,7 @@ def test_unitary_to_pauligate(self): # U_2Q is 4x4 unitary matrix operating on isolated two-qubit space (CX(pi) rotation) op_2Q = ot.unitary_to_pauligate(U_2Q) - op_2Q_inv = ot.process_mx_to_unitary(bt.change_basis(op_2Q, 'pp', 'std')) + op_2Q_inv = ot.std_process_mx_to_unitary(bt.change_basis(op_2Q, 'pp', 'std')) self.assertArraysAlmostEqual(U_2Q, op_2Q_inv) def test_decompose_gate_matrix(self): @@ -379,6 +379,13 @@ def setUp(self): [-0.35432747-0.27939404j, -0.02266757+0.71502652j, -0.27452307+0.07511567j, 0.35432747+0.27939404j], [ 0.71538573+0.j, 0.2680266 +0.36300238j, 0.2680266 -0.36300238j, 0.28461427+0.j]]) + def test_frobenius_distance(self): + self.assertAlmostEqual(ot.frobeniusdist(self.A, self.A), 0.0) + self.assertAlmostEqual(ot.frobeniusdist(self.A, self.B), 0.6204836823) + + self.assertAlmostEqual(ot.frobeniusdist_squared(self.A, self.A), 0.0) + self.assertAlmostEqual(ot.frobeniusdist_squared(self.A, self.B), 0.385) + def test_jtrace_distance(self): val = ot.jtracedist(self.A_TP, self.A_TP, mx_basis="pp") self.assertAlmostEqual(val, 0.0) @@ -412,3 +419,9 @@ def test_fidelity_upper_bound(self): bad_superoperator = ot.unitary_to_superop(Q) upperBound, _ = ot.fidelity_upper_bound(bad_superoperator) self.assertAlmostEqual(upperBound, 0.75) + np.random.seed(0) + Q = np.linalg.qr(np.random.randn(4,4) + 1j*np.random.randn(4,4))[0] + Q[:, 0] = 0.0 # zero out the first column + bad_superoperator = ot.unitary_to_superop(Q) + upperBound, _ = ot.fidelity_upper_bound(bad_superoperator) + self.assertAlmostEqual(upperBound, 0.75) From 2a8959ebd03d238362f17467fe35f97e6e9fe293 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Wed, 21 Feb 2024 19:20:18 -0500 Subject: [PATCH 05/51] update implementation of state fidelity --- pygsti/tools/optools.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pygsti/tools/optools.py b/pygsti/tools/optools.py index e64401b21..36cce2a57 100644 --- a/pygsti/tools/optools.py +++ b/pygsti/tools/optools.py @@ -570,7 +570,7 @@ def average_gate_fidelity(a, b, mx_basis='pp', is_tp=None, is_unitary=None): d = int(round(_np.sqrt(a.shape[0]))) PF = entanglement_fidelity(a, b, mx_basis, is_tp, is_unitary) AGF = (d * PF + 1) / (1 + d) - return float(AGF) + return AGF def average_gate_infidelity(a, b, mx_basis='pp', is_tp=None, is_unitary=None): @@ -660,7 +660,7 @@ def entanglement_infidelity(a, b, mx_basis='pp', is_tp=None, is_unitary=None): EI : float The EI of a to b. """ - return 1 - float(entanglement_fidelity(a, b, mx_basis, is_tp, is_unitary)) + return 1 - entanglement_fidelity(a, b, mx_basis, is_tp, is_unitary) def gateset_infidelity(model, target_model, itype='EI', From 0bcc11690a1501ff6af384befd9432a538068bf3 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Wed, 28 Feb 2024 13:53:34 -0500 Subject: [PATCH 06/51] comments. About to remove. --- pygsti/tools/optools.py | 102 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 102 insertions(+) diff --git a/pygsti/tools/optools.py b/pygsti/tools/optools.py index 36cce2a57..078750948 100644 --- a/pygsti/tools/optools.py +++ b/pygsti/tools/optools.py @@ -518,6 +518,108 @@ def is_tp_fn(x): return fidelity(JA, JB) +def simple_entanglement_fidelity(a_, b_, mx_basis, n_leak=0): + dim = int(_np.sqrt(a_.shape[0])) + assert a_.shape == (dim**2, dim**2) + assert b_.shape == (dim**2, dim**2) + + # set psi = (sum_i |ii>)/sqrt(dim), + # where |ii> is the tensor product of the i^th basis vector in C^dim with itself. + I = _np.eye(dim, dtype=_np.complex128) + summands = [] + for i in range(dim - n_leak): + s = _np.outer(I[i], I[i]) + # s = |i> \tensor \tensor |i> === |ii> + # ... unclear if that's certain to be the case. + summands.append(s_superket) + psi = _np.sum(summands, axis=0) / _np.sqrt(dim - n_leak) + proj_psi = _np.outer(psi, psi) + # ^ That is absolutely positively correct! + proj_psi_superket = _bt.stdmx_to_ppvec(proj_psi).ravel() + # ^ Robin thinks things are breaking here. + # The ordering of basis vectors isn't respecting tensor + # product structure. That is, proj_psi_superket + # doesn't have the same vectorized tensor product + # structure convention as kron (or kron with the order changed). + # + # TODO: figure out how to order basis elements for stdvec + # so that tensor product structure in stdmx is preserved + # in stdvec. + # + # X = [X11, X12] , size(X11) == size(X22) + # [X21, X22] + # + # basis(space(X)) = [basis(space(X11)), basis(space(X12)),..., basis(space(X22))] + # ^ Concatenation. + # ^ The vectorization of the matrix commutes with tensor product. + # + # ^ The tensor product of the bases for matrix units of X11, X12. + # + # If we just convert X from stdmx to stdvec, then we're doing something like + # basis(space(X)) = [basis(space(row(X, 1))), basis(space(row(X, 2))), ..., basis(space(row(X, n)))]. + # ^ Vectorization of the matrix does not commute with tensor product. + # + # + """ + We are currently taking two 2-dimensional vector spaces (U and V) and: + Tensoring them together to get a 4-d space W + Constructing a basis of matrix units for the 16-d space B(W) + Instead we should take U and V and: + Construct a basis B1 of matrix units for the 4-d space B(U) + Construct a basis B2 of matrix units for the 4-d space B(V) + Tensor B1 \otimes B2 to get a 16-d space. + + ^ Super easy to handle with powers of 2 and Pauli-Products. + ... NOT super easy to handle in ... literally anything else. + (but also not crazy hard. Just be careful.) + """ + + a = _mt.change_basis(a_, mx_basis, 'pp') + b = _mt.change_basis(b_, mx_basis, 'pp') + idle_gate = _np.eye(dim**2, dtype=_np.complex128) + a_tensor_idle = _np.kron(a, idle_gate) + b_tensor_idle = _np.kron(b, idle_gate) + temp1 = a_tensor_idle @ proj_psi_superket + temp2 = b_tensor_idle @ proj_psi_superket + ent_fid = _np.real(temp1.conj() @ temp2) + return ent_fid + + +def simple_entanglement_fidelity2(a_, b_, mx_basis, n_leak=0): + import numpy.linalg as la + dim = int(_np.sqrt(a_.shape[0])) + assert a_.shape == (dim**2, dim**2) + assert b_.shape == (dim**2, dim**2) + a = _mt.change_basis(a_, mx_basis, 'std') + b = _mt.change_basis(b_, mx_basis, 'std') + # set psi = (sum_i |ii>)/sqrt(dim), + # where |ii> is the tensor product of the i^th basis vector in C^dim with itself. + I = _np.eye(dim, dtype=_np.complex128) + summands = [] + for i in range(dim - n_leak): + s = _np.outer(I[i], I[i]) + s_superket = _bt.stdmx_to_stdvec(s).ravel() + summands.append(s_superket) + psi = _np.sum(summands, axis=0) / _np.sqrt(dim - n_leak) + proj_psi = _np.outer(psi, psi) + proj_psi_superket = _bt.stdmx_to_stdvec(proj_psi).ravel() + + idle_gate = _np.eye(dim**2, dtype=_np.complex128) + a_tensor_idle = _np.kron(a, idle_gate) + invb_tensor_idle = _np.kron(_np.linalg.inv(b), idle_gate) + temp1 = a_tensor_idle @ proj_psi_superket + temp2 = invb_tensor_idle @ temp1 + ent_fid = proj_psi_superket.conj() @ temp2 + + temp2_as_density = _bt.stdvec_to_stdmx(temp2) + b_tensor_idle = _np.kron(b, idle_gate) + ent_fid2 = fidelity(proj_psi, temp2_as_density) + + return ent_fid + + def average_gate_fidelity(a, b, mx_basis='pp', is_tp=None, is_unitary=None): """ Computes the average gate fidelity (AGF) between two gates. From 9bf8bf6ffaa892800bf0f2c6257b89e29bd703a5 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Wed, 28 Feb 2024 15:45:05 -0500 Subject: [PATCH 07/51] extend implementation to support bases other than pauli-product. Add lots of in-line comments explaining the logic. --- pygsti/tools/optools.py | 110 +++++++++++++++++++--------------------- 1 file changed, 53 insertions(+), 57 deletions(-) diff --git a/pygsti/tools/optools.py b/pygsti/tools/optools.py index 078750948..37bb0cdcc 100644 --- a/pygsti/tools/optools.py +++ b/pygsti/tools/optools.py @@ -23,6 +23,7 @@ from pygsti.tools import jamiolkowski as _jam from pygsti.tools import lindbladtools as _lt from pygsti.tools import matrixtools as _mt +from pygsti.baseobjs import basis as _pgb from pygsti.baseobjs.basis import Basis as _Basis, ExplicitBasis as _ExplicitBasis, DirectSumBasis as _DirectSumBasis from pygsti.baseobjs.label import Label as _Label from pygsti.baseobjs.errorgenlabel import LocalElementaryErrorgenLabel as _LocalElementaryErrorgenLabel @@ -518,72 +519,62 @@ def is_tp_fn(x): return fidelity(JA, JB) -def simple_entanglement_fidelity(a_, b_, mx_basis, n_leak=0): - dim = int(_np.sqrt(a_.shape[0])) - assert a_.shape == (dim**2, dim**2) - assert b_.shape == (dim**2, dim**2) +def simple_entanglement_fidelity(op_a, op_b, mx_basis, n_leak=0): + dim = int(_np.sqrt(op_a.shape[0])) + assert op_a.shape == (dim**2, dim**2) + assert op_b.shape == (dim**2, dim**2) - # set psi = (sum_i |ii>)/sqrt(dim), - # where |ii> is the tensor product of the i^th basis vector in C^dim with itself. - I = _np.eye(dim, dtype=_np.complex128) - summands = [] - for i in range(dim - n_leak): - s = _np.outer(I[i], I[i]) - # s = |i> \tensor \tensor |i> === |ii> - # ... unclear if that's certain to be the case. - summands.append(s_superket) - psi = _np.sum(summands, axis=0) / _np.sqrt(dim - n_leak) - proj_psi = _np.outer(psi, psi) - # ^ That is absolutely positively correct! - proj_psi_superket = _bt.stdmx_to_ppvec(proj_psi).ravel() - # ^ Robin thinks things are breaking here. - # The ordering of basis vectors isn't respecting tensor - # product structure. That is, proj_psi_superket - # doesn't have the same vectorized tensor product - # structure convention as kron (or kron with the order changed). - # - # TODO: figure out how to order basis elements for stdvec - # so that tensor product structure in stdmx is preserved - # in stdvec. + # op_a and op_b act on a space "S" that's isomorphic to density matrices of a dim-level system. + # + # We care about op_a and op_b only up to their action on the subspace + # U = {rho in S : = 0 for all i >= dim - n_leak }. # - # X = [X11, X12] , size(X11) == size(X22) - # [X21, X22] + # It's easier to talk about this subspace (and related subspaces) if op_a and op_b are in + # the standard basis. So the first thing we do is convert to that basis. + std_basis = _pgb.BuiltinBasis('std', dim**2) + op_a = _mt.change_basis(op_a, mx_basis, std_basis) + op_b = _mt.change_basis(op_b, mx_basis, std_basis) + + # Our next step is to construct lifted operators "lift_op_a" and "lift_op_b" that act on the + # tensor product space S2 = (S \otimes S) according to the identities # - # basis(space(X)) = [basis(space(X11)), basis(space(X12)),..., basis(space(X22))] - # ^ Concatenation. - # ^ The vectorization of the matrix commutes with tensor product. + # lift_op_a( sigma \otimes rho ) = op_a(sigma) \otimes rho + # lift_op_b( sigma \otimes rho ) = op_b(sigma) \otimes rho # - # ^ The tensor product of the bases for matrix units of X11, X12. + # for all sigma, rho in S. The way we do this implicitly fixes a basis for S2 as the + # tensor product basis (std_basis \otimes std_basis). We'll make that explicit later on. + idle_gate = _np.eye(dim**2, dtype=_np.complex128) + lift_op_a = _np.kron(op_a, idle_gate) + lift_op_b = _np.kron(op_b, idle_gate) + + # Now we'll compare these lifted operators by how they act on specific state in S2. + # That state is rho_mm = |psi> = (|00> + |11> + ... + |dim - n_leak - 1>) / sqrt(dim - n_leak). # + # The "mm" in "rho_mm" stands for "maximally mixed." + I = _np.eye(dim, dtype=_np.complex128) + summands = [] + for i in range(dim - n_leak): + ket_i = I[i] + temp = _np.outer(ket_i, ket_i) + ket_ii = _bt.stdmx_to_stdvec(temp).ravel() + summands.append(ket_ii) + psi = _np.sum(summands, axis=0) / _np.sqrt(dim - n_leak) + rho_mm = _np.outer(psi, psi) + + # Of course, lift_op_a and lift_op_b only act on states in their superket representations. + # We need the superket representation of rho_mm in terms of the tensor product basis for S2. # - """ - We are currently taking two 2-dimensional vector spaces (U and V) and: - Tensoring them together to get a 4-d space W - Constructing a basis of matrix units for the 16-d space B(W) - Instead we should take U and V and: - Construct a basis B1 of matrix units for the 4-d space B(U) - Construct a basis B2 of matrix units for the 4-d space B(V) - Tensor B1 \otimes B2 to get a 16-d space. + # Luckily, pyGSTi has a class for generating bases for a tensor-product space given + # bases for the constituent spaces appearing in the tensor product. + ten_basis = _pgb.TensorProdBasis((std_basis, std_basis)) + rho_mm_superket = _bt.stdmx_to_vec(rho_mm, ten_basis).ravel() - ^ Super easy to handle with powers of 2 and Pauli-Products. - ... NOT super easy to handle in ... literally anything else. - (but also not crazy hard. Just be careful.) - """ - - a = _mt.change_basis(a_, mx_basis, 'pp') - b = _mt.change_basis(b_, mx_basis, 'pp') - idle_gate = _np.eye(dim**2, dtype=_np.complex128) - a_tensor_idle = _np.kron(a, idle_gate) - b_tensor_idle = _np.kron(b, idle_gate) - temp1 = a_tensor_idle @ proj_psi_superket - temp2 = b_tensor_idle @ proj_psi_superket + temp1 = lift_op_a @ rho_mm_superket + temp2 = lift_op_b @ rho_mm_superket ent_fid = _np.real(temp1.conj() @ temp2) + return ent_fid @@ -604,6 +595,11 @@ def simple_entanglement_fidelity2(a_, b_, mx_basis, n_leak=0): summands.append(s_superket) psi = _np.sum(summands, axis=0) / _np.sqrt(dim - n_leak) proj_psi = _np.outer(psi, psi) + + # I have a density matrix representation I want. + # I need to convert it into a vector whose coefficients + # are interpreted w.r.t. a tensor product basis. + proj_psi_superket = _bt.stdmx_to_stdvec(proj_psi).ravel() idle_gate = _np.eye(dim**2, dtype=_np.complex128) From 87232fce79c5b4466b7d5dc809113102d713f2d9 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Wed, 28 Feb 2024 15:51:23 -0500 Subject: [PATCH 08/51] remove unused and broken function that I only added for debugging purposes --- pygsti/tools/optools.py | 38 -------------------------------------- 1 file changed, 38 deletions(-) diff --git a/pygsti/tools/optools.py b/pygsti/tools/optools.py index 37bb0cdcc..7727c191a 100644 --- a/pygsti/tools/optools.py +++ b/pygsti/tools/optools.py @@ -578,44 +578,6 @@ def simple_entanglement_fidelity(op_a, op_b, mx_basis, n_leak=0): return ent_fid -def simple_entanglement_fidelity2(a_, b_, mx_basis, n_leak=0): - import numpy.linalg as la - dim = int(_np.sqrt(a_.shape[0])) - assert a_.shape == (dim**2, dim**2) - assert b_.shape == (dim**2, dim**2) - a = _mt.change_basis(a_, mx_basis, 'std') - b = _mt.change_basis(b_, mx_basis, 'std') - # set psi = (sum_i |ii>)/sqrt(dim), - # where |ii> is the tensor product of the i^th basis vector in C^dim with itself. - I = _np.eye(dim, dtype=_np.complex128) - summands = [] - for i in range(dim - n_leak): - s = _np.outer(I[i], I[i]) - s_superket = _bt.stdmx_to_stdvec(s).ravel() - summands.append(s_superket) - psi = _np.sum(summands, axis=0) / _np.sqrt(dim - n_leak) - proj_psi = _np.outer(psi, psi) - - # I have a density matrix representation I want. - # I need to convert it into a vector whose coefficients - # are interpreted w.r.t. a tensor product basis. - - proj_psi_superket = _bt.stdmx_to_stdvec(proj_psi).ravel() - - idle_gate = _np.eye(dim**2, dtype=_np.complex128) - a_tensor_idle = _np.kron(a, idle_gate) - invb_tensor_idle = _np.kron(_np.linalg.inv(b), idle_gate) - temp1 = a_tensor_idle @ proj_psi_superket - temp2 = invb_tensor_idle @ temp1 - ent_fid = proj_psi_superket.conj() @ temp2 - - temp2_as_density = _bt.stdvec_to_stdmx(temp2) - b_tensor_idle = _np.kron(b, idle_gate) - ent_fid2 = fidelity(proj_psi, temp2_as_density) - - return ent_fid - - def average_gate_fidelity(a, b, mx_basis='pp', is_tp=None, is_unitary=None): """ Computes the average gate fidelity (AGF) between two gates. From f39dcfbd1d5535b7eab90040d8e1603bf53c350b Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Wed, 28 Feb 2024 15:54:02 -0500 Subject: [PATCH 09/51] comment clarification --- pygsti/tools/optools.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pygsti/tools/optools.py b/pygsti/tools/optools.py index 7727c191a..c0d026a56 100644 --- a/pygsti/tools/optools.py +++ b/pygsti/tools/optools.py @@ -524,7 +524,8 @@ def simple_entanglement_fidelity(op_a, op_b, mx_basis, n_leak=0): assert op_a.shape == (dim**2, dim**2) assert op_b.shape == (dim**2, dim**2) - # op_a and op_b act on a space "S" that's isomorphic to density matrices of a dim-level system. + # op_a and op_b act on the smallest real-linear space "S" that contains density matrices + # for a dim-level system. # # We care about op_a and op_b only up to their action on the subspace # U = {rho in S : = 0 for all i >= dim - n_leak }. From 2a5d7370a4ba765404512831f4e598805274f18c Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Mon, 4 Mar 2024 09:38:00 -0500 Subject: [PATCH 10/51] clean up implementation of ExplicitOpModelCalc.residuals --- pygsti/models/explicitcalc.py | 67 ++++++++++++++--------------------- 1 file changed, 26 insertions(+), 41 deletions(-) diff --git a/pygsti/models/explicitcalc.py b/pygsti/models/explicitcalc.py index 8cd54681b..11a93a9c8 100644 --- a/pygsti/models/explicitcalc.py +++ b/pygsti/models/explicitcalc.py @@ -241,47 +241,32 @@ def residuals(self, other_calc, transform_mx=None, item_weights=None): sqrt_itemWeights = {k: _np.sqrt(v) for k, v in item_weights.items()} opWeight = sqrt_itemWeights.get('gates', 1.0) spamWeight = sqrt_itemWeights.get('spam', 1.0) - - if T is not None: - Ti = _np.linalg.inv(T) # TODO: generalize inverse op (call T.inverse() if T were a "transform" object?) - for opLabel, gate in self.operations.items(): - wt = sqrt_itemWeights.get(opLabel, opWeight) - resids.append( - wt * gate.residuals( - other_calc.operations[opLabel], T, Ti)) - nSummands += wt**2 * (gate.dim)**2 - - for lbl, rhoV in self.preps.items(): - wt = sqrt_itemWeights.get(lbl, spamWeight) - resids.append( - wt * rhoV.residuals(other_calc.preps[lbl], T, Ti)) - nSummands += wt**2 * rhoV.dim - - for lbl, Evec in self.effects.items(): - wt = sqrt_itemWeights.get(lbl, spamWeight) - resids.append( - wt * Evec.residuals(other_calc.effects[lbl], T, Ti)) - - nSummands += wt**2 * Evec.dim - - else: - for opLabel, gate in self.operations.items(): - wt = sqrt_itemWeights.get(opLabel, opWeight) - resids.append( - wt * gate.residuals(other_calc.operations[opLabel])) - nSummands += wt**2 * (gate.dim)**2 - - for lbl, rhoV in self.preps.items(): - wt = sqrt_itemWeights.get(lbl, spamWeight) - resids.append( - wt * rhoV.residuals(other_calc.preps[lbl])) - nSummands += wt**2 * rhoV.dim - - for lbl, Evec in self.effects.items(): - wt = sqrt_itemWeights.get(lbl, spamWeight) - resids.append( - wt * Evec.residuals(other_calc.effects[lbl])) - nSummands += wt**2 * Evec.dim + # if T is None: + # T = _np.eye(self.dim) + Ti = None if T is None else _np.linalg.pinv(T) + # ^ TODO: generalize inverse op (call T.inverse() if T were a "transform" object?) + + for opLabel, gate in self.operations.items(): + wt = sqrt_itemWeights.get(opLabel, opWeight) + other_gate = other_calc.operations[opLabel] + resid = wt * gate.residuals(other_gate, T, Ti) + resids.append(resid) + nSummands += wt**2 * (gate.dim)**2 + + for lbl, rhoV in self.preps.items(): + wt = sqrt_itemWeights.get(lbl, spamWeight) + other_prep = other_calc.preps[lbl] + resid = wt * rhoV.residuals(other_prep, T, Ti) + resids.append(resid) + nSummands += wt**2 * rhoV.dim + + for lbl, Evec in self.effects.items(): + wt = sqrt_itemWeights.get(lbl, spamWeight) + other_effect = other_calc.effects[lbl] + resid = wt * Evec.residuals(other_effect, T, Ti) + resids.append(resid) + + nSummands += wt**2 * Evec.dim resids = [r.ravel() for r in resids] resids = _np.concatenate(resids) From 2a7019d5e870380515c132f4ee7b2234a602035c Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Mon, 4 Mar 2024 09:40:25 -0500 Subject: [PATCH 11/51] slightly simplify minimize(...) function. There was an unnecessary branch that omitted jac as a keyword argument to scipy.minimize, even if jac was None. --- pygsti/optimize/optimize.py | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/pygsti/optimize/optimize.py b/pygsti/optimize/optimize.py index 7411eb4d1..a6a6a982f 100644 --- a/pygsti/optimize/optimize.py +++ b/pygsti/optimize/optimize.py @@ -159,11 +159,7 @@ def _basin_callback(x, f, accept): elif method == "L-BFGS-B": opts['gtol'] = opts['ftol'] = tol # gradient norm and fractional y-tolerance elif method == "Nelder-Mead": opts['maxfev'] = maxfev # max fn evals (note: ftol and xtol can also be set) - if method in ("BFGS", "CG", "Newton-CG", "L-BFGS-B", "TNC", "SLSQP", "dogleg", "trust-ncg"): # use jacobian - solution = _spo.minimize(fn, x0, options=opts, method=method, tol=tol, callback=callback, jac=jac) - else: - solution = _spo.minimize(fn, x0, options=opts, method=method, tol=tol, callback=callback) - + solution = _spo.minimize(fn, x0, options=opts, method=method, tol=tol, callback=callback, jac=jac) return solution From e38f9a9184b0700749875f8b23c6d809fdd7f1ae Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Mon, 4 Mar 2024 09:42:51 -0500 Subject: [PATCH 12/51] factor out the logic needed to set up calculation of leaky entanglement fidelity from the calculation of that metric itself. This makes it easier to reuse the setup code in other metrics. While making this change, also make the setup more efficient by using TensorProduct bases from the beginning. --- pygsti/tools/optools.py | 24 ++++++++++++++---------- 1 file changed, 14 insertions(+), 10 deletions(-) diff --git a/pygsti/tools/optools.py b/pygsti/tools/optools.py index c0d026a56..aba8f6cc2 100644 --- a/pygsti/tools/optools.py +++ b/pygsti/tools/optools.py @@ -519,7 +519,9 @@ def is_tp_fn(x): return fidelity(JA, JB) -def simple_entanglement_fidelity(op_a, op_b, mx_basis, n_leak=0): +def lift_and_act_on_maxmixed_state(op_a, op_b, mx_basis, n_leak=0): + # Note: this function is only really useful for gates on a single system (qubit, qutrit, qudit); + # not tensor products of such systems. dim = int(_np.sqrt(op_a.shape[0])) assert op_a.shape == (dim**2, dim**2) assert op_b.shape == (dim**2, dim**2) @@ -554,14 +556,11 @@ def simple_entanglement_fidelity(op_a, op_b, mx_basis, n_leak=0): # |psi> = (|00> + |11> + ... + |dim - n_leak - 1>) / sqrt(dim - n_leak). # # The "mm" in "rho_mm" stands for "maximally mixed." - I = _np.eye(dim, dtype=_np.complex128) - summands = [] - for i in range(dim - n_leak): - ket_i = I[i] - temp = _np.outer(ket_i, ket_i) - ket_ii = _bt.stdmx_to_stdvec(temp).ravel() - summands.append(ket_ii) - psi = _np.sum(summands, axis=0) / _np.sqrt(dim - n_leak) + temp = _np.eye(dim, dtype=_np.complex128) + if n_leak > 0: + temp[-n_leak:,-n_leak:] = 0.0 + temp /= _np.sqrt(dim - n_leak) + psi = _bt.stdmx_to_stdvec(temp).ravel() rho_mm = _np.outer(psi, psi) # Of course, lift_op_a and lift_op_b only act on states in their superket representations. @@ -574,8 +573,13 @@ def simple_entanglement_fidelity(op_a, op_b, mx_basis, n_leak=0): temp1 = lift_op_a @ rho_mm_superket temp2 = lift_op_b @ rho_mm_superket - ent_fid = _np.real(temp1.conj() @ temp2) + return temp1, temp2, ten_basis + + +def leaky_entanglement_fidelity(op_a, op_b, mx_basis, n_leak=0): + temp1, temp2, _ = lift_and_act_on_maxmixed_state(op_a, op_b, mx_basis, n_leak) + ent_fid = _np.real(temp1.conj() @ temp2) return ent_fid From 7cbc59f8bcac838109393c540359d1bb6488cc4d Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Mon, 4 Mar 2024 09:49:19 -0500 Subject: [PATCH 13/51] add implementation of leaky jtracedist --- pygsti/tools/optools.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/pygsti/tools/optools.py b/pygsti/tools/optools.py index aba8f6cc2..6fa74adca 100644 --- a/pygsti/tools/optools.py +++ b/pygsti/tools/optools.py @@ -583,6 +583,14 @@ def leaky_entanglement_fidelity(op_a, op_b, mx_basis, n_leak=0): return ent_fid +def leaky_jtracedist(op_a, op_b, mx_basis, n_leak=0): + temp1, temp2, ten_basis = lift_and_act_on_maxmixed_state(op_a, op_b, mx_basis, n_leak) + temp1_std = _bt.vec_to_stdmx(temp1, ten_basis, keep_complex=True) + temp2_std = _bt.basistools.vec_to_stdmx(temp2, ten_basis, keep_complex=True) + j_dist = tracedist(temp1_std, temp2_std) + return j_dist + + def average_gate_fidelity(a, b, mx_basis='pp', is_tp=None, is_unitary=None): """ Computes the average gate fidelity (AGF) between two gates. From 8e9c5dfcbb43056aa01bca1320e8fbedbb33a60e Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Mon, 4 Mar 2024 09:52:43 -0500 Subject: [PATCH 14/51] enable gauge optimization with leakage-aware metrics using non-LS optimizers. Right now only entanglement fidelity and jtracedist have implementations that can consider leakage. (So there are no leakage-aware metrics for SPAM.) --- pygsti/algorithms/gaugeopt.py | 96 ++++++++++++++++++----------------- 1 file changed, 49 insertions(+), 47 deletions(-) diff --git a/pygsti/algorithms/gaugeopt.py b/pygsti/algorithms/gaugeopt.py index b2137e99b..7cdeed1a6 100644 --- a/pygsti/algorithms/gaugeopt.py +++ b/pygsti/algorithms/gaugeopt.py @@ -29,7 +29,7 @@ def gaugeopt_to_target(model, target_model, item_weights=None, gauge_group=None, method='auto', maxiter=100000, maxfev=None, tol=1e-8, oob_check_interval=0, convert_model_to=None, return_all=False, comm=None, - verbosity=0, check_jac=False): + verbosity=0, check_jac=False, n_leak=0): """ Optimize the gauge degrees of freedom of a model to that of a target. @@ -170,7 +170,7 @@ def gaugeopt_to_target(model, target_model, item_weights=None, objective_fn, jacobian_fn = _create_objective_fn( model, target_model, item_weights, cptp_penalty_factor, spam_penalty_factor, - gates_metric, spam_metric, method, comm, check_jac) + gates_metric, spam_metric, method, comm, check_jac, n_leak) result = gaugeopt_custom(model, objective_fn, gauge_group, method, maxiter, maxfev, tol, oob_check_interval, @@ -303,9 +303,6 @@ def _call_jacobian_fn(gauge_group_el_vec): printer.log("--- Gauge Optimization (%s method, %s) ---" % (method, str(type(gauge_group))), 2) if method == 'ls': - #minSol = _opt.least_squares(_call_objective_fn, x0, #jac=_call_jacobian_fn, - # max_nfev=maxfev, ftol=tol) - #solnX = minSol.x assert(_call_jacobian_fn is not None), "Cannot use 'ls' method unless jacobian is available" ralloc = _baseobjs.ResourceAllocation(comm) # FUTURE: plumb up a resource alloc object? test_f = _call_objective_fn(x0) @@ -353,7 +350,7 @@ def _call_jacobian_fn(gauge_group_el_vec): def _create_objective_fn(model, target_model, item_weights=None, cptp_penalty_factor=0, spam_penalty_factor=0, gates_metric="frobenius", spam_metric="frobenius", - method=None, comm=None, check_jac=False): + method=None, comm=None, check_jac=False, n_leak=0): """ Creates the objective function and jacobian (if available) for gaugeopt_to_target @@ -387,9 +384,9 @@ def _transform_with_oob_check(mdl, gauge_group_el, oob_check): assert(gates_metric.startswith("frobenius") and spam_metric.startswith("frobenius")), \ "Only 'frobenius' and 'frobeniustt' metrics can be used when `method='ls'`!" assert(gates_metric == spam_metric) - frobenius_transform_target = bool(gates_metric == 'frobeniustt') # tt = "transform target" + frobenius_transform_target = bool(gates_metric == 'frobeniustt') # tt = "transform target"; this is rarely True - if frobenius_transform_target: + if frobenius_transform_target: # Riley note: this is rare! full_target_model = target_model.copy() full_target_model.convert_members_inplace("full") # so we can gauge-transform the target model. else: @@ -404,7 +401,7 @@ def _objective_fn(gauge_group_el, oob_check): transformed = _transform_with_oob_check(model, gauge_group_el, oob_check) other = target_model - residuals, _ = transformed.residuals(other, None, item_weights) + residuals, _ = transformed.residuals(other, None, item_weights) # < Metrics will be computed. # We still the non-target model to be transformed and checked for these penalties if cptp_penalty_factor > 0 or spam_penalty_factor > 0: @@ -456,31 +453,33 @@ def _jacobian_fn(gauge_group_el): jacMx = _np.zeros((L, N)) - #Overview of terms: - # objective: op_term = (S_inv * gate * S - target_op) - # jac: d(op_term) = (d (S_inv) * gate * S + S_inv * gate * dS ) - # d(op_term) = (-(S_inv * dS * S_inv) * gate * S + S_inv * gate * dS ) + """ + Overview of terms: + objective: op_term = (S_inv * gate * S - target_op) + jac: d(op_term) = (d (S_inv) * gate * S + S_inv * gate * dS ) + d(op_term) = (-(S_inv * dS * S_inv) * gate * S + S_inv * gate * dS ) - # objective: rho_term = (S_inv * rho - target_rho) - # jac: d(rho_term) = d (S_inv) * rho - # d(rho_term) = -(S_inv * dS * S_inv) * rho + objective: rho_term = (S_inv * rho - target_rho) + jac: d(rho_term) = d (S_inv) * rho + d(rho_term) = -(S_inv * dS * S_inv) * rho - # objective: ET_term = (E.T * S - target_E.T) - # jac: d(ET_term) = E.T * dS + objective: ET_term = (E.T * S - target_E.T) + jac: d(ET_term) = E.T * dS - #Overview of terms when frobenius_transform_target == True). Note that the objective - #expressions are identical to the above except for an additional overall minus sign and S <=> S_inv. + Overview of terms when frobenius_transform_target == True). Note that the objective + expressions are identical to the above except for an additional overall minus sign and S <=> S_inv. - # objective: op_term = (gate - S * target_op * S_inv) - # jac: d(op_term) = -(dS * target_op * S_inv + S * target_op * -(S_inv * dS * S_inv) ) - # d(op_term) = (-dS * target_op * S_inv + S * target_op * (S_inv * dS * S_inv) ) + objective: op_term = (gate - S * target_op * S_inv) + jac: d(op_term) = -(dS * target_op * S_inv + S * target_op * -(S_inv * dS * S_inv) ) + d(op_term) = (-dS * target_op * S_inv + S * target_op * (S_inv * dS * S_inv) ) - # objective: rho_term = (rho - S * target_rho) - # jac: d(rho_term) = - dS * target_rho + objective: rho_term = (rho - S * target_rho) + jac: d(rho_term) = - dS * target_rho - # objective: ET_term = (E.T - target_E.T * S_inv) - # jac: d(ET_term) = - target_E.T * -(S_inv * dS * S_inv) - # d(ET_term) = target_E.T * (S_inv * dS * S_inv) + objective: ET_term = (E.T - target_E.T * S_inv) + jac: d(ET_term) = - target_E.T * -(S_inv * dS * S_inv) + d(ET_term) = target_E.T * (S_inv * dS * S_inv) + """ #Distribute computation across processors allDerivColSlice = slice(0, N) @@ -590,6 +589,8 @@ def _mock_objective_fn(v): else: # non-least-squares case where objective function returns a single float # and (currently) there's no analytic jacobian + if n_leak > 0: + pass def _objective_fn(gauge_group_el, oob_check): mdl = _transform_with_oob_check(model, gauge_group_el, oob_check) @@ -636,14 +637,16 @@ def _objective_fn(gauge_group_el, oob_check): elif gates_metric == "fidelity": for opLbl in mdl.operations: wt = item_weights.get(opLbl, opWeight) - ret += wt * (1.0 - _tools.entanglement_fidelity( - target_model.operations[opLbl], mdl.operations[opLbl]))**2 + top = target_model.operations[opLbl].to_dense() + mop = mdl.operations[opLbl].to_dense() + ret += wt * (1.0 - _tools.leaky_entanglement_fidelity(top, mop, mxBasis, n_leak))**2 elif gates_metric == "tracedist": for opLbl in mdl.operations: wt = item_weights.get(opLbl, opWeight) - ret += opWeight * _tools.jtracedist( - target_model.operations[opLbl], mdl.operations[opLbl]) + top = target_model.operations[opLbl].to_dense() + mop = mdl.operations[opLbl].to_dense() + ret += wt * _tools.leaky_jtracedist(top, mop, mxBasis, n_leak) else: raise ValueError("Invalid gates_metric: %s" % gates_metric) @@ -664,30 +667,29 @@ def _objective_fn(gauge_group_el, oob_check): ret += transformed_target.frobeniusdist(model, None, wts) elif spam_metric == "fidelity": - for preplabel, prep in mdl.preps.items(): + for preplabel, m_prep in mdl.preps.items(): wt = item_weights.get(preplabel, spamWeight) - rhoMx1 = _tools.vec_to_stdmx(prep, mxBasis) - rhoMx2 = _tools.vec_to_stdmx( - target_model.preps[preplabel], mxBasis) + rhoMx1 = _tools.vec_to_stdmx(m_prep.to_dense(), mxBasis) + t_prep = target_model.preps[preplabel] + rhoMx2 = _tools.vec_to_stdmx(t_prep.to_dense(), mxBasis) ret += wt * (1.0 - _tools.fidelity(rhoMx1, rhoMx2))**2 - - for povmlabel, povm in mdl.povms.items(): + + for povmlabel in mdl.povms.keys(): wt = item_weights.get(povmlabel, spamWeight) - ret += wt * (1.0 - _tools.povm_fidelity( - mdl, target_model, povmlabel))**2 + fidelity = _tools.povm_fidelity(mdl, target_model, povmlabel) + ret += wt * (1.0 - fidelity)**2 elif spam_metric == "tracedist": - for preplabel, prep in mdl.preps.items(): + for preplabel, m_prep in mdl.preps.items(): wt = item_weights.get(preplabel, spamWeight) - rhoMx1 = _tools.vec_to_stdmx(prep, mxBasis) - rhoMx2 = _tools.vec_to_stdmx( - target_model.preps[preplabel], mxBasis) + rhoMx1 = _tools.vec_to_stdmx(m_prep.to_dense(), mxBasis) + t_prep = target_model.preps[preplabel] + rhoMx2 = _tools.vec_to_stdmx(t_prep.to_dense(), mxBasis) ret += wt * _tools.tracedist(rhoMx1, rhoMx2) - for povmlabel, povm in mdl.povms.items(): + for povmlabel in mdl.povms.keys(): wt = item_weights.get(povmlabel, spamWeight) - ret += wt * (1.0 - _tools.povm_jtracedist( - mdl, target_model, povmlabel))**2 + ret += wt * _tools.povm_jtracedist(mdl, target_model, povmlabel) else: raise ValueError("Invalid spam_metric: %s" % spam_metric) From 0113730751fb0727b428f4dba04e557d0a2790b9 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Mon, 4 Mar 2024 12:03:42 -0500 Subject: [PATCH 15/51] support for leakage-aware Frobenius distance with non-LS optimizer. Added annotations to all gauage optimization objectives used in the non-LS optimizer, indicating if that particular objective has support for leakage-aware metrics. --- pygsti/algorithms/gaugeopt.py | 88 ++++++++++++++++++++--------------- pygsti/models/explicitcalc.py | 60 ++++++++++++++---------- 2 files changed, 87 insertions(+), 61 deletions(-) diff --git a/pygsti/algorithms/gaugeopt.py b/pygsti/algorithms/gaugeopt.py index 7cdeed1a6..b19c9aa15 100644 --- a/pygsti/algorithms/gaugeopt.py +++ b/pygsti/algorithms/gaugeopt.py @@ -590,19 +590,32 @@ def _mock_objective_fn(v): # non-least-squares case where objective function returns a single float # and (currently) there's no analytic jacobian if n_leak > 0: - pass + dim = int(_np.sqrt(mxBasis.dim)) + B = _tools.leading_dxd_submatrix_basis_vectors(dim - n_leak, dim, mxBasis) + P = B @ B.T.conj() + if _np.linalg.norm(P.imag) > 1e-12: + raise ValueError() + else: + P = P.real + transform_mx_arg = (P, None) + # ^ The semantics of this tuple are defined by the frobeniusdist function + # in the ExplicitOpModelCalc class. There are intended semantics for + # the second element of the tuple, but those aren't implemented yet so + # for now I'm setting the second entry to None. -- Riley + else: + transform_mx_arg = None + + assert gates_metric != "frobeniustt" + assert spam_metric != "frobeniustt" + # ^ Erik and Corey said these are rarely used. I've removed support for + # them in this codepath (non-LS optimizer) in order to make it easier to + # read my updated code for leakage-aware metrics. It wouldn't be hard to + # add support back, but I just want to keep things simple. -- Riley def _objective_fn(gauge_group_el, oob_check): mdl = _transform_with_oob_check(model, gauge_group_el, oob_check) ret = 0 - if gates_metric == "frobeniustt" or spam_metric == "frobeniustt": - full_target_model = target_model.copy() - full_target_model.convert_members_inplace("full") # so we can gauge-transform the target model. - transformed_target = _transform_with_oob_check(full_target_model, gauge_group_el.inverse(), oob_check) - else: - transformed_target = None - if cptp_penalty_factor > 0: mdl.basis = mxBasis # set basis for jamiolkowski iso cpPenaltyVec = _cptp_penalty(mdl, cptp_penalty_factor, mdl.basis) @@ -614,27 +627,24 @@ def _objective_fn(gauge_group_el, oob_check): ret += _np.sum(spamPenaltyVec) if target_model is not None: - if gates_metric == "frobenius": - if spam_metric == "frobenius": - ret += mdl.frobeniusdist(target_model, None, item_weights) - else: - wts = item_weights.copy(); wts['spam'] = 0.0 - for k in wts: - if k in mdl.preps or \ - k in mdl.povms: wts[k] = 0.0 - ret += mdl.frobeniusdist(target_model, None, wts) - - elif gates_metric == "frobeniustt": - if spam_metric == "frobeniustt": - ret += transformed_target.frobeniusdist(model, None, item_weights) + # Leakage-aware metric supported, per implementation in mdl.frobeniusdist. + # Refer to how mdl.frobeniusdist handles the case when transform_mx_arg + # is a tuple in order to understand how the leakage-aware metric is defined. + if "frobenius" in gates_metric: + if spam_metric == gates_metric: + val = mdl.frobeniusdist(target_model, transform_mx_arg, item_weights) else: wts = item_weights.copy(); wts['spam'] = 0.0 for k in wts: if k in mdl.preps or \ k in mdl.povms: wts[k] = 0.0 - ret += transformed_target.frobeniusdist(model, None, wts) + val = mdl.frobeniusdist(target_model, transform_mx_arg, wts) + if "squared" in gates_metric: + val = val ** 2 + ret += val elif gates_metric == "fidelity": + # Leakage-aware metric supported, using leaky_entanglement_fidelity. for opLbl in mdl.operations: wt = item_weights.get(opLbl, opWeight) top = target_model.operations[opLbl].to_dense() @@ -642,6 +652,7 @@ def _objective_fn(gauge_group_el, oob_check): ret += wt * (1.0 - _tools.leaky_entanglement_fidelity(top, mop, mxBasis, n_leak))**2 elif gates_metric == "tracedist": + # Leakage-aware metric supported, using leaky_jtracedist. for opLbl in mdl.operations: wt = item_weights.get(opLbl, opWeight) top = target_model.operations[opLbl].to_dense() @@ -650,23 +661,25 @@ def _objective_fn(gauge_group_el, oob_check): else: raise ValueError("Invalid gates_metric: %s" % gates_metric) - if spam_metric == "frobenius": - if gates_metric != "frobenius": # otherwise handled above to match normalization in frobeniusdist - wts = item_weights.copy(); wts['gates'] = 0.0 - for k in wts: - if k in mdl.operations or \ - k in mdl.instruments: wts[k] = 0.0 - ret += mdl.frobeniusdist(target_model, None, wts) - - elif spam_metric == "frobeniustt": - if gates_metric != "frobeniustt": # otherwise handled above to match normalization in frobeniusdist - wts = item_weights.copy(); wts['gates'] = 0.0 - for k in wts: - if k in mdl.operations or \ - k in mdl.instruments: wts[k] = 0.0 - ret += transformed_target.frobeniusdist(model, None, wts) + if "frobenius" in spam_metric and gates_metric == spam_metric: + # We already handled SPAM error in this case. Just return. + return ret + + elif "frobenius" in spam_metric: + # Leakage-aware metric supported in principle via implementation in + # mdl.frobeniusdist (check implementation to see how it handles the + # case when transform_mx_arg is a tuple). + wts = item_weights.copy(); wts['gates'] = 0.0 + for k in wts: + if k in mdl.operations or \ + k in mdl.instruments: wts[k] = 0.0 + val = mdl.frobeniusdist(target_model, transform_mx_arg, wts) + if "squared" in spam_metric: + val = val ** 2 + ret += val elif spam_metric == "fidelity": + # Leakage-aware metrics NOT available for preplabel, m_prep in mdl.preps.items(): wt = item_weights.get(preplabel, spamWeight) rhoMx1 = _tools.vec_to_stdmx(m_prep.to_dense(), mxBasis) @@ -680,6 +693,7 @@ def _objective_fn(gauge_group_el, oob_check): ret += wt * (1.0 - fidelity)**2 elif spam_metric == "tracedist": + # Leakage-aware metrics NOT available. for preplabel, m_prep in mdl.preps.items(): wt = item_weights.get(preplabel, spamWeight) rhoMx1 = _tools.vec_to_stdmx(m_prep.to_dense(), mxBasis) diff --git a/pygsti/models/explicitcalc.py b/pygsti/models/explicitcalc.py index 11a93a9c8..466669d35 100644 --- a/pygsti/models/explicitcalc.py +++ b/pygsti/models/explicitcalc.py @@ -120,12 +120,18 @@ def frobeniusdist(self, other_calc, transform_mx=None, other_calc : ForwardSimulator the other gate calculator to difference against. - transform_mx : numpy array, optional - if not None, transform this model by - G => inv(transform_mx) * G * transform_mx, for each operation matrix G - (and similar for rho and E vectors) before taking the difference. - This transformation is applied only for the difference and does - not alter the values stored in this model. + transform_mx : numpy array or tuple, optional + if transform_mx is a numpy array, then for each operation matrix + G we implicitly consider the transformed quantity + G => inv(transform_mx) * G * transform_mx + Similar transformations are applied for effect vectors. + This transformation is applied only for the difference and does not + alter the values stored in this model. + + if transform_mx is a tuple then it should consist of two numpy arrays, + the first of which will be interpreted as transform_mx in the usual + sense and the second of which will either be None or will be syntactically + subtituted for inv(transform_mx). item_weights : dict, optional Dictionary of weighting factors for individual gates and spam @@ -146,18 +152,22 @@ def frobeniusdist(self, other_calc, transform_mx=None, ------- float """ - d = 0; T = transform_mx + if isinstance(transform_mx, tuple): + T, Ti = transform_mx + else: + T = transform_mx + Ti = None if T is None else _np.linalg.pinv(T) + d = 0 nSummands = 0.0 if item_weights is None: item_weights = {} opWeight = item_weights.get('gates', 1.0) spamWeight = item_weights.get('spam', 1.0) - if T is not None: - Ti = _np.linalg.inv(T) # TODO: generalize inverse op (call T.inverse() if T were a "transform" object?) + if (T is None and Ti is None) or isinstance(transform_mx, _np.ndarray): + # use the original implementation. for opLabel, gate in self.operations.items(): wt = item_weights.get(opLabel, opWeight) - d += wt * gate.frobeniusdist_squared( - other_calc.operations[opLabel], T, Ti) + d += wt * gate.frobeniusdist_squared(other_calc.operations[opLabel], T, Ti) nSummands += wt * (gate.dim)**2 for lbl, rhoV in self.preps.items(): @@ -169,28 +179,32 @@ def frobeniusdist(self, other_calc, transform_mx=None, wt = item_weights.get(lbl, spamWeight) d += wt * Evec.frobeniusdist_squared(other_calc.effects[lbl], T, Ti) nSummands += wt * Evec.dim - else: + # we're in the special case that I'm creating. for opLabel, gate in self.operations.items(): wt = item_weights.get(opLabel, opWeight) - d += wt * gate.frobeniusdist_squared(other_calc.operations[opLabel]) + gate_mx = gate.to_dense() + other_mx = other_calc.operations[opLabel].to_dense() + delta = gate_mx - other_mx + if T is not None: + delta = delta @ T + if Ti is not None: + delta = Ti @ delta + val = _np.linalg.norm(delta.flatten()) + d += wt * val**2 nSummands += wt * (gate.dim)**2 for lbl, rhoV in self.preps.items(): + # keep the original implementation, for now. wt = item_weights.get(lbl, spamWeight) - d += wt * rhoV.frobeniusdist_squared(other_calc.preps[lbl]) + d += wt * rhoV.frobeniusdist_squared(other_calc.preps[lbl], None, None) nSummands += wt * rhoV.dim for lbl, Evec in self.effects.items(): + # keep the original implementation, for now. wt = item_weights.get(lbl, spamWeight) - d += wt * Evec.frobeniusdist_squared(other_calc.effects[lbl]) - nSummands += wt * Evec.dim - - #Temporary: check that this function can be computed by - # calling residuals - replace with this later. - resids, chk_nSummands = self.residuals(other_calc, transform_mx, item_weights) - assert(_np.isclose(_np.sum(resids**2), d)) - assert(_np.isclose(chk_nSummands, nSummands)) + d += wt * Evec.frobeniusdist_squared(other_calc.effects[lbl], None, None) + nSummands += wt * Evec.dimxw if normalize and nSummands > 0: return _np.sqrt(d / nSummands) @@ -241,8 +255,6 @@ def residuals(self, other_calc, transform_mx=None, item_weights=None): sqrt_itemWeights = {k: _np.sqrt(v) for k, v in item_weights.items()} opWeight = sqrt_itemWeights.get('gates', 1.0) spamWeight = sqrt_itemWeights.get('spam', 1.0) - # if T is None: - # T = _np.eye(self.dim) Ti = None if T is None else _np.linalg.pinv(T) # ^ TODO: generalize inverse op (call T.inverse() if T were a "transform" object?) From 885856aca76c2fe2bd4bf9cb29debe17d6151f02 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Mon, 4 Mar 2024 12:23:56 -0500 Subject: [PATCH 16/51] left out a function needed for code in the last commit to work --- pygsti/tools/optools.py | 81 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 81 insertions(+) diff --git a/pygsti/tools/optools.py b/pygsti/tools/optools.py index 6fa74adca..c1bf2ba70 100644 --- a/pygsti/tools/optools.py +++ b/pygsti/tools/optools.py @@ -591,6 +591,87 @@ def leaky_jtracedist(op_a, op_b, mx_basis, n_leak=0): return j_dist +def leading_dxd_submatrix_basis_vectors(d: int, n: int, current_basis): + """ + Let "H" denote n^2 dimensional Hilbert-Schdmit space, and let "U" denote the d^2 + dimensional subspace of H spanned by vectors whose Hermitian matrix representations + are zero outside the leading d-by-d submatrix. + + This function returns a column-unitary matrix "B" where P = B B^{\dagger} is the + orthogonal projector from H to U with respect to current_basis. We return B rather + than P only because it's simpler to get P from B than it is to get B from P. + + See below for this function's original use-case. + + Raison d'etre + ------------- + Suppose we canonically measure the distance between two process matrices (M1, M2) by + + D(M1, M2; H) = max || (M1 - M2) v || + v is in H, (Eq. 1) + tr(v) = 1, + v is positive + + for some norm || * ||. Suppose also that we want an analog of this distance when + (M1, M2) are restricted to the linear subspace U consisting of all vectors in H + whose matrix representations are zero outside of their leading d-by-d submatrix. + + One natural way to do this is via the function D(M1, M2; U) -- i.e., just replace + H in (Eq. 1) with the subspace U. Using P to denote the orthogonal projector onto U, + we claim that we can evaluate this function via the identity + + D(M1, M2; U) = D(M1 P, M2 P; H). (Eq. 2) + + To see why this is the case, consider a positive vector v and its projection u = P v. + Since a vector is called positive whenever its Hermitian matrix representation is positive + semidefinite (PSD), we need to show that u is positive. This can be seen by considering + block 2-by-2 partitions of the matrix representations of (u,v), where the leading block + is d-by-d: + + mat(v) = [x11, x12] and mat(u) = [x11, 0] + [x21, x22] [ 0, 0]. + + In particular, u is positive if and only if x11 is PSD, and x11 must be PSD for v + to be positive. Furthermore, positivity of v requires that x22 is PSD, which implies + + 0 <= tr(u) = tr(x11) <= tr(v). + + Given this, it is easy to establish (Eq 2.) by considering how the following pair + of problems have the same optimal objective function value + + max || (M1 - M2) P v || and max || (M1 - M2) P v || + mat(v) = [x11, x12] mat(v) = [x11, x12] + [x21, x22] [x21, x22] + mat(v) is PSD x11 is PSD + tr(x11) + tr(x22) = 1 tr(x11) <= 1. + + In fact, this can be taken a little further! The whole argument goes through unchanged + if, instead of starting with the objective function || (M1 - M2) v ||, we started with + f((M1 - M2) v) and f satisfied the property that f(c v) >= f(v) whenever c is a scalar + greater than or equal to one. + + Interesting idea: + Set M2 = 0. + Use || (I - P) M1 P || as a metric for leakage. + Use || P M1 (I - P) || as a metric for seepage + """ + assert d <= n + current_basis = _pgb.Basis.cast(current_basis, dim=n**2) + std_to_current = current_basis.create_transform_matrix('std') + if d == n: + return std_to_current + # we have to select a proper subset of columns in current_basis + std_basis = _pgb.BuiltinBasis(name='std', dim_or_statespace=n**2) + label2ind = {ell: idx for idx,ell in enumerate(std_basis.labels)} + basis_ind = [] + for i in range(d): + for j in range(d): + basis_ind.append(label2ind[f"({i},{j})"]) + basis_ind = _np.array(basis_ind) + submatrix_basis_vectors = std_to_current[:, basis_ind] + return submatrix_basis_vectors + + def average_gate_fidelity(a, b, mx_basis='pp', is_tp=None, is_unitary=None): """ Computes the average gate fidelity (AGF) between two gates. From 040970540dd9b81bb2ba62bc97686098234540d1 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Mon, 4 Mar 2024 12:27:53 -0500 Subject: [PATCH 17/51] tests for the leaky_entanglement_fidelity function that Ive had for a while. The tests only consider when no leakage is modeled. New tests wil be needed for when theres a leakage dimension --- test/unit/tools/test_optools.py | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/test/unit/tools/test_optools.py b/test/unit/tools/test_optools.py index 23c8c36fa..f69286b7c 100644 --- a/test/unit/tools/test_optools.py +++ b/test/unit/tools/test_optools.py @@ -412,6 +412,19 @@ def test_entanglement_fidelity(self): self.assertAlmostEqual(fidelity_TP_unitary_jam, expect) self.assertAlmostEqual(fidelity_TP_unitary_std, expect) + def test_leaky_entanglement_fidelity(self): + fidelity_TP_unitary= ot.leaky_entanglement_fidelity(self.A_TP, self.B_unitary, 'pp') + fidelity_TP_unitary_no_flag= ot.leaky_entanglement_fidelity(self.A_TP, self.B_unitary, 'pp') + fidelity_TP_unitary_jam= ot.leaky_entanglement_fidelity(self.A_TP, self.B_unitary, 'pp') + fidelity_TP_unitary_std= ot.leaky_entanglement_fidelity(self.A_TP_std, self.B_unitary_std, mx_basis='std') + + expect = 0.4804724656092404 + self.assertAlmostEqual(fidelity_TP_unitary, expect) + self.assertAlmostEqual(fidelity_TP_unitary_no_flag, expect) + self.assertAlmostEqual(fidelity_TP_unitary_jam, expect) + self.assertAlmostEqual(fidelity_TP_unitary_std, expect) + pass + def test_fidelity_upper_bound(self): np.random.seed(0) Q = np.linalg.qr(np.random.randn(4,4) + 1j*np.random.randn(4,4))[0] From 0e8eea7350188908ae76a9e7532f252d2285668c Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Mon, 4 Mar 2024 13:44:30 -0500 Subject: [PATCH 18/51] fix syntax mistake --- pygsti/tools/optools.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pygsti/tools/optools.py b/pygsti/tools/optools.py index c1bf2ba70..b9aacb1db 100644 --- a/pygsti/tools/optools.py +++ b/pygsti/tools/optools.py @@ -586,7 +586,7 @@ def leaky_entanglement_fidelity(op_a, op_b, mx_basis, n_leak=0): def leaky_jtracedist(op_a, op_b, mx_basis, n_leak=0): temp1, temp2, ten_basis = lift_and_act_on_maxmixed_state(op_a, op_b, mx_basis, n_leak) temp1_std = _bt.vec_to_stdmx(temp1, ten_basis, keep_complex=True) - temp2_std = _bt.basistools.vec_to_stdmx(temp2, ten_basis, keep_complex=True) + temp2_std = _bt.vec_to_stdmx(temp2, ten_basis, keep_complex=True) j_dist = tracedist(temp1_std, temp2_std) return j_dist From 44dcc94387a597cd3803e47724cc9f8ccf97e571 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Mon, 4 Mar 2024 17:19:47 -0500 Subject: [PATCH 19/51] rename test_gauageopt.py file to indicate that its contents only test that various gauge optimization functions don`t raise errors. --- .../algorithms/{test_gaugeopt.py => test_gaugeopt_noerrors.py} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename test/unit/algorithms/{test_gaugeopt.py => test_gaugeopt_noerrors.py} (100%) diff --git a/test/unit/algorithms/test_gaugeopt.py b/test/unit/algorithms/test_gaugeopt_noerrors.py similarity index 100% rename from test/unit/algorithms/test_gaugeopt.py rename to test/unit/algorithms/test_gaugeopt_noerrors.py From c3818d2113e58b4ff1da80cd80e106b75fa69628 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Tue, 19 Mar 2024 15:13:59 -0400 Subject: [PATCH 20/51] correctness test for gauge optimization (ls-based and L-BFGS-based) --- .../algorithms/test_gaugeopt_correctness.py | 223 ++++++++++++++++++ 1 file changed, 223 insertions(+) create mode 100644 test/unit/algorithms/test_gaugeopt_correctness.py diff --git a/test/unit/algorithms/test_gaugeopt_correctness.py b/test/unit/algorithms/test_gaugeopt_correctness.py new file mode 100644 index 000000000..af4b39c70 --- /dev/null +++ b/test/unit/algorithms/test_gaugeopt_correctness.py @@ -0,0 +1,223 @@ + +from pygsti.models.gaugegroup import TPGaugeGroup, TPGaugeGroupElement, \ + UnitaryGaugeGroup, UnitaryGaugeGroupElement,\ + FullGaugeGroup, FullGaugeGroupElement +from ..util import BaseCase + +import time +from collections import OrderedDict +import numpy as np +from pygsti.modelpacks import smq1Q_XYI +import scipy.linalg as la +import pygsti.tools.optools as pgo +import pygsti.baseobjs.label as pgl +import pygsti.baseobjs.basisconstructors as pgbc +import pygsti.algorithms.gaugeopt as gop + + + +def gate_metrics_dict(model, target): + metrics = {'infids': OrderedDict(), 'frodists': OrderedDict(), 'tracedists': OrderedDict()} + for lbl in model.operations.keys(): + model_gate = model[lbl] + target_gate = target[lbl] + metrics['infids'][lbl] = pgo.entanglement_infidelity(model_gate, target_gate) + metrics['frodists'][lbl] = pgo.frobeniusdist(model_gate, target_gate) + metrics['tracedists'][lbl] = pgo.tracedist(model_gate, target_gate) + return metrics + + +def check_gate_metrics_are_nontrivial(metrics, tol): + """ + "metrics" is a dict of the kind produced by gate_metrics_dict(model, target). + + This function checks if all values in that dict are strictly greater than tol, + with the exception of those corresponding to the idle gate. + + This makes sure that gauge optimization for "model" matching "target" is a + nontrivial problem. + """ + idle_label = pgl.Label(tuple()) + for lbl, val in metrics['infids'].items(): + if lbl == idle_label: + continue + assert val > tol, f"{val} is at most {tol}, failure for gate {lbl} w.r.t. infidelity." + for lbl, val in metrics['frodists'].items(): + if lbl == idle_label: + continue + assert val > tol, f"{val} is at most {tol}, failure for gate {lbl} w.r.t. frobenius distance." + for lbl, val in metrics['tracedists'].items(): + if lbl == idle_label: + continue + assert val > tol, f"{val} is at most {tol}, failure for gate {lbl} w.r.r. trace distance." + return + + +def check_gate_metrics_near_zero(metrics, tol): + """ + "metrics" is a dict of the kind produced by gate_metrics_dict(model, target). + + This function should be called to check if gauge optimization for "model" + matching "target" succeeded. + """ + for lbl, val in metrics['infids'].items(): + assert val <= tol, f"{val} exceeds {tol}, failure for gate {lbl} w.r.t infidelity." + for lbl, val in metrics['frodists'].items(): + assert val <= tol, f"{val} exceeds {tol}, failure for gate {lbl} w.r.t. frobenius distance." + for lbl, val in metrics['tracedists'].items(): + assert val <= tol, f"{val} exceeds {tol}, failure for gate {lbl} w.r.t. trace distance." + return + + + +class sm1QXYI_FindPerfectGauge: + + @property + def default_gauge_group(self): + """ + Instance of an object that subclasses GaugeGroup. + The state space for this gauge group is set based + on that of "self.target". + """ + return self._default_gauge_group + + @property + def gauge_grp_el_class(self): + """ + Handle for some subclass of GaugeGroupElement. + Must be compatible with self.default_gauge_group. + """ + return self._gauge_grp_el_class + + @property + def gauge_space_name(self) -> str: + """ + Only used for printing purposes in "self._main_tester". + """ + return self._gauge_space_name + + def setUp(self): + np.random.seed(0) + target = smq1Q_XYI.target_model() + self.target = target.depolarize(op_noise=0.01, spam_noise=0.001) + self.model = None + self._default_gauge_group = None + self._gauge_grp_el_class = None + self._gauge_space_name = '' + self.N_REPS = 1 + raise NotImplementedError() + + def _gauge_transform_model(self, seed): + self.model = self.target.copy() + self.model.default_gauge_group = self.default_gauge_group + np.random.seed(seed) + self.U = la.expm(np.random.randn()/2 * -1j * (pgbc.sigmax + pgbc.sigmaz)/np.sqrt(2)) + self.gauge_grp_el = self.gauge_grp_el_class(pgo.unitary_to_pauligate(self.U)) + self.model.transform_inplace(self.gauge_grp_el) + self.metrics_before = gate_metrics_dict(self.model, self.target) + check_gate_metrics_are_nontrivial(self.metrics_before, tol=1e-2) + return + + def _main_tester(self, method, seed, test_tol, alg_tol, gop_objective): + assert gop_objective in {'frobenius', 'frobenius_squared', 'tracedist', 'fidelity'} + self._gauge_transform_model(seed) + tic = time.time() + newmodel = gop.gaugeopt_to_target(self.model, self.target, method=method, tol=alg_tol, spam_metric=gop_objective, gates_metric=gop_objective) + toc = time.time() + dt = toc - tic + metrics_after = gate_metrics_dict(newmodel, self.target) + check_gate_metrics_near_zero(metrics_after, test_tol) + return dt + + def test_lbfgs_frodist(self): + self.setUp() + times = [] + for seed in range(self.N_REPS): + dt = self._main_tester('L-BFGS-B', seed, test_tol=1e-6, alg_tol=1e-8, gop_objective='frobenius') + times.append(dt) + print(f'GaugeOpt over {self.gauge_space_name} w.r.t. Frobenius dist, using L-BFGS: {times}.') + return + + def test_lbfgs_tracedist(self): + self.setUp() + times = [] + for seed in range(self.N_REPS): + dt = self._main_tester('L-BFGS-B', seed, test_tol=1e-6, alg_tol=1e-8, gop_objective='tracedist') + times.append(dt) + print(f'GaugeOpt over {self.gauge_space_name} w.r.t. trace dist, using L-BFGS: {times}.') + return + + def test_ls(self): + self.setUp() + times = [] + for seed in range(self.N_REPS): + dt = self._main_tester('ls', seed, test_tol=1e-6, alg_tol=1e-15, gop_objective='frobenius') + times.append(dt) + print(f'GaugeOpt over {self.gauge_space_name} w.r.t. Frobenius dist, using LS : {times}.') + return + + pass + + +""" +This class' test for minimizing trace distance with L-BFGS takes excessively long (about 5 seconds). +""" +class FindPerfectGauge_TPGroup_Tester(BaseCase, sm1QXYI_FindPerfectGauge): + + def setUp(self): + np.random.seed(0) + target = smq1Q_XYI.target_model() + self.target = target.depolarize(op_noise=0.01, spam_noise=0.001) + self.model = None + self._default_gauge_group = TPGaugeGroup(self.target.state_space, self.target.basis) + self._gauge_grp_el_class = TPGaugeGroupElement + self._gauge_space_name = 'TPGaugeGroup' + self.N_REPS = 1 + return + + def gauge_transform_model(self, seed): + sm1QXYI_FindPerfectGauge.gauge_transform_model(self, seed) + # ^ That sets self.model + return + + +class FindPerfectGauge_UnitaryGroup_Tester(BaseCase, sm1QXYI_FindPerfectGauge): + + def setUp(self): + np.random.seed(0) + target = smq1Q_XYI.target_model() + self.target = target.depolarize(op_noise=0.01, spam_noise=0.001) + self.model = None + self._default_gauge_group = UnitaryGaugeGroup(self.target.state_space, self.target.basis) + self._gauge_grp_el_class = UnitaryGaugeGroupElement + self._gauge_space_name = 'UnitaryGaugeGroup' + self.N_REPS = 3 + # ^ This is fast, so we can afford more tests. + return + + def gauge_transform_model(self, seed): + sm1QXYI_FindPerfectGauge.gauge_transform_model(self, seed) + # ^ That sets self.model + return + + +""" +This class' test for minimizing trace distance with L-BFGS takes excessively long (about 5 seconds). +""" +class FindPerfectGauge_FullGroup_Tester(BaseCase, sm1QXYI_FindPerfectGauge): + + def setUp(self): + np.random.seed(0) + target = smq1Q_XYI.target_model() + self.target = target.depolarize(op_noise=0.01, spam_noise=0.001) + self.model = None + self._default_gauge_group = FullGaugeGroup(self.target.state_space, self.target.basis) + self._gauge_grp_el_class = FullGaugeGroupElement + self._gauge_space_name = 'FullGaugeGroup' + self.N_REPS = 1 + return + + def gauge_transform_model(self, seed): + sm1QXYI_FindPerfectGauge.gauge_transform_model(self, seed) + # ^ That sets self.model + return From 89d8c63d39759caeb88f8b147e86b3f696d23134 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Tue, 19 Mar 2024 15:53:46 -0400 Subject: [PATCH 21/51] fix typo --- pygsti/models/explicitcalc.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pygsti/models/explicitcalc.py b/pygsti/models/explicitcalc.py index 466669d35..530197cf3 100644 --- a/pygsti/models/explicitcalc.py +++ b/pygsti/models/explicitcalc.py @@ -204,7 +204,7 @@ def frobeniusdist(self, other_calc, transform_mx=None, # keep the original implementation, for now. wt = item_weights.get(lbl, spamWeight) d += wt * Evec.frobeniusdist_squared(other_calc.effects[lbl], None, None) - nSummands += wt * Evec.dimxw + nSummands += wt * Evec.dim if normalize and nSummands > 0: return _np.sqrt(d / nSummands) From 5f25fa2fe82e21f503aa2be74381e3f42d5cb23d Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Tue, 19 Mar 2024 15:55:32 -0400 Subject: [PATCH 22/51] remove some awkward casting that it turns out wasnt needed --- pygsti/tools/basistools.py | 10 ---------- 1 file changed, 10 deletions(-) diff --git a/pygsti/tools/basistools.py b/pygsti/tools/basistools.py index a380f2dfd..2ec896cd5 100644 --- a/pygsti/tools/basistools.py +++ b/pygsti/tools/basistools.py @@ -141,13 +141,6 @@ def change_basis(mx, from_basis, to_basis): The given operation matrix converted to the `to_basis` basis. Array size is the same as `mx`. """ - if isinstance(mx, _np.ndarray) and mx.ndim == 0: - # mx is probably a wrapper around a pyGSTi object. - mx = mx.item() - if hasattr(mx, 'to_dense'): - # Question: how do we know the returned representation is - # in fact in the "from_basis"? - mx = mx.to_dense() if len(mx.shape) not in (1, 2): raise ValueError("Invalid dimension of object - must be 1 or 2, i.e. a vector or matrix") @@ -279,9 +272,6 @@ def create_basis_for_matrix(mx, basis): ------- Basis """ - if isinstance(mx, _np.ndarray) and mx.ndim == 0: - # ^ mx is probably just holding a pyGSTi object - mx = mx.item() dim = mx.shape[0] if isinstance(basis, _basis.Basis): assert(basis.dim == dim), "Supplied Basis has wrong dimension!" From 1f05402fc417b41721a64c456591dbc41a3c9128 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Tue, 2 Apr 2024 10:13:31 -0400 Subject: [PATCH 23/51] remove unnecessary comment changes --- pygsti/algorithms/gaugeopt.py | 46 +++++++++++++++++------------------ 1 file changed, 22 insertions(+), 24 deletions(-) diff --git a/pygsti/algorithms/gaugeopt.py b/pygsti/algorithms/gaugeopt.py index b19c9aa15..d60f75c76 100644 --- a/pygsti/algorithms/gaugeopt.py +++ b/pygsti/algorithms/gaugeopt.py @@ -384,9 +384,9 @@ def _transform_with_oob_check(mdl, gauge_group_el, oob_check): assert(gates_metric.startswith("frobenius") and spam_metric.startswith("frobenius")), \ "Only 'frobenius' and 'frobeniustt' metrics can be used when `method='ls'`!" assert(gates_metric == spam_metric) - frobenius_transform_target = bool(gates_metric == 'frobeniustt') # tt = "transform target"; this is rarely True + frobenius_transform_target = bool(gates_metric == 'frobeniustt') # tt = "transform target" - if frobenius_transform_target: # Riley note: this is rare! + if frobenius_transform_target: full_target_model = target_model.copy() full_target_model.convert_members_inplace("full") # so we can gauge-transform the target model. else: @@ -401,7 +401,7 @@ def _objective_fn(gauge_group_el, oob_check): transformed = _transform_with_oob_check(model, gauge_group_el, oob_check) other = target_model - residuals, _ = transformed.residuals(other, None, item_weights) # < Metrics will be computed. + residuals, _ = transformed.residuals(other, None, item_weights) # We still the non-target model to be transformed and checked for these penalties if cptp_penalty_factor > 0 or spam_penalty_factor > 0: @@ -453,33 +453,31 @@ def _jacobian_fn(gauge_group_el): jacMx = _np.zeros((L, N)) - """ - Overview of terms: - objective: op_term = (S_inv * gate * S - target_op) - jac: d(op_term) = (d (S_inv) * gate * S + S_inv * gate * dS ) - d(op_term) = (-(S_inv * dS * S_inv) * gate * S + S_inv * gate * dS ) + #Overview of terms: + # objective: op_term = (S_inv * gate * S - target_op) + # jac: d(op_term) = (d (S_inv) * gate * S + S_inv * gate * dS ) + # d(op_term) = (-(S_inv * dS * S_inv) * gate * S + S_inv * gate * dS ) - objective: rho_term = (S_inv * rho - target_rho) - jac: d(rho_term) = d (S_inv) * rho - d(rho_term) = -(S_inv * dS * S_inv) * rho + # objective: rho_term = (S_inv * rho - target_rho) + # jac: d(rho_term) = d (S_inv) * rho + # d(rho_term) = -(S_inv * dS * S_inv) * rho - objective: ET_term = (E.T * S - target_E.T) - jac: d(ET_term) = E.T * dS + # objective: ET_term = (E.T * S - target_E.T) + # jac: d(ET_term) = E.T * dS - Overview of terms when frobenius_transform_target == True). Note that the objective - expressions are identical to the above except for an additional overall minus sign and S <=> S_inv. + #Overview of terms when frobenius_transform_target == True). Note that the objective + #expressions are identical to the above except for an additional overall minus sign and S <=> S_inv. - objective: op_term = (gate - S * target_op * S_inv) - jac: d(op_term) = -(dS * target_op * S_inv + S * target_op * -(S_inv * dS * S_inv) ) - d(op_term) = (-dS * target_op * S_inv + S * target_op * (S_inv * dS * S_inv) ) + # objective: op_term = (gate - S * target_op * S_inv) + # jac: d(op_term) = -(dS * target_op * S_inv + S * target_op * -(S_inv * dS * S_inv) ) + # d(op_term) = (-dS * target_op * S_inv + S * target_op * (S_inv * dS * S_inv) ) - objective: rho_term = (rho - S * target_rho) - jac: d(rho_term) = - dS * target_rho + # objective: rho_term = (rho - S * target_rho) + # jac: d(rho_term) = - dS * target_rho - objective: ET_term = (E.T - target_E.T * S_inv) - jac: d(ET_term) = - target_E.T * -(S_inv * dS * S_inv) - d(ET_term) = target_E.T * (S_inv * dS * S_inv) - """ + # objective: ET_term = (E.T - target_E.T * S_inv) + # jac: d(ET_term) = - target_E.T * -(S_inv * dS * S_inv) + # d(ET_term) = target_E.T * (S_inv * dS * S_inv) #Distribute computation across processors allDerivColSlice = slice(0, N) From 943e1119b2a15c8fc6c355abeb892b7d4dad8c09 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Mon, 15 Apr 2024 11:54:49 -0400 Subject: [PATCH 24/51] change how leakage affects distances between SPAM model members for Frobenius distance --- pygsti/models/explicitcalc.py | 31 +++++++++++++++++-------------- 1 file changed, 17 insertions(+), 14 deletions(-) diff --git a/pygsti/models/explicitcalc.py b/pygsti/models/explicitcalc.py index 530197cf3..15e15ed23 100644 --- a/pygsti/models/explicitcalc.py +++ b/pygsti/models/explicitcalc.py @@ -153,31 +153,31 @@ def frobeniusdist(self, other_calc, transform_mx=None, float """ if isinstance(transform_mx, tuple): - T, Ti = transform_mx + P, invP = transform_mx else: - T = transform_mx - Ti = None if T is None else _np.linalg.pinv(T) + P = transform_mx + invP = None if P is None else _np.linalg.pinv(P) d = 0 nSummands = 0.0 if item_weights is None: item_weights = {} opWeight = item_weights.get('gates', 1.0) spamWeight = item_weights.get('spam', 1.0) - if (T is None and Ti is None) or isinstance(transform_mx, _np.ndarray): + if (P is None and invP is None) or isinstance(transform_mx, _np.ndarray): # use the original implementation. for opLabel, gate in self.operations.items(): wt = item_weights.get(opLabel, opWeight) - d += wt * gate.frobeniusdist_squared(other_calc.operations[opLabel], T, Ti) + d += wt * gate.frobeniusdist_squared(other_calc.operations[opLabel], P, invP) nSummands += wt * (gate.dim)**2 for lbl, rhoV in self.preps.items(): wt = item_weights.get(lbl, spamWeight) - d += wt * rhoV.frobeniusdist_squared(other_calc.preps[lbl], T, Ti) + d += wt * rhoV.frobeniusdist_squared(other_calc.preps[lbl], P, invP) nSummands += wt * rhoV.dim for lbl, Evec in self.effects.items(): wt = item_weights.get(lbl, spamWeight) - d += wt * Evec.frobeniusdist_squared(other_calc.effects[lbl], T, Ti) + d += wt * Evec.frobeniusdist_squared(other_calc.effects[lbl], P, invP) nSummands += wt * Evec.dim else: # we're in the special case that I'm creating. @@ -186,24 +186,27 @@ def frobeniusdist(self, other_calc, transform_mx=None, gate_mx = gate.to_dense() other_mx = other_calc.operations[opLabel].to_dense() delta = gate_mx - other_mx - if T is not None: - delta = delta @ T - if Ti is not None: - delta = Ti @ delta + if P is not None: + delta = delta @ P + if invP is not None: + delta = invP @ delta val = _np.linalg.norm(delta.flatten()) d += wt * val**2 nSummands += wt * (gate.dim)**2 for lbl, rhoV in self.preps.items(): - # keep the original implementation, for now. wt = item_weights.get(lbl, spamWeight) d += wt * rhoV.frobeniusdist_squared(other_calc.preps[lbl], None, None) nSummands += wt * rhoV.dim for lbl, Evec in self.effects.items(): - # keep the original implementation, for now. wt = item_weights.get(lbl, spamWeight) - d += wt * Evec.frobeniusdist_squared(other_calc.effects[lbl], None, None) + evec = Evec.to_dense() + other = other_calc.effects[lbl].to_dense() + delta = evec - other + delta = delta @ P + val = _np.linalg.norm(delta.flatten()) + d += wt * val**2 nSummands += wt * Evec.dim if normalize and nSummands > 0: From 6e9de5c8885d4c25ee0080ede77cfb937b8a3efb Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Mon, 15 Apr 2024 12:18:59 -0400 Subject: [PATCH 25/51] changes in README for MFT_20230125 from Corey --- pygsti/baseobjs/basisconstructors.py | 8 ++++++-- pygsti/protocols/gst.py | 3 +++ 2 files changed, 9 insertions(+), 2 deletions(-) diff --git a/pygsti/baseobjs/basisconstructors.py b/pygsti/baseobjs/basisconstructors.py index 5fa8f3a19..0706ef6bd 100644 --- a/pygsti/baseobjs/basisconstructors.py +++ b/pygsti/baseobjs/basisconstructors.py @@ -862,8 +862,12 @@ def _is_integer(x): nQubits = _np.log2(matrix_dim) if not _is_integer(nQubits): - raise ValueError( - "Dimension for Pauli tensor product matrices must be an integer *power of 2* (not %d)" % matrix_dim) + nQutrits = _np.log(matrix_dim) / _np.log(3) + if _is_integer(nQutrits): + return gm_matrices(matrix_dim) + else: + raise ValueError( + "Dimension for Pauli tensor product matrices must be an integer *power of 2* (not %d)" % matrix_dim) nQubits = int(round(nQubits)) if nQubits == 0: # special case: return single 1x1 identity mx diff --git a/pygsti/protocols/gst.py b/pygsti/protocols/gst.py index 64c015635..17323bc3f 100644 --- a/pygsti/protocols/gst.py +++ b/pygsti/protocols/gst.py @@ -2453,6 +2453,9 @@ def _compute_1d_reference_values_and_name(estimate, badfit_options, gaugeopt_sui spamdd[key] = 0.5 * _tools.optools.povm_diamonddist(gaugeopt_model, target_model, key) dd[lbl]['SPAM'] = sum(spamdd.values()) + for k,v in spamdd.items(): + dd[lbl][k] = v + return dd, 'diamond distance' else: raise ValueError("Invalid wildcard1d_reference value (%s) in bad-fit options!" From 703c3ff20dc66fb5871c0cccba46d79224f1b07e Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Fri, 19 Apr 2024 07:27:44 -0400 Subject: [PATCH 26/51] tiny docstring change --- pygsti/report/factory.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/pygsti/report/factory.py b/pygsti/report/factory.py index 8d2f675d7..a4d573b5c 100644 --- a/pygsti/report/factory.py +++ b/pygsti/report/factory.py @@ -1177,8 +1177,7 @@ def construct_standard_report(results, title="auto", Returns ------- - Workspace - The workspace object used to create the report + Report """ printer = _VerbosityPrinter.create_printer(verbosity, comm=comm) From 3cbfda5a397b355295e6c3b8c0c6bb6e5259b50c Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Tue, 7 May 2024 16:22:28 -0400 Subject: [PATCH 27/51] very messy support for subspace-restricted error metrics --- pygsti/report/factory.py | 10 +++++++++- pygsti/report/reportables.py | 35 ++++++++++++++++++++++++++++++++++ pygsti/report/section/gauge.py | 6 +++++- 3 files changed, 49 insertions(+), 2 deletions(-) diff --git a/pygsti/report/factory.py b/pygsti/report/factory.py index a4d573b5c..04e49c434 100644 --- a/pygsti/report/factory.py +++ b/pygsti/report/factory.py @@ -1184,6 +1184,10 @@ def construct_standard_report(results, title="auto", ws = ws or _ws.Workspace() advanced_options = advanced_options or {} + n_leak = advanced_options.get('n_leak', 0) + # ^ Not sure if this is a good place to pass this parameter. + # A better implementation would be to store this info in model.basis + # (for the various models). linlogPercentile = advanced_options.get('linlog percentile', 5) nmthreshold = advanced_options.get('nmthreshold', DEFAULT_NONMARK_ERRBAR_THRESHOLD) embed_figures = advanced_options.get('embed_figures', True) @@ -1316,8 +1320,12 @@ def construct_standard_report(results, title="auto", 'gauge_opt_labels': tuple(gauge_opt_labels), 'max_lengths': tuple(Ls), 'switchbd_maxlengths': tuple(swLs), - 'show_unmodeled_error': bool('ShowUnmodeledError' in flags) + 'show_unmodeled_error': bool('ShowUnmodeledError' in flags), + 'n_leak' : n_leak } + # ^ Not sure if this is a good place to pass n_leak. + # A better implementation would be to store this info in model.basis + # (for the various models). templates = dict( html='~standard_html_report', diff --git a/pygsti/report/reportables.py b/pygsti/report/reportables.py index 1a0cf4152..d4d501187 100644 --- a/pygsti/report/reportables.py +++ b/pygsti/report/reportables.py @@ -985,6 +985,14 @@ def maximum_trace_dist(gate, mx_basis): Maximum_trace_dist = _modf.opfn_factory(maximum_trace_dist) # init args == (model, op_label) +def leaky_maximum_trace_dist(gate, mx_basis): + closestUOpMx = _alg.find_closest_unitary_opmx(gate) + _tools.jamiolkowski_iso(closestUOpMx, mx_basis, mx_basis) + n_leak = 1 + return _tools.leaky_jtracedist(gate, closestUOpMx, mx_basis, n_leak) + +Leaky_maximum_trace_dist = _modf.opfn_factory(leaky_maximum_trace_dist) + def angles_btwn_rotn_axes(model): """ @@ -1056,6 +1064,12 @@ def entanglement_fidelity(a, b, mx_basis): Entanglement_fidelity = _modf.opsfn_factory(entanglement_fidelity) # init args == (model1, model2, op_label) +def leaky_entanglement_fidelity(a, b, mx_basis): + n_leak = 1 + return _tools.leaky_entanglement_fidelity(a, b, mx_basis, n_leak) + +Leaky_entanglement_fidelity = _modf.opsfn_factory(leaky_entanglement_fidelity) + def entanglement_infidelity(a, b, mx_basis): """ @@ -1082,6 +1096,11 @@ def entanglement_infidelity(a, b, mx_basis): Entanglement_infidelity = _modf.opsfn_factory(entanglement_infidelity) # init args == (model1, model2, op_label) +def leaky_entanglement_infidelity(a, b, mx_basis): + return 1 - leaky_entanglement_fidelity(a, b, mx_basis) + +Leaky_entanglement_infidelity = _modf.opsfn_factory(leaky_entanglement_infidelity) + def closest_unitary_fidelity(a, b, mx_basis): # assume vary model1, model2 fixed """ @@ -1173,6 +1192,12 @@ def jtrace_diff(a, b, mx_basis): # assume vary model1, model2 fixed Jt_diff = _modf.opsfn_factory(jtrace_diff) # init args == (model1, model2, op_label) +def leaky_jtrace_diff(a, b, mx_basis): + n_leak = 1 + return _tools.leaky_jtracedist(a, b, mx_basis, n_leak) + +Leaky_Jt_diff = _modf.opsfn_factory(leaky_jtrace_diff) + if _CVXPY_AVAILABLE: @@ -2393,10 +2418,14 @@ def info_of_opfn_by_name(name): info = { "inf": ("Entanglement|Infidelity", "1.0 - "), + "la-inf": ("Entanglement|Infidelity (subspace)", + "TO-WRITE"), "agi": ("Avg. Gate|Infidelity", "d/(d+1) (entanglement infidelity)"), "trace": ("1/2 Trace|Distance", "0.5 | Chi(A) - Chi(B) |_tr"), + "la-trace" : ("1/2 Trace|Distance (subspace)", + "TO-WRITE"), "diamond": ("1/2 Diamond-Dist", "0.5 sup | (1 x (A-B))(rho) |_tr"), "nuinf": ("Non-unitary|Ent. Infidelity", @@ -2466,12 +2495,18 @@ def evaluate_opfn_by_name(name, model, target_model, op_label_or_string, if name == "inf": fn = Entanglement_infidelity if b else \ Circuit_entanglement_infidelity + elif name == "la-inf": + assert b + fn = Leaky_entanglement_infidelity elif name == "agi": fn = Avg_gate_infidelity if b else \ Circuit_avg_gate_infidelity elif name == "trace": fn = Jt_diff if b else \ Circuit_jt_diff + elif name == "la-trace": + assert b + fn = Leaky_Jt_diff elif name == "diamond": fn = HalfDiamondNorm if b else \ CircuitHalfDiamondNorm diff --git a/pygsti/report/section/gauge.py b/pygsti/report/section/gauge.py index 84a66908e..c9c642e81 100644 --- a/pygsti/report/section/gauge.py +++ b/pygsti/report/section/gauge.py @@ -168,9 +168,13 @@ def final_model_spam_vs_target_table(workspace, switchboard=None, confidence_lev @_Section.figure_factory(4) def final_gates_vs_target_table_gauge_var(workspace, switchboard=None, confidence_level=None, ci_brevity=1, **kwargs): + if kwargs.get('n_leak', 0) == 0: + display = ('inf', 'agi', 'trace', 'diamond', 'nuinf', 'nuagi') + else: + display = ('inf', 'la-inf', 'agi', 'trace', 'la-trace', 'diamond', 'nuinf', 'nuagi') return workspace.GatesVsTargetTable( switchboard.mdl_final, switchboard.mdl_target, _cri(1, switchboard, confidence_level, ci_brevity), - display=('inf', 'agi', 'trace', 'diamond', 'nuinf', 'nuagi') + display=display ) @_Section.figure_factory(3) From 43df6156027b199287414585d76191d3dbe8cbc0 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Tue, 24 Sep 2024 16:29:08 -0400 Subject: [PATCH 28/51] de-nest two lines --- pygsti/tools/jamiolkowski.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/pygsti/tools/jamiolkowski.py b/pygsti/tools/jamiolkowski.py index cafabb2d9..598889902 100644 --- a/pygsti/tools/jamiolkowski.py +++ b/pygsti/tools/jamiolkowski.py @@ -215,11 +215,13 @@ def fast_jamiolkowski_iso_std(operation_mx, op_mx_basis): #first, get operation matrix into std basis operation_mx = _np.asarray(operation_mx) op_mx_basis = _bt.create_basis_for_matrix(operation_mx, op_mx_basis) - opMxInStdBasis = _bt.change_basis(operation_mx, op_mx_basis, op_mx_basis.create_equivalent('std')) + temp = op_mx_basis.create_equivalent('std') + opMxInStdBasis = _bt.change_basis(operation_mx, op_mx_basis, temp) #expand operation matrix so it acts on entire space of dmDim x dmDim density matrices - opMxInStdBasis = _bt.resize_std_mx(opMxInStdBasis, 'expand', op_mx_basis.create_equivalent('std'), - op_mx_basis.create_simple_equivalent('std')) + temp1 = op_mx_basis.create_equivalent('std') + temp2 = op_mx_basis.create_simple_equivalent('std') + opMxInStdBasis = _bt.resize_std_mx(opMxInStdBasis, 'expand',temp1, temp2) #Shuffle indices to go from process matrix to Jamiolkowski matrix (they vectorize differently) N2 = opMxInStdBasis.shape[0]; N = int(_np.sqrt(N2)) From fa2aae42b3022688172434d06f83ff4c0e93bdb8 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Fri, 27 Sep 2024 16:57:13 -0400 Subject: [PATCH 29/51] Add utilities for constructing CVXPY models involving (convex) diamond-norm expressions or (concave) root-fidelity expressions. Add a "normalize" keyword argument to the various functions for computing(/uncomputing) Jamiolkowski isomorphisms of operation matrices. Added option to return all optimization model variables used in computing diamonddist. --- pygsti/tools/jamiolkowski.py | 82 ++++++++++++---- pygsti/tools/optools.py | 108 +++++---------------- pygsti/tools/sdptools.py | 167 ++++++++++++++++++++++++++++++++ test/unit/tools/test_optools.py | 44 +++++++-- 4 files changed, 288 insertions(+), 113 deletions(-) create mode 100644 pygsti/tools/sdptools.py diff --git a/pygsti/tools/jamiolkowski.py b/pygsti/tools/jamiolkowski.py index 598889902..0c9e95249 100644 --- a/pygsti/tools/jamiolkowski.py +++ b/pygsti/tools/jamiolkowski.py @@ -64,9 +64,10 @@ # J(Phi) = sum_(0 + < J(phi).dag, X.dag > ) - # Subject to [[ I otimes rho0, X ], - # [ X.dag , I otimes rho1]] >> 0 - # rho0, rho1 are density matrices - # X is linear operator - - import cvxpy as _cp - - rho0 = _cp.Variable((smallDim, smallDim), name='rho0', hermitian=True) - rho1 = _cp.Variable((smallDim, smallDim), name='rho1', hermitian=True) - X = _cp.Variable((dim, dim), name='X', complex=True) - Y = _cp.real(X) - Z = _cp.imag(X) - - K = J.real - L = J.imag - if hasattr(_cp, 'scalar_product'): - objective_expr = _cp.scalar_product(K, Y) + _cp.scalar_product(L, Z) + if return_all_vars: + return objective_val, varvals + elif return_x: + return objective_val, varvals[0] else: - Kf = K.flatten(order='F') - Yf = Y.flatten(order='F') - Lf = L.flatten(order='F') - Zf = Z.flatten(order='F') - objective_expr = Kf @ Yf + Lf @ Zf - - objective = _cp.Maximize(objective_expr) - - ident = _np.identity(smallDim, 'd') - kr_tau0 = _cp.kron(ident, _cp.imag(rho0)) - kr_tau1 = _cp.kron(ident, _cp.imag(rho1)) - kr_sig0 = _cp.kron(ident, _cp.real(rho0)) - kr_sig1 = _cp.kron(ident, _cp.real(rho1)) - - block_11 = _cp.bmat([[kr_sig0 , Y ], - [ Y.T , kr_sig1]]) - block_21 = _cp.bmat([[kr_tau0 , Z ], - [ -Z.T , kr_tau1]]) - block_12 = block_21.T - mat_joint = _cp.bmat([[block_11, block_12], - [block_21, block_11]]) - constraints = [ - mat_joint >> 0, - rho0 >> 0, - rho1 >> 0, - _cp.trace(rho0) == 1., - _cp.trace(rho1) == 1. - ] - prob = _cp.Problem(objective, constraints) - return prob, [X, rho0, rho1] + return objective_val def jtracedist(a, b, mx_basis='pp'): # Jamiolkowski trace distance: Tr(|J(a)-J(b)|) diff --git a/pygsti/tools/sdptools.py b/pygsti/tools/sdptools.py new file mode 100644 index 000000000..6652369e0 --- /dev/null +++ b/pygsti/tools/sdptools.py @@ -0,0 +1,167 @@ +""" +Functions for constructing semidefinite programming models +""" +#*************************************************************************************************** +# Copyright 2024 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. +#*************************************************************************************************** + +from __future__ import annotations +from typing import TYPE_CHECKING, List, Tuple +if TYPE_CHECKING: + import cvxpy as cp + +import numpy as np +import scipy.linalg as la + +from pygsti.tools import basistools as bt +from pygsti.tools import jamiolkowski as jam +from pygsti.tools import lindbladtools as lt +from pygsti.tools import matrixtools as mt +from pygsti.baseobjs import basis as pgb +from pygsti.baseobjs.basis import Basis, ExplicitBasis, DirectSumBasis + + +try: + import cvxpy as cp + old_cvxpy = bool(tuple(map(int, cp.__version__.split('.'))) < (1, 0)) + CVXPY_ENABLED = not old_cvxpy +except: + CVXPY_ENABLED = False + + +def diamond_norm_model_jamiolkowski(J): + # return a model for computing the diamond norm. + # + # Uses the primal SDP from arXiv:1207.5726v2, Sec 3.2 + # + # Maximize 1/2 ( < J, X > + < J.dag, X.dag > ) + # Subject to [[ I otimes rho0, X ], + # [ X.dag , I otimes rho1]] >> 0 + # rho0, rho1 are density matrices + # X is linear operator + # + # ".dag" returns the adjoint. + dim = J.shape[0] + smallDim = int(np.sqrt(dim)) + assert dim == smallDim**2 + + rho0 = cp.Variable((smallDim, smallDim), name='rho0', hermitian=True) + rho1 = cp.Variable((smallDim, smallDim), name='rho1', hermitian=True) + X = cp.Variable((dim, dim), name='X', complex=True) + Y = cp.real(X) + Z = cp.imag(X) + # = J.dag.ravel() @ X.ravel() + # = J.ravel() @ X.dag.ravel() = conj() + # + # ---> real() = 1/2 ( + ) + # ---> can skip the factor 1/2 if we just form real() directly. + # + + K = J.real + L = J.imag + if hasattr(cp, 'scalar_product'): + objective_expr = cp.scalar_product(K, Y) + cp.scalar_product(L, Z) + else: + Kf = K.flatten(order='F') + Yf = Y.flatten(order='F') + Lf = L.flatten(order='F') + Zf = Z.flatten(order='F') + objective_expr = Kf @ Yf + Lf @ Zf + + objective = cp.Maximize(objective_expr) + + ident = np.identity(smallDim, 'd') + kr_tau0 = cp.kron(ident, cp.imag(rho0)) + kr_tau1 = cp.kron(ident, cp.imag(rho1)) + kr_sig0 = cp.kron(ident, cp.real(rho0)) + kr_sig1 = cp.kron(ident, cp.real(rho1)) + + block_11 = cp.bmat([[kr_sig0 , Y ], + [ Y.T , kr_sig1]]) + block_21 = cp.bmat([[kr_tau0 , Z ], + [ -Z.T , kr_tau1]]) + block_12 = block_21.T + mat_joint = cp.bmat([[block_11, block_12], + [block_21, block_11]]) + constraints = [ + mat_joint >> 0, + rho0 >> 0, + rho1 >> 0, + cp.trace(rho0) == 1., + cp.trace(rho1) == 1. + ] + prob = cp.Problem(objective, constraints) + return prob, [X, rho0, rho1] + + +def diamond_norm_canon(arg : cp.Expression, basis) -> Tuple[cp.Expression, List[cp.Constraint]]: + """ + This more or less implements canonicalization of the nonlinear expression + \|arg\|_{\diamond} into CVXPY Constraints and a representation of its epigraph. + The canonicalization isn't quite "complete" in CVXPY's usual sense, which would + require that the epigraph is affine and that no structured variables (like + Hermitian matrices) are used. + """ + constraints = [] + d = arg.shape[0] + small_d = int(np.sqrt(d)) + assert d == small_d**2 + assert arg.shape == (d, d) + Jarg = jam.jamiolkowski_iso(arg, basis, basis, normalized=False) + Y0 = cp.Variable(shape=(d, d), hermitian=True) + Y1 = cp.Variable(shape=(d, d), hermitian=True) + bmat = cp.bmat([ + [ Y0 , -Jarg], + [-Jarg.T.conj(), Y1 ] + ]) + constraints.append(bmat >> 0) + TrX_Y0 = cp.partial_trace(Y0, [small_d, small_d], 0) + TrX_Y1 = cp.partial_trace(Y1, [small_d, small_d], 0) + expr0 = cp.lambda_max(TrX_Y0) + expr1 = cp.lambda_max(TrX_Y1) + epi = (expr0 + expr1)/2 + return epi, constraints + + +def root_fidelity_canon(sigma: cp.Expression, rho: cp.Expression) -> Tuple[cp.Expression, List[cp.Constraint]]: + """ + pyGSTi defines fidelity as + + F(sigma, rho) = tr([sigma^{1/2} rho sigma^{1/2}]^{1/2})^2. + + Others (including Neilson and Chuang, Sect. 9.2.2) define it without the + square on the trace. We'll call the unsquared version the *root fidelity,* + and denote it by + + \sqrt{F}(sigma, rho) = (F(sigma, rho))^{1/2}. + + The root fidelity is jointly concave (Neilson and Chuang, Exercise 9.19). + In fact, it admits the following semidefinite programming characterization + + \sqrt{F}(sigma, rho) = Maximize real(tr(X)) + s.t. [[sigma, X],[X.T.conj(), rho]] >> 0 + + -- see Section 7.1.3 of Killoran's PhD thesis, "Entanglement quantification + and quantum benchmarking of optical communication devices." + + This function returns a pair (expr, constraints) where expr is the hypograph + variable for \sqrt{F}(sigma, rho) and constraints is a list of CVXPY Constraint + objects used in the semidefinite representation of the hypograph. + """ + t = cp.Variable() + d = sigma.shape[0] + X = cp.Variable(shape=(d, d), complex=True) + bmat = cp.hermitian_wrap(cp.bmat([ + [ sigma, X ], + [ X.T.conj(), rho ] + ])) + constraints = [ + bmat >> 0, + cp.trace(cp.real(X)) >= t + ] + return t, constraints diff --git a/test/unit/tools/test_optools.py b/test/unit/tools/test_optools.py index f69286b7c..a94a4c280 100644 --- a/test/unit/tools/test_optools.py +++ b/test/unit/tools/test_optools.py @@ -4,12 +4,14 @@ import sys import numpy as np import scipy +import scipy.linalg as la from pygsti.baseobjs.basis import Basis from pygsti.baseobjs.errorgenlabel import LocalElementaryErrorgenLabel as LEEL import pygsti.tools.basistools as bt import pygsti.tools.lindbladtools as lt import pygsti.tools.optools as ot +import pygsti.tools.sdptools as sdps from pygsti.modelmembers.operations.lindbladcoefficients import LindbladCoefficientBlock from pygsti.modelpacks.legacy import std2Q_XXYYII from ..util import BaseCase, needs_cvxpy @@ -392,13 +394,45 @@ def test_jtrace_distance(self): val = ot.jtracedist(self.A_TP, self.B_unitary, mx_basis="pp") self.assertGreaterEqual(val, 0.5) + def validate_primal_diamondnorm(self, objective_val, operator, rho0, rho1, places=4): + d = rho0.shape[0] + assert rho0.shape == (d, d) + assert rho1.shape == (d, d) + assert operator.shape == (d**2, d**2) + left = np.kron(np.eye(d), la.sqrtm(rho0)) + right = np.kron(np.eye(d), la.sqrtm(rho1)) + Jop = ot._jam.jamiolkowski_iso(operator, 'pp', 'std', normalized=False) + arg = left @ Jop @ right + s = la.svdvals(arg) + expected_val = np.sum(s) + self.assertAlmostEqual(objective_val, expected_val, places=places) + return + @needs_cvxpy def test_diamond_distance(self): if SKIP_DIAMONDIST_ON_WIN and sys.platform.startswith('win'): return val = ot.diamonddist(self.A_TP, self.A_TP, mx_basis="pp") self.assertAlmostEqual(val, 0.0) - val = ot.diamonddist(self.A_TP, self.B_unitary, mx_basis="pp") - self.assertGreaterEqual(val, 0.7) + objective_val, modelvars = ot.diamonddist(self.A_TP, self.B_unitary, mx_basis="pp", return_all_vars=True) + rho0, rho1 = modelvars[1:] + self.validate_primal_diamondnorm(objective_val, self.A_TP - self.B_unitary, rho0, rho1) + self.assertGreaterEqual(objective_val, 0.7) + return + + @needs_cvxpy + def test_diamond_norm_epigraph(self): + if SKIP_DIAMONDIST_ON_WIN and sys.platform.startswith('win'): return + val0 = ot.diamonddist(self.A_TP, self.B_unitary, mx_basis="pp") + delta0 = self.A_TP - self.B_unitary + delta1 = bt.change_basis(delta0, 'pp', 'std') + import cvxpy as cp + obj_expr, constraints = sdps.diamond_norm_canon(delta1, 'std') + objective = cp.Minimize(obj_expr) + problem = cp.Problem(objective, constraints) + val1 = problem.solve(verbose=True, solver='MOSEK') + self.assertGreaterEqual(val0, 0.7) + self.assertAlmostEqual(val0, val1, places=4) + return def test_entanglement_fidelity(self): fidelity_TP_unitary= ot.entanglement_fidelity(self.A_TP, self.B_unitary, is_tp=True, is_unitary=True) @@ -432,9 +466,3 @@ def test_fidelity_upper_bound(self): bad_superoperator = ot.unitary_to_superop(Q) upperBound, _ = ot.fidelity_upper_bound(bad_superoperator) self.assertAlmostEqual(upperBound, 0.75) - np.random.seed(0) - Q = np.linalg.qr(np.random.randn(4,4) + 1j*np.random.randn(4,4))[0] - Q[:, 0] = 0.0 # zero out the first column - bad_superoperator = ot.unitary_to_superop(Q) - upperBound, _ = ot.fidelity_upper_bound(bad_superoperator) - self.assertAlmostEqual(upperBound, 0.75) From 9db6fe3dfb555b9078d82349a20583eaaf1f9056 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Fri, 27 Sep 2024 17:50:26 -0400 Subject: [PATCH 30/51] fix bugs in changes from last commit --- pygsti/tools/jamiolkowski.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/pygsti/tools/jamiolkowski.py b/pygsti/tools/jamiolkowski.py index 0c9e95249..485c399c0 100644 --- a/pygsti/tools/jamiolkowski.py +++ b/pygsti/tools/jamiolkowski.py @@ -125,20 +125,19 @@ def jamiolkowski_iso(operation_mx, op_mx_basis='pp', choi_mx_basis='pp', normali M = len(BVec) # can be < N if basis has multiple block dims assert(M == N), 'Expected {}, got {}'.format(M, N) - opMxInStdBasis_vec = opMxInStdBasis.reshape((N*N, 1)) - # ^ An explicit column vector. - choiMx_cols = [] + opMxInStdBasis_vec = opMxInStdBasis.ravel() + choiMx_rows = [] for i in range(M): rows = [] for j in range(M): BiBj = _np.kron(BVec[i], _np.conjugate(BVec[j])) BiBj /= _np.linalg.norm(BiBj) ** 2 rows.append(BiBj.conj().ravel()) - choiMx_cols.append( _np.array(rows) @ opMxInStdBasis_vec ) + choiMx_rows.append( _np.array(rows) @ opMxInStdBasis_vec ) if is_cvxpy_expression: - choiMx = cp.hstack(choiMx_cols) + choiMx = cp.vstack(choiMx_rows) else: - choiMx = _np.hstack(choiMx_cols) + choiMx = _np.vstack(choiMx_rows) # This construction results in a Jmx with trace == dim(H) = sqrt(operation_mx.shape[0]) # (dimension of density matrix) but we'd like a Jmx with trace == 1, so normalize: if normalized: @@ -302,7 +301,8 @@ def fast_jamiolkowski_iso_std_inv(choi_mx, op_mx_basis, normalized=True): assert(N * N == N2) # make sure N2 is a perfect square opMxInStdBasis = choi_mx.reshape((N, N, N, N)) if normalized: - opMxInStdBasis *= N + opMxInStdBasis = N * opMxInStdBasis + # ^ Don't do this in-place. opMxInStdBasis = _np.swapaxes(opMxInStdBasis, 1, 2).ravel() opMxInStdBasis = opMxInStdBasis.reshape((N2, N2)) op_mx_basis = _bt.create_basis_for_matrix(opMxInStdBasis, op_mx_basis) From e35167c48c6b9472cb4f5720401c221daa15e260 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Sat, 28 Sep 2024 09:59:14 -0400 Subject: [PATCH 31/51] add ability to return labels from leading_dxd_submatrix_basis_vectors --- pygsti/tools/optools.py | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/pygsti/tools/optools.py b/pygsti/tools/optools.py index 14c0ded58..6a151848c 100644 --- a/pygsti/tools/optools.py +++ b/pygsti/tools/optools.py @@ -531,7 +531,7 @@ def leaky_jtracedist(op_a, op_b, mx_basis, n_leak=0): return j_dist -def leading_dxd_submatrix_basis_vectors(d: int, n: int, current_basis): +def leading_dxd_submatrix_basis_vectors(d: int, n: int, current_basis, return_labels=False): """ Let "H" denote n^2 dimensional Hilbert-Schdmit space, and let "U" denote the d^2 dimensional subspace of H spanned by vectors whose Hermitian matrix representations @@ -604,12 +604,18 @@ def leading_dxd_submatrix_basis_vectors(d: int, n: int, current_basis): std_basis = _pgb.BuiltinBasis(name='std', dim_or_statespace=n**2) label2ind = {ell: idx for idx,ell in enumerate(std_basis.labels)} basis_ind = [] + basis_labels = [] for i in range(d): for j in range(d): - basis_ind.append(label2ind[f"({i},{j})"]) + ell = f"({i},{j})" + basis_ind.append(label2ind[ell]) + basis_labels.append(ell) basis_ind = _np.array(basis_ind) submatrix_basis_vectors = std_to_current[:, basis_ind] - return submatrix_basis_vectors + if return_labels: + return submatrix_basis_vectors, basis_labels + if not return_labels: + return submatrix_basis_vectors def average_gate_fidelity(a, b, mx_basis='pp', is_tp=None, is_unitary=None): From 611fc26404992b3fea17f3b9b0f46ec76d31173d Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Sat, 28 Sep 2024 10:04:19 -0400 Subject: [PATCH 32/51] define new OuterProductBasis class. Unclear if we'll want to keep this.' --- pygsti/baseobjs/basis.py | 351 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 351 insertions(+) diff --git a/pygsti/baseobjs/basis.py b/pygsti/baseobjs/basis.py index bd2bdf7de..87deebcbb 100644 --- a/pygsti/baseobjs/basis.py +++ b/pygsti/baseobjs/basis.py @@ -1760,3 +1760,354 @@ def create_simple_equivalent(self, builtin_basis_name=None): if all([c.name == first_comp_name for c in self.component_bases]): builtin_basis_name = first_comp_name # if all components have the same name return BuiltinBasis(builtin_basis_name, self.elsize, sparse=self.sparse) + +class OuterProdBasis(TensorProdBasis): + + def __init__(self, component_bases, name=None, longname=None, squeeze=True): + """ + Let ⨂ denote an infix operator where (a ⨂ b) := numpy.tensordot(a, b, axes=0). + + Suppose component_bases has length k. Under the default setting of squeeze=True, + the elements of this OuterProdBasis are all arrays of the form + + v1.squeeze() ⨂ v2.squeeze() ⨂ ... ⨂ vk.squeeze(), + + where vi belongs to component_bases[i].elements. The definition is changed in the + natural way if squeeze=False. + """ + TensorProdBasis.__init__(self, component_bases, name, longname) + if self.sparse: + raise NotImplementedError() + cnames = [c.name for c in self.component_bases] + if name is None: + self.name = "⨂".join(cnames) + if longname is None: + self.longname = "Outer-product basis with components " + ", ".join(cnames) + self._squeeze = squeeze + return + + @property + def elshape(self): + shape = [] + if self._squeeze: + for c in self.component_bases: + shape.extend(d for d in c.elshape if d > 1) + else: + for c in self.component_bases: + shape.extend(c.elshape) + return tuple(shape) + + def _lazy_build_elements(self): + compMxs = _np.zeros((self.size,) + self.elshape, 'complex') + comp_els = [c.elements for c in self.component_bases] + for i, factors in enumerate(_itertools.product(*comp_els)): + M = factors[0].squeeze() if self._squeeze else factors[0] + for f in factors[1:]: + if self._squeeze: + f = f.squeeze() + M = _np.tensordot(M, f, axes=0) + compMxs[i] = M + self._elements = compMxs + + def _lazy_build_labels(self): + self._labels = [] + comp_lbls = [c.labels for c in self.component_bases] + for i, factor_lbls in enumerate(_itertools.product(*comp_lbls)): + self._labels.append('⨂'.join(factor_lbls)) + + def is_equivalent(self, other, sparseness_must_match=True): + if not sparseness_must_match: + raise NotImplementedError() + otherIsBasis = isinstance(other, OuterProdBasis) + if not otherIsBasis: return False # can't be equal to a non-DirectSumBasis + return all([c1.is_equivalent(c2, True) + for (c1, c2) in zip(self.component_bases, other.component_bases)]) + + def _copy_with_toggled_sparsity(self): + raise NotImplementedError() + + def create_equivalent(self, builtin_basis_name): + raise NotImplementedError() + + def create_simple_equivalent(self, builtin_basis_name=None): + raise NotImplementedError() + + +class EmbeddedBasis(LazyBasis): + """ + A basis that embeds a basis for a smaller state space within a larger state space. + + The elements of an EmbeddedBasis are therefore just embedded versions + of the elements of the basis that is embedded. + + Parameters + ---------- + basis_to_embed : Basis + The basis being embedded. + + state_space_labels : StateSpaceLabels + An object describing the struture of the entire state space. + + target_labels : list or tuple + The labels contained in `stateSpaceLabels` which demarcate the + portions of the state space acted on by `basis_to_embed`. + + name : str, optional + The name of this basis. If `None`, the names of `basis_to_embed` + is joined with ':' characters to the elements of `target_labels`. + + longname : str, optional + A longer description of this basis. If `None`, then a long name is + automatically generated. + """ + + @classmethod + def embed_label(cls, lbl, target_labels): + """ + Gets the EmbeddedBasis label for `lbl`. + + Convenience method that gives the EmbeddedBasis label for `lbl` + without needing to construct the `EmbeddedBasis`. E.g. `"XX:1,2"`. + + Parameters + ---------- + lbl : str + Un-embedded basis element label, e.g. `"XX"`. + + target_labels : tuple + The target state space labels upon which this basis element + will be embedded, e.g. `(1,2)` + + Returns + ------- + str + The embedded-basis-element label as an EmbeddedBasis would + assign it. E.g. `"XX:1,2"`. + """ + return "%s:%s" % (lbl, ",".join(map(str, target_labels))) + + @classmethod + def unembed_label(cls, lbl, target_labels): + """ + Convenience method that performs the reverse of :meth:`embed_label` + + Parameters + ---------- + lbl : str + Embedded basis element label, e.g. `"XX:1,2"`. + + target_labels : tuple + The target state space labels upon which this basis element + will be embedded, e.g. `(1,2)` + + Returns + ------- + str + The un-embedded label, e.g. `"XX"`. + """ + suffix = ":" + ",".join(map(str, target_labels)) + if lbl.endswith(suffix): + return lbl[:-len(suffix)] + else: + raise ValueError("Cannot unembed '%s' - doesn't end in '%s'!" % (lbl, suffix)) + + def __init__(self, basis_to_embed, state_space, target_labels, name=None, longname=None): + ''' + Create a new EmbeddedBasis. + + Parameters + ---------- + basis_to_embed : Basis + The basis being embedded. + + state_space : StateSpace + An object describing the struture of the entire state space. + + target_labels : list or tuple + The labels contained in `stateSpaceLabels` which demarcate the + portions of the state space acted on by `basis_to_embed`. + + name : str, optional + The name of this basis. If `None`, the names of `basis_to_embed` + is joined with ':' characters to the elements of `target_labels`. + + longname : str, optional + A longer description of this basis. If `None`, then a long name is + automatically generated. + ''' + from pygsti.baseobjs.statespace import StateSpace as _StateSpace + self.embedded_basis = basis_to_embed + self.target_labels = target_labels + self.state_space = _StateSpace.cast(state_space) + + if name is None: + name = ':'.join((basis_to_embed.name,) + tuple(map(str, target_labels))) + if longname is None: + longname = "Embedded %s basis as %s within %s" % \ + (basis_to_embed.name, ':'.join(map(str, target_labels)), str(self.state_space)) + + real = basis_to_embed.real + sparse = basis_to_embed.sparse + + super(EmbeddedBasis, self).__init__(name, longname, real, sparse) + + def _to_nice_serialization(self): + state = super()._to_nice_serialization() + state.update({'name': self.name, + 'longname': self.longname, + 'state_space': self.state_space.to_nice_serialization(), + 'embedded_basis': self.embedded_basis.to_nice_serialization() + }) + return state + + @classmethod + def _from_nice_serialization(cls, state): + basis_to_embed = Basis.from_nice_serialization(state['embedded_basis']) + state_space = _StateSpace.from_nice_serialization(state['state_space']) + return cls(basis_to_embed, state_space, state['target_labels'], state['name'], state['longname']) + + @property + def dim(self): + """ + The dimension of the vector space this basis fully or partially + spans. Equivalently, the length of the `vector_elements` of the + basis. + """ + return self.state_space.dim + + @property + def size(self): + """ + The number of elements (or vector-elements) in the basis. + """ + return self.embedded_basis.size + + @property + def elshape(self): + """ + The shape of each element. Typically either a length-1 or length-2 + tuple, corresponding to vector or matrix elements, respectively. + Note that *vector elements* always have shape `(dim,)` (or `(dim,1)` + in the sparse case). + """ + elndim = self.embedded_basis.elndim + if elndim == 2: # a "matrix" basis + d = int(_np.sqrt(self.dim)) + assert(d**2 == self.dim), \ + "Dimension of state_space must be a perfect square when embedding a matrix basis" + elshape = (d, d) + elif elndim == 1: + elshape = (self.dim,) + else: + raise ValueError("Can only embed bases with .elndim == 1 or 2 (received %d)!" % elndim) + return elshape + + def __hash__(self): + return hash(tuple(hash(self.embedded_basis), self.target_labels, self.state_space)) + + def _lazy_build_elements(self): + """ Take a dense or sparse basis matrix and embed it. """ + #LAZY building of elements (in case we never need them) + if self.elndim == 2: # then use EmbeddedOp to do matrix + from ..modelmembers.operations import StaticArbitraryOp + from ..modelmembers.operations import EmbeddedOp + sslbls = self.state_space.copy() + sslbls.reduce_dims_densitymx_to_state_inplace() # because we're working with basis matrices not gates + + if self.sparse: + self._elements = [] + for spmx in self.embedded_basis.elements: + mxAsOp = StaticArbitraryOp(spmx.to_dense(), evotype='statevec') + self._elements.append(EmbeddedOp(sslbls, self.target_labels, + mxAsOp).to_sparse()) + else: + self._elements = _np.zeros((self.size,) + self.elshape, 'complex') + for i, mx in enumerate(self.embedded_basis.elements): + self._elements[i] = EmbeddedOp( + sslbls, self.target_labels, StaticArbitraryOp(mx, evotype='statevec') + ).to_dense(on_space='HilbertSchmidt') + else: + # we need to perform embedding using vectors rather than matrices - doable, but + # not needed yet, so defer implementation to later. + raise NotImplementedError("Embedding *vector*-type bases not implemented yet") + + def _lazy_build_labels(self): + self._labels = [EmbeddedBasis.embed_label(lbl, self.target_labels) + for lbl in self.embedded_basis.labels] + + def _copy_with_toggled_sparsity(self): + return EmbeddedBasis(self.embedded_basis._copy_with_toggled_sparsity(), + self.state_space, + self.target_labels, + self.name, self.longname) + + def is_equivalent(self, other, sparseness_must_match=True): + """ + Tests whether this basis is equal to another basis, optionally ignoring sparseness. + + Parameters + ----------- + other : Basis or str + The basis to compare with. + + sparseness_must_match : bool, optional + If `False` then comparison ignores differing sparseness, and this function + returns `True` when the two bases are equal except for their `.sparse` values. + + Returns + ------- + bool + """ + otherIsBasis = isinstance(other, EmbeddedBasis) + if not otherIsBasis: return False # can't be equal to a non-EmbeddedBasis + if self.target_labels != other.target_labels or self.state_space != other.state_space: + return False + return self.embedded_basis.is_equivalent(other.embedded_basis, sparseness_must_match) + + def create_equivalent(self, builtin_basis_name): + """ + Create an equivalent basis with components of type `builtin_basis_name`. + + Create a Basis that is equivalent in structure & dimension to this + basis but whose simple components (perhaps just this basis itself) is + of the builtin basis type given by `builtin_basis_name`. + + Parameters + ---------- + builtin_basis_name : str + The name of a builtin basis, e.g. `"pp"`, `"gm"`, or `"std"`. Used to + construct the simple components of the returned basis. + + Returns + ------- + EmbeddedBasis + """ + equiv_embedded = self.embedded_basis.create_equivalent(builtin_basis_name) + return EmbeddedBasis(equiv_embedded, self.state_space, self.target_labels) + + def create_simple_equivalent(self, builtin_basis_name=None): + """ + Create a basis of type `builtin_basis_name` whose elements are compatible with this basis. + + Create a simple basis *and* one without components (e.g. a + :class:`TensorProdBasis`, is a simple basis w/components) of the + builtin type specified whose dimension is compatible with the + *elements* of this basis. This function might also be named + "element_equivalent", as it returns the `builtin_basis_name`-analogue + of the standard basis that this basis's elements are expressed in. + + Parameters + ---------- + builtin_basis_name : str, optional + The name of the built-in basis to use. If `None`, then a + copy of this basis is returned (if it's simple) or this + basis's name is used to try to construct a simple and + component-free version of the same builtin-basis type. + + Returns + ------- + Basis + """ + if builtin_basis_name is None: + builtin_basis_name = self.embedded_basis.name # default + return BuiltinBasis(builtin_basis_name, self.elsize, sparse=self.sparse) From 26eb4b2d8e33a3488ecc8ea6f68e3c99b304ada9 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Sat, 28 Sep 2024 10:59:24 -0400 Subject: [PATCH 33/51] rollback change in 1c7d5aaf9d2e45888132c1d6995e3d868a039ba2 --- pygsti/baseobjs/basisconstructors.py | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/pygsti/baseobjs/basisconstructors.py b/pygsti/baseobjs/basisconstructors.py index 0706ef6bd..5fa8f3a19 100644 --- a/pygsti/baseobjs/basisconstructors.py +++ b/pygsti/baseobjs/basisconstructors.py @@ -862,12 +862,8 @@ def _is_integer(x): nQubits = _np.log2(matrix_dim) if not _is_integer(nQubits): - nQutrits = _np.log(matrix_dim) / _np.log(3) - if _is_integer(nQutrits): - return gm_matrices(matrix_dim) - else: - raise ValueError( - "Dimension for Pauli tensor product matrices must be an integer *power of 2* (not %d)" % matrix_dim) + raise ValueError( + "Dimension for Pauli tensor product matrices must be an integer *power of 2* (not %d)" % matrix_dim) nQubits = int(round(nQubits)) if nQubits == 0: # special case: return single 1x1 identity mx From a87177bd4b363959c7b3a91bc78526f471dde81a Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Sat, 28 Sep 2024 11:24:08 -0400 Subject: [PATCH 34/51] remove explicit call to MOSEK --- test/unit/tools/test_optools.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/unit/tools/test_optools.py b/test/unit/tools/test_optools.py index a94a4c280..17c7641a9 100644 --- a/test/unit/tools/test_optools.py +++ b/test/unit/tools/test_optools.py @@ -429,7 +429,7 @@ def test_diamond_norm_epigraph(self): obj_expr, constraints = sdps.diamond_norm_canon(delta1, 'std') objective = cp.Minimize(obj_expr) problem = cp.Problem(objective, constraints) - val1 = problem.solve(verbose=True, solver='MOSEK') + val1 = problem.solve(verbose=True, solver='CLARABEL') self.assertGreaterEqual(val0, 0.7) self.assertAlmostEqual(val0, val1, places=4) return From d386741f5c220c53196a4dab9200adad17214993 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Mon, 18 Nov 2024 10:11:42 -0500 Subject: [PATCH 35/51] update add_gauge_opt so it can accept a dict specifying a gaugeopt suite rather than a literal gaugeoptsuite object. --- pygsti/protocols/gst.py | 25 ++++++++++++++----------- 1 file changed, 14 insertions(+), 11 deletions(-) diff --git a/pygsti/protocols/gst.py b/pygsti/protocols/gst.py index 17323bc3f..61cae7445 100644 --- a/pygsti/protocols/gst.py +++ b/pygsti/protocols/gst.py @@ -872,14 +872,13 @@ def cast(cls, obj): def __init__(self, gaugeopt_suite_names=None, gaugeopt_argument_dicts=None, gaugeopt_target=None): super().__init__() - if gaugeopt_suite_names is not None: - if gaugeopt_suite_names == 'none': - self.gaugeopt_suite_names = None - else: - self.gaugeopt_suite_names = (gaugeopt_suite_names,) \ - if isinstance(gaugeopt_suite_names, str) else tuple(gaugeopt_suite_names) - else: + if gaugeopt_suite_names is None or gaugeopt_suite_names == 'none': self.gaugeopt_suite_names = None + elif isinstance(gaugeopt_suite_names, str): + self.gaugeopt_suite_names = gaugeopt_suite_names + else: + self.gaugeopt_suite_names = tuple(gaugeopt_suite_names) + if gaugeopt_argument_dicts is not None: self.gaugeopt_argument_dicts = gaugeopt_argument_dicts.copy() @@ -2039,8 +2038,8 @@ def _add_gaugeopt_and_badfit(results, estlbl, target_model, gaugeopt_suite, return results -def _add_gauge_opt(results, base_est_label, gaugeopt_suite, starting_model, - unreliable_ops, comm=None, verbosity=0): +def _add_gauge_opt(results, base_est_label, gaugeopt_suite_or_gsdict, starting_model, + unreliable_ops=(), comm=None, verbosity=0): """ Add a gauge optimization to an estimate. @@ -2085,8 +2084,12 @@ def _add_gauge_opt(results, base_est_label, gaugeopt_suite, starting_model, printer = _baseobjs.VerbosityPrinter.create_printer(verbosity, comm) #Get gauge optimization dictionary - gaugeopt_suite_dict = gaugeopt_suite.to_dictionary(starting_model, - unreliable_ops, printer - 1) + if isinstance(gaugeopt_suite_or_gsdict, dict): + gaugeopt_suite_dict = gaugeopt_suite_or_gsdict + else: + gaugeopt_suite_dict = gaugeopt_suite_or_gsdict.to_dictionary( + starting_model, unreliable_ops, printer - 1 + ) #Gauge optimize to list of gauge optimization parameters for go_label, goparams in gaugeopt_suite_dict.items(): From 504827733d809b809fd59b75c83327954e3288e9 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Mon, 18 Nov 2024 10:15:24 -0500 Subject: [PATCH 36/51] REVERT ME --- pygsti/circuits/circuit.py | 22 ++++++++++++---------- pygsti/data/dataset.py | 2 +- pygsti/objectivefns/objectivefns.py | 2 +- 3 files changed, 14 insertions(+), 12 deletions(-) diff --git a/pygsti/circuits/circuit.py b/pygsti/circuits/circuit.py index 7624d1f98..f13b0bf32 100644 --- a/pygsti/circuits/circuit.py +++ b/pygsti/circuits/circuit.py @@ -801,12 +801,13 @@ def str(self, value): self._str = value def __hash__(self): - if not self._static: - _warnings.warn(("Editable circuit is being converted to read-only" - " mode in order to hash it. You should call" - " circuit.done_editing() beforehand.")) - self.done_editing() - return self._hash + return hash(self.tup) + # if not self._static: + # _warnings.warn(("Editable circuit is being converted to read-only" + # " mode in order to hash it. You should call" + # " circuit.done_editing() beforehand.")) + # self.done_editing() + # return self._hash def __len__(self): return len(self._labels) @@ -978,10 +979,11 @@ def __eq__(self, x): if len(self) != len(x): return False else: - if self._static and x._static: - return self._hash == x._hash - else: - return self.tup == x.tup + return self.tup == x.tup + # if self._static and x._static: + # return self._hash == x._hash + # else: + # return self.tup == x.tup elif x is None: return False else: diff --git a/pygsti/data/dataset.py b/pygsti/data/dataset.py index ce7bb52c6..02ef47582 100644 --- a/pygsti/data/dataset.py +++ b/pygsti/data/dataset.py @@ -1176,7 +1176,7 @@ def _get_row(self, circuit): # needed because name-only Labels don't hash the same as strings # so key lookups need to be done at least with tuples of Labels. circuit = _cir.Circuit.cast(circuit) - + circuit._static = False #Note: cirIndex value is either an int (non-static) or a slice (static) cirIndex = self.cirIndex[circuit] repData = self.repData[cirIndex] if (self.repData is not None) else None diff --git a/pygsti/objectivefns/objectivefns.py b/pygsti/objectivefns/objectivefns.py index 0c44fe74c..97b6cf4e2 100644 --- a/pygsti/objectivefns/objectivefns.py +++ b/pygsti/objectivefns/objectivefns.py @@ -871,7 +871,7 @@ def __init__(self, model, dataset, circuits=None, resource_alloc=None, array_typ #If a matrix layout then we have some precached circuit structures we can #grab to speed up store generation. - if isinstance(self.layout, _MatrixCOPALayout): + if False: #Grab the split_circuit_cache and down select to those in #self.circuits self.split_circuit_cache = self.layout.split_circuit_cache From 328a5bce70656d6556de6902cd9f75a2ae2f6432 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Mon, 18 Nov 2024 10:17:59 -0500 Subject: [PATCH 37/51] comments pointing to possible bugs --- pygsti/models/gaugegroup.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/pygsti/models/gaugegroup.py b/pygsti/models/gaugegroup.py index ca644a72b..11b2ee232 100644 --- a/pygsti/models/gaugegroup.py +++ b/pygsti/models/gaugegroup.py @@ -669,6 +669,8 @@ def __init__(self, state_space, evotype='default'): rtrans = _np.identity(dim, 'd') baseMx = _np.identity(dim, 'd') parameterArray = _np.zeros(dim, 'd') + # ^ ... shouldn't that be ones? + # parameterArray = _np.ones(dim, 'd') parameterToBaseIndicesMap = {i: [(i, i)] for i in range(dim)} operation = _op.LinearlyParamArbitraryOp(baseMx, parameterArray, parameterToBaseIndicesMap, ltrans, rtrans, real=True, evotype=evotype, state_space=state_space) @@ -725,6 +727,7 @@ def __init__(self, state_space, evotype='default'): rtrans = _np.identity(dim, 'd') baseMx = _np.identity(dim, 'd') parameterArray = _np.zeros(dim - 1, 'd') + # parameterArray = _np.ones(dim - 1, 'd') parameterToBaseIndicesMap = {i: [(i + 1, i + 1)] for i in range(dim - 1)} operation = _op.LinearlyParamArbitraryOp(baseMx, parameterArray, parameterToBaseIndicesMap, ltrans, rtrans, real=True, evotype=evotype, state_space=state_space) @@ -835,6 +838,7 @@ def __init__(self, state_space, evotype='default'): rtrans = _np.identity(dim, 'd') baseMx = _np.identity(dim, 'd') parameterArray = _np.zeros(2, 'd') + # parameterArray = _np.ones(2, 'd') parameterToBaseIndicesMap = {0: [(0, 0)], 1: [(i, i) for i in range(1, dim)]} operation = _op.LinearlyParamArbitraryOp(baseMx, parameterArray, parameterToBaseIndicesMap, ltrans, rtrans, @@ -893,6 +897,7 @@ def __init__(self, state_space, evotype='default'): rtrans = _np.identity(dim, 'd') baseMx = _np.identity(dim, 'd') parameterArray = _np.zeros(1, 'd') + # parameterArray = _np.ones(1, 'd') parameterToBaseIndicesMap = {0: [(i, i) for i in range(1, dim)]} operation = _op.LinearlyParamArbitraryOp(baseMx, parameterArray, parameterToBaseIndicesMap, ltrans, rtrans, real=True, evotype=evotype, state_space=state_space) From 018d5e848e38961d7a416d78d80ab20bd7908106 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Mon, 18 Nov 2024 10:23:06 -0500 Subject: [PATCH 38/51] note --- pygsti/report/reportables.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/pygsti/report/reportables.py b/pygsti/report/reportables.py index d4d501187..b4c5eeed6 100644 --- a/pygsti/report/reportables.py +++ b/pygsti/report/reportables.py @@ -2144,6 +2144,11 @@ def general_decomposition(model_a, model_b): target_evals = _np.linalg.eigvals(targetOp) if _np.any(_np.isclose(target_evals, -1.0)): + # It's not clear if this codepath is required for correctness. + # In my experience it can lead to runtime errors when basis conversions + # end up with a matrix with imaginary part when the matrix should only + # be real. So I'd like to wrap this in a try-catch. + # target_logG = _tools.unitary_superoperator_matrix_log(targetOp, mxBasis) logG = _tools.approximate_matrix_log(gate, target_logG) else: From 34fc1dc364f09e8dc223af45e0dbf9a133f13631 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Mon, 18 Nov 2024 14:52:56 -0500 Subject: [PATCH 39/51] leave note for myself --- pygsti/modelmembers/operations/fulltpop.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/pygsti/modelmembers/operations/fulltpop.py b/pygsti/modelmembers/operations/fulltpop.py index b8ce2251e..58722be55 100644 --- a/pygsti/modelmembers/operations/fulltpop.py +++ b/pygsti/modelmembers/operations/fulltpop.py @@ -48,6 +48,12 @@ class FullTPOp(_DenseOperator, _Torchable): of this state as a super-operator. If None, certain functionality, such as access to Kraus operators, will be unavailable. + NOTE: this function might raise a weird error message if the passed + mx isn't in the pp or gm bases. I noticed this when using an + ExplicitBasis that is "leakage-friendly." + The implicit underlying assumption is that the first element + of this basis is the identity matrix. + evotype : Evotype or str, optional The evolution type. The special value `"default"` is equivalent to specifying the value of `pygsti.evotypes.Evotype.default_evotype`. From d4fcaee3d3c4dac8e2e2ea5fae282f14194ffd83 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Thu, 21 Nov 2024 07:17:09 -0500 Subject: [PATCH 40/51] reduce number of options available in non-LS gauge optimization --- pygsti/algorithms/gaugeopt.py | 47 ++++++++++++++++++++++++++++------- 1 file changed, 38 insertions(+), 9 deletions(-) diff --git a/pygsti/algorithms/gaugeopt.py b/pygsti/algorithms/gaugeopt.py index d60f75c76..5993c152a 100644 --- a/pygsti/algorithms/gaugeopt.py +++ b/pygsti/algorithms/gaugeopt.py @@ -587,9 +587,23 @@ def _mock_objective_fn(v): else: # non-least-squares case where objective function returns a single float # and (currently) there's no analytic jacobian + dim = int(_np.sqrt(mxBasis.dim)) if n_leak > 0: - dim = int(_np.sqrt(mxBasis.dim)) B = _tools.leading_dxd_submatrix_basis_vectors(dim - n_leak, dim, mxBasis) + """ + ^ Need to do something else. + + In a leakage-friendly basis the operation matrices can be partitioned as + [comp , semileak1] + [semileak2, fulleak ]. + Instead of projecting only on to comp, we want to keep everything + except fullleak. ... But I don't think we can implement that just by + pre-multiplying by a projector. + + I think this distinction might not be significant under certain idealized + assumptions, but small deviations from those conditions (which are certain + to happen in practice) might make it matter. + """ P = B @ B.T.conj() if _np.linalg.norm(P.imag) > 1e-12: raise ValueError() @@ -601,10 +615,11 @@ def _mock_objective_fn(v): # the second element of the tuple, but those aren't implemented yet so # for now I'm setting the second entry to None. -- Riley else: - transform_mx_arg = None + transform_mx_arg = (_np.eye(mxBasis.dim), None) assert gates_metric != "frobeniustt" assert spam_metric != "frobeniustt" + # assert spam_metric == gates_metric # ^ Erik and Corey said these are rarely used. I've removed support for # them in this codepath (non-LS optimizer) in order to make it easier to # read my updated code for leakage-aware metrics. It wouldn't be hard to @@ -625,18 +640,32 @@ def _objective_fn(gauge_group_el, oob_check): ret += _np.sum(spamPenaltyVec) if target_model is not None: - # Leakage-aware metric supported, per implementation in mdl.frobeniusdist. - # Refer to how mdl.frobeniusdist handles the case when transform_mx_arg - # is a tuple in order to understand how the leakage-aware metric is defined. + """ + Leakage-aware metric supported, per implementation in mdl.frobeniusdist. + Refer to how mdl.frobeniusdist handles the case when transform_mx_arg + is a tuple in order to understand how the leakage-aware metric is defined. + + Idea: raise an error if mxBasis isn't leakage-friendly. + PROBLEM: if we're deep in the code then we end up needing to be + working in a basis that has the identity matrix as its + first element. The leakage-friendly basis doesn't even + have the identity matrix as an element, let alone the first element + + TODO: dig into function calls below and see where we can + access the mxBasis object. (I'm sure we can.) + Looks like it isn't accessible within ExplicitOpModelCalc objects. + It's most definitely available in ExplicitOpModel. + """ if "frobenius" in gates_metric: if spam_metric == gates_metric: val = mdl.frobeniusdist(target_model, transform_mx_arg, item_weights) else: - wts = item_weights.copy(); wts['spam'] = 0.0 + wts = item_weights.copy() + wts['spam'] = 0.0 for k in wts: - if k in mdl.preps or \ - k in mdl.povms: wts[k] = 0.0 - val = mdl.frobeniusdist(target_model, transform_mx_arg, wts) + if k in mdl.preps or k in mdl.povms: + wts[k] = 0.0 + val = mdl.frobeniusdist(target_model, transform_mx_arg, wts, n_leak) if "squared" in gates_metric: val = val ** 2 ret += val From 229f3fc4fe82b4215a8b2967d2fb9113f6789d36 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Thu, 21 Nov 2024 07:22:27 -0500 Subject: [PATCH 41/51] last commit actually just included comments. Material changes are in this commit --- pygsti/algorithms/gaugeopt.py | 164 ++++++++++++++-------------------- 1 file changed, 68 insertions(+), 96 deletions(-) diff --git a/pygsti/algorithms/gaugeopt.py b/pygsti/algorithms/gaugeopt.py index 5993c152a..823d98c99 100644 --- a/pygsti/algorithms/gaugeopt.py +++ b/pygsti/algorithms/gaugeopt.py @@ -619,6 +619,8 @@ def _mock_objective_fn(v): assert gates_metric != "frobeniustt" assert spam_metric != "frobeniustt" + assert spam_metric == gates_metric + metric = spam_metric # assert spam_metric == gates_metric # ^ Erik and Corey said these are rarely used. I've removed support for # them in this codepath (non-LS optimizer) in order to make it easier to @@ -639,103 +641,73 @@ def _objective_fn(gauge_group_el, oob_check): spamPenaltyVec = _spam_penalty(mdl, spam_penalty_factor, mdl.basis) ret += _np.sum(spamPenaltyVec) - if target_model is not None: - """ - Leakage-aware metric supported, per implementation in mdl.frobeniusdist. - Refer to how mdl.frobeniusdist handles the case when transform_mx_arg - is a tuple in order to understand how the leakage-aware metric is defined. - - Idea: raise an error if mxBasis isn't leakage-friendly. - PROBLEM: if we're deep in the code then we end up needing to be - working in a basis that has the identity matrix as its - first element. The leakage-friendly basis doesn't even - have the identity matrix as an element, let alone the first element - - TODO: dig into function calls below and see where we can - access the mxBasis object. (I'm sure we can.) - Looks like it isn't accessible within ExplicitOpModelCalc objects. - It's most definitely available in ExplicitOpModel. - """ - if "frobenius" in gates_metric: - if spam_metric == gates_metric: - val = mdl.frobeniusdist(target_model, transform_mx_arg, item_weights) - else: - wts = item_weights.copy() - wts['spam'] = 0.0 - for k in wts: - if k in mdl.preps or k in mdl.povms: - wts[k] = 0.0 - val = mdl.frobeniusdist(target_model, transform_mx_arg, wts, n_leak) - if "squared" in gates_metric: - val = val ** 2 - ret += val - - elif gates_metric == "fidelity": - # Leakage-aware metric supported, using leaky_entanglement_fidelity. - for opLbl in mdl.operations: - wt = item_weights.get(opLbl, opWeight) - top = target_model.operations[opLbl].to_dense() - mop = mdl.operations[opLbl].to_dense() - ret += wt * (1.0 - _tools.leaky_entanglement_fidelity(top, mop, mxBasis, n_leak))**2 - - elif gates_metric == "tracedist": - # Leakage-aware metric supported, using leaky_jtracedist. - for opLbl in mdl.operations: - wt = item_weights.get(opLbl, opWeight) - top = target_model.operations[opLbl].to_dense() - mop = mdl.operations[opLbl].to_dense() - ret += wt * _tools.leaky_jtracedist(top, mop, mxBasis, n_leak) - - else: raise ValueError("Invalid gates_metric: %s" % gates_metric) - - if "frobenius" in spam_metric and gates_metric == spam_metric: - # We already handled SPAM error in this case. Just return. - return ret - - elif "frobenius" in spam_metric: - # Leakage-aware metric supported in principle via implementation in - # mdl.frobeniusdist (check implementation to see how it handles the - # case when transform_mx_arg is a tuple). - wts = item_weights.copy(); wts['gates'] = 0.0 - for k in wts: - if k in mdl.operations or \ - k in mdl.instruments: wts[k] = 0.0 - val = mdl.frobeniusdist(target_model, transform_mx_arg, wts) - if "squared" in spam_metric: - val = val ** 2 - ret += val - - elif spam_metric == "fidelity": - # Leakage-aware metrics NOT available - for preplabel, m_prep in mdl.preps.items(): - wt = item_weights.get(preplabel, spamWeight) - rhoMx1 = _tools.vec_to_stdmx(m_prep.to_dense(), mxBasis) - t_prep = target_model.preps[preplabel] - rhoMx2 = _tools.vec_to_stdmx(t_prep.to_dense(), mxBasis) - ret += wt * (1.0 - _tools.fidelity(rhoMx1, rhoMx2))**2 - - for povmlabel in mdl.povms.keys(): - wt = item_weights.get(povmlabel, spamWeight) - fidelity = _tools.povm_fidelity(mdl, target_model, povmlabel) - ret += wt * (1.0 - fidelity)**2 - - elif spam_metric == "tracedist": - # Leakage-aware metrics NOT available. - for preplabel, m_prep in mdl.preps.items(): - wt = item_weights.get(preplabel, spamWeight) - rhoMx1 = _tools.vec_to_stdmx(m_prep.to_dense(), mxBasis) - t_prep = target_model.preps[preplabel] - rhoMx2 = _tools.vec_to_stdmx(t_prep.to_dense(), mxBasis) - ret += wt * _tools.tracedist(rhoMx1, rhoMx2) - - for povmlabel in mdl.povms.keys(): - wt = item_weights.get(povmlabel, spamWeight) - ret += wt * _tools.povm_jtracedist(mdl, target_model, povmlabel) - - else: raise ValueError("Invalid spam_metric: %s" % spam_metric) - - return ret + assert target_model is not None + """ + Leakage-aware metric supported, per implementation in mdl.frobeniusdist. + Refer to how mdl.frobeniusdist handles the case when transform_mx_arg + is a tuple in order to understand how the leakage-aware metric is defined. + + Idea: raise an error if mxBasis isn't leakage-friendly. + PROBLEM: if we're deep in the code then we end up needing to be + working in a basis that has the identity matrix as its + first element. The leakage-friendly basis doesn't even + have the identity matrix as an element, let alone the first element + + TODO: dig into function calls below and see where we can + access the mxBasis object. (I'm sure we can.) + Looks like it isn't accessible within ExplicitOpModelCalc objects. + It's most definitely available in ExplicitOpModel. + """ + if "frobenius" in metric: + val = mdl.frobeniusdist(target_model, transform_mx_arg, item_weights) + if "squared" in metric: + val = val ** 2 + ret += val + return ret + + elif metric == "fidelity": + # Leakage-aware metric supported for gates, using leaky_entanglement_fidelity. + for opLbl in mdl.operations: + wt = item_weights.get(opLbl, opWeight) + top = target_model.operations[opLbl].to_dense() + mop = mdl.operations[opLbl].to_dense() + ret += wt * (1.0 - _tools.leaky_entanglement_fidelity(top, mop, mxBasis, n_leak))**2 + # Leakage-aware metrics NOT available for SPAM + for preplabel, m_prep in mdl.preps.items(): + wt = item_weights.get(preplabel, spamWeight) + rhoMx1 = _tools.vec_to_stdmx(m_prep.to_dense(), mxBasis) + t_prep = target_model.preps[preplabel] + rhoMx2 = _tools.vec_to_stdmx(t_prep.to_dense(), mxBasis) + ret += wt * (1.0 - _tools.fidelity(rhoMx1, rhoMx2))**2 + for povmlabel in mdl.povms.keys(): + wt = item_weights.get(povmlabel, spamWeight) + fidelity = _tools.povm_fidelity(mdl, target_model, povmlabel) + ret += wt * (1.0 - fidelity)**2 + return ret + + elif metric == "tracedist": + # Leakage-aware metric supported, using leaky_jtracedist. + for opLbl in mdl.operations: + wt = item_weights.get(opLbl, opWeight) + top = target_model.operations[opLbl].to_dense() + mop = mdl.operations[opLbl].to_dense() + ret += wt * _tools.leaky_jtracedist(top, mop, mxBasis, n_leak) + # Leakage-aware metrics NOT available for SPAM + for preplabel, m_prep in mdl.preps.items(): + wt = item_weights.get(preplabel, spamWeight) + rhoMx1 = _tools.vec_to_stdmx(m_prep.to_dense(), mxBasis) + t_prep = target_model.preps[preplabel] + rhoMx2 = _tools.vec_to_stdmx(t_prep.to_dense(), mxBasis) + ret += wt * _tools.tracedist(rhoMx1, rhoMx2) + for povmlabel in mdl.povms.keys(): + wt = item_weights.get(povmlabel, spamWeight) + ret += wt * _tools.povm_jtracedist(mdl, target_model, povmlabel) + return ret + else: + raise ValueError("Invalid metric: %s" % metric) + # end _objective_fn + _jacobian_fn = None return _objective_fn, _jacobian_fn From f8b816a21bcff561ed35f8ef993ec9c08200d02e Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Thu, 21 Nov 2024 08:01:45 -0500 Subject: [PATCH 42/51] simplify gaugeopt in the non-LS case --- pygsti/algorithms/gaugeopt.py | 91 +++++++++++++++-------------------- 1 file changed, 40 insertions(+), 51 deletions(-) diff --git a/pygsti/algorithms/gaugeopt.py b/pygsti/algorithms/gaugeopt.py index 823d98c99..51225a83a 100644 --- a/pygsti/algorithms/gaugeopt.py +++ b/pygsti/algorithms/gaugeopt.py @@ -587,45 +587,24 @@ def _mock_objective_fn(v): else: # non-least-squares case where objective function returns a single float # and (currently) there's no analytic jacobian - dim = int(_np.sqrt(mxBasis.dim)) - if n_leak > 0: - B = _tools.leading_dxd_submatrix_basis_vectors(dim - n_leak, dim, mxBasis) - """ - ^ Need to do something else. - - In a leakage-friendly basis the operation matrices can be partitioned as - [comp , semileak1] - [semileak2, fulleak ]. - Instead of projecting only on to comp, we want to keep everything - except fullleak. ... But I don't think we can implement that just by - pre-multiplying by a projector. - - I think this distinction might not be significant under certain idealized - assumptions, but small deviations from those conditions (which are certain - to happen in practice) might make it matter. - """ - P = B @ B.T.conj() - if _np.linalg.norm(P.imag) > 1e-12: - raise ValueError() - else: - P = P.real - transform_mx_arg = (P, None) - # ^ The semantics of this tuple are defined by the frobeniusdist function - # in the ExplicitOpModelCalc class. There are intended semantics for - # the second element of the tuple, but those aren't implemented yet so - # for now I'm setting the second entry to None. -- Riley - else: - transform_mx_arg = (_np.eye(mxBasis.dim), None) + assert target_model is not None assert gates_metric != "frobeniustt" assert spam_metric != "frobeniustt" - assert spam_metric == gates_metric - metric = spam_metric - # assert spam_metric == gates_metric # ^ Erik and Corey said these are rarely used. I've removed support for # them in this codepath (non-LS optimizer) in order to make it easier to # read my updated code for leakage-aware metrics. It wouldn't be hard to # add support back, but I just want to keep things simple. -- Riley + assert spam_metric == gates_metric + metric = spam_metric + + dim = int(_np.sqrt(mxBasis.dim)) + B = _tools.leading_dxd_submatrix_basis_vectors(dim - n_leak, dim, mxBasis) + P = B @ B.T.conj() + assert _np.linalg.norm(P.imag) <= 1e-12 + P = P.real + + log = [] def _objective_fn(gauge_group_el, oob_check): mdl = _transform_with_oob_check(model, gauge_group_el, oob_check) @@ -641,25 +620,35 @@ def _objective_fn(gauge_group_el, oob_check): spamPenaltyVec = _spam_penalty(mdl, spam_penalty_factor, mdl.basis) ret += _np.sum(spamPenaltyVec) - assert target_model is not None - """ - Leakage-aware metric supported, per implementation in mdl.frobeniusdist. - Refer to how mdl.frobeniusdist handles the case when transform_mx_arg - is a tuple in order to understand how the leakage-aware metric is defined. - - Idea: raise an error if mxBasis isn't leakage-friendly. - PROBLEM: if we're deep in the code then we end up needing to be - working in a basis that has the identity matrix as its - first element. The leakage-friendly basis doesn't even - have the identity matrix as an element, let alone the first element - - TODO: dig into function calls below and see where we can - access the mxBasis object. (I'm sure we can.) - Looks like it isn't accessible within ExplicitOpModelCalc objects. - It's most definitely available in ExplicitOpModel. - """ if "frobenius" in metric: - val = mdl.frobeniusdist(target_model, transform_mx_arg, item_weights) + d = 0 + nSummands = 0.0 + for opLabel, gate in mdl.operations.items(): + wt = item_weights.get(opLabel, opWeight) + gate_mx = gate.to_dense() + other_mx = target_model.operations[opLabel].to_dense() + delta = gate_mx - other_mx + delta = delta @ P + val = _np.linalg.norm(delta.flatten()) + d += wt * val**2 + nSummands += wt * (gate.dim)**2 + + for lbl, rhoV in mdl.preps.items(): + wt = item_weights.get(lbl, spamWeight) + d += wt * rhoV.frobeniusdist_squared(target_model.preps[lbl], None, None) + nSummands += wt * rhoV.dim + + for lbl, Evec in mdl.effects.items(): + wt = item_weights.get(lbl, spamWeight) + evec = Evec.to_dense() + other = target_model.effects[lbl].to_dense() + delta = evec - other + delta = delta @ P + val = _np.linalg.norm(delta.flatten()) + d += wt * val**2 + nSummands += wt * Evec.dim + + val = _np.sqrt(d / nSummands) if "squared" in metric: val = val ** 2 ret += val @@ -707,7 +696,7 @@ def _objective_fn(gauge_group_el, oob_check): else: raise ValueError("Invalid metric: %s" % metric) # end _objective_fn - + _jacobian_fn = None return _objective_fn, _jacobian_fn From 2d8dc7ef243c06ddfd914d6a8509888ec08fb03d Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Thu, 21 Nov 2024 08:07:21 -0500 Subject: [PATCH 43/51] remove non-Frobenius objectives from gauge optimization --- pygsti/algorithms/gaugeopt.py | 106 ++++++++++------------------------ 1 file changed, 32 insertions(+), 74 deletions(-) diff --git a/pygsti/algorithms/gaugeopt.py b/pygsti/algorithms/gaugeopt.py index 51225a83a..1b9e9dda3 100644 --- a/pygsti/algorithms/gaugeopt.py +++ b/pygsti/algorithms/gaugeopt.py @@ -619,82 +619,40 @@ def _objective_fn(gauge_group_el, oob_check): mdl.basis = mxBasis spamPenaltyVec = _spam_penalty(mdl, spam_penalty_factor, mdl.basis) ret += _np.sum(spamPenaltyVec) + + d = 0 + nSummands = 0.0 + for opLabel, gate in mdl.operations.items(): + wt = item_weights.get(opLabel, opWeight) + gate_mx = gate.to_dense() + other_mx = target_model.operations[opLabel].to_dense() + delta = gate_mx - other_mx + delta = delta @ P + val = _np.linalg.norm(delta.flatten()) + d += wt * val**2 + nSummands += wt * (gate.dim)**2 + + for lbl, rhoV in mdl.preps.items(): + wt = item_weights.get(lbl, spamWeight) + d += wt * rhoV.frobeniusdist_squared(target_model.preps[lbl], None, None) + nSummands += wt * rhoV.dim - if "frobenius" in metric: - d = 0 - nSummands = 0.0 - for opLabel, gate in mdl.operations.items(): - wt = item_weights.get(opLabel, opWeight) - gate_mx = gate.to_dense() - other_mx = target_model.operations[opLabel].to_dense() - delta = gate_mx - other_mx - delta = delta @ P - val = _np.linalg.norm(delta.flatten()) - d += wt * val**2 - nSummands += wt * (gate.dim)**2 - - for lbl, rhoV in mdl.preps.items(): - wt = item_weights.get(lbl, spamWeight) - d += wt * rhoV.frobeniusdist_squared(target_model.preps[lbl], None, None) - nSummands += wt * rhoV.dim - - for lbl, Evec in mdl.effects.items(): - wt = item_weights.get(lbl, spamWeight) - evec = Evec.to_dense() - other = target_model.effects[lbl].to_dense() - delta = evec - other - delta = delta @ P - val = _np.linalg.norm(delta.flatten()) - d += wt * val**2 - nSummands += wt * Evec.dim - - val = _np.sqrt(d / nSummands) - if "squared" in metric: - val = val ** 2 - ret += val - return ret - - elif metric == "fidelity": - # Leakage-aware metric supported for gates, using leaky_entanglement_fidelity. - for opLbl in mdl.operations: - wt = item_weights.get(opLbl, opWeight) - top = target_model.operations[opLbl].to_dense() - mop = mdl.operations[opLbl].to_dense() - ret += wt * (1.0 - _tools.leaky_entanglement_fidelity(top, mop, mxBasis, n_leak))**2 - # Leakage-aware metrics NOT available for SPAM - for preplabel, m_prep in mdl.preps.items(): - wt = item_weights.get(preplabel, spamWeight) - rhoMx1 = _tools.vec_to_stdmx(m_prep.to_dense(), mxBasis) - t_prep = target_model.preps[preplabel] - rhoMx2 = _tools.vec_to_stdmx(t_prep.to_dense(), mxBasis) - ret += wt * (1.0 - _tools.fidelity(rhoMx1, rhoMx2))**2 - for povmlabel in mdl.povms.keys(): - wt = item_weights.get(povmlabel, spamWeight) - fidelity = _tools.povm_fidelity(mdl, target_model, povmlabel) - ret += wt * (1.0 - fidelity)**2 - return ret - - elif metric == "tracedist": - # Leakage-aware metric supported, using leaky_jtracedist. - for opLbl in mdl.operations: - wt = item_weights.get(opLbl, opWeight) - top = target_model.operations[opLbl].to_dense() - mop = mdl.operations[opLbl].to_dense() - ret += wt * _tools.leaky_jtracedist(top, mop, mxBasis, n_leak) - # Leakage-aware metrics NOT available for SPAM - for preplabel, m_prep in mdl.preps.items(): - wt = item_weights.get(preplabel, spamWeight) - rhoMx1 = _tools.vec_to_stdmx(m_prep.to_dense(), mxBasis) - t_prep = target_model.preps[preplabel] - rhoMx2 = _tools.vec_to_stdmx(t_prep.to_dense(), mxBasis) - ret += wt * _tools.tracedist(rhoMx1, rhoMx2) - for povmlabel in mdl.povms.keys(): - wt = item_weights.get(povmlabel, spamWeight) - ret += wt * _tools.povm_jtracedist(mdl, target_model, povmlabel) - return ret + for lbl, Evec in mdl.effects.items(): + wt = item_weights.get(lbl, spamWeight) + evec = Evec.to_dense() + other = target_model.effects[lbl].to_dense() + delta = evec - other + delta = delta @ P + val = _np.linalg.norm(delta.flatten()) + d += wt * val**2 + nSummands += wt * Evec.dim + + val = _np.sqrt(d / nSummands) + if "squared" in metric: + val = val ** 2 + ret += val + return ret - else: - raise ValueError("Invalid metric: %s" % metric) # end _objective_fn _jacobian_fn = None From 805c9ee96ac29a6f05ed7326f6a33aed372a2f5d Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Thu, 21 Nov 2024 16:36:20 -0500 Subject: [PATCH 44/51] status for notebook 2024-11-21a.ipynb --- pygsti/algorithms/gaugeopt.py | 159 ++++++++++++++++++++++++--------- pygsti/optimize/optimize.py | 2 +- pygsti/report/reportables.py | 16 +++- pygsti/report/section/gauge.py | 2 +- pygsti/tools/optools.py | 13 +++ 5 files changed, 147 insertions(+), 45 deletions(-) diff --git a/pygsti/algorithms/gaugeopt.py b/pygsti/algorithms/gaugeopt.py index 1b9e9dda3..fc5982057 100644 --- a/pygsti/algorithms/gaugeopt.py +++ b/pygsti/algorithms/gaugeopt.py @@ -23,6 +23,39 @@ from pygsti.models.gaugegroup import TrivialGaugeGroupElement as _TrivialGaugeGroupElement +class LoggingFriendlyGaugeoptObjective: + + def __init__(self, callme, loglist): + self.callme = callme + self.loglist = loglist + self.calls_since_callabackcall = 0 + self.log = None + self.num_callbackcalls = 0 + self.first_call = None + + + def __call__(self, gaugeGroupEl, oob_check): + self.calls_since_callabackcall += 1 + val = self.callme(gaugeGroupEl, oob_check) + if self.first_call is None: + self.first_call = self.loglist[0] + return val + + def callback(self, arg): + keep = self.loglist[-1] + self.num_callbackcalls += 1 + for _ in range(self.calls_since_callabackcall): + self.loglist.pop() + self.loglist.append(keep) + self.calls_since_callabackcall = 0 + + def finalize(self): + self.loglist.insert(0, self.first_call) + temp = _np.array(self.loglist) + temp = _np.swapaxes(_np.swapaxes(temp,0,2),1,2) + self.log = temp + + def gaugeopt_to_target(model, target_model, item_weights=None, cptp_penalty_factor=0, spam_penalty_factor=0, gates_metric="frobenius", spam_metric="frobenius", @@ -167,20 +200,30 @@ def gaugeopt_to_target(model, target_model, item_weights=None, else: raise ValueError("Invalid `convert_model_to` arguments: %s" % str(args)) - objective_fn, jacobian_fn = _create_objective_fn( + objective_fn, jacobian_fn, iter_logs = _create_objective_fn( model, target_model, item_weights, cptp_penalty_factor, spam_penalty_factor, gates_metric, spam_metric, method, comm, check_jac, n_leak) + objective_fn = LoggingFriendlyGaugeoptObjective(objective_fn, iter_logs) + result = gaugeopt_custom(model, objective_fn, gauge_group, method, maxiter, maxfev, tol, oob_check_interval, return_all, jacobian_fn, comm, verbosity) + objective_fn.finalize() + _np.set_printoptions(linewidth=200, edgeitems=20) + print(objective_fn.log) + # if len(iter_logs) > 0: + # iter_logs = _np.array(iter_logs) + # _np.set_printoptions(linewidth=200, edgeitems=20) + # print(iter_logs) + #If we've gauge optimized to a target model, declare that the # resulting model is now in the same basis as the target. - if target_model is not None: - newModel = result[-1] if return_all else result - newModel.basis = target_model.basis.copy() + # if target_model is not None: + # newModel = result[-1] if return_all else result + # newModel.basis = target_model.basis.copy() return result @@ -331,7 +374,7 @@ def _call_jacobian_fn(gauge_group_el_vec): minSol = _opt.minimize(_call_objective_fn, x0, method=method, maxiter=maxiter, maxfev=maxfev, tol=tol, jac=_call_jacobian_fn, - callback=print_obj_func) + callback=objective_fn.callback) solnX = minSol.x solnF = minSol.fun @@ -355,7 +398,8 @@ def _create_objective_fn(model, target_model, item_weights=None, Creates the objective function and jacobian (if available) for gaugeopt_to_target """ - if item_weights is None: item_weights = {} + if item_weights is None: + item_weights = {} opWeight = item_weights.get('gates', 1.0) spamWeight = item_weights.get('spam', 1.0) mxBasis = model.basis @@ -377,6 +421,8 @@ def _transform_with_oob_check(mdl, gauge_group_el, oob_check): mdl.transform_inplace(gauge_group_el) return mdl + iteration_log = [] + if method == "ls": # least-squares case where objective function returns an array of # the before-they're-squared difference terms and there's an analytic jacobian @@ -599,65 +645,94 @@ def _mock_objective_fn(v): metric = spam_metric dim = int(_np.sqrt(mxBasis.dim)) - B = _tools.leading_dxd_submatrix_basis_vectors(dim - n_leak, dim, mxBasis) - P = B @ B.T.conj() - assert _np.linalg.norm(P.imag) <= 1e-12 - P = P.real + if dim % 3 == 0: + B = _tools.leading_dxd_submatrix_basis_vectors(dim - 1, dim, mxBasis) + P = B @ B.T.conj() + assert _np.linalg.norm(P.imag) <= 1e-12 + P = P.real + else: + P = _np.eye(dim) + - log = [] + elements_per_iter_log = 1 + len(model.operations) + if spamWeight > 0: + elements_per_iter_log += len(model.preps) + elements_per_iter_log += len(model.effects) + if cptp_penalty_factor > 0: + elements_per_iter_log += 1 + if spam_penalty_factor > 0: + elements_per_iter_log += 1 def _objective_fn(gauge_group_el, oob_check): mdl = _transform_with_oob_check(model, gauge_group_el, oob_check) - ret = 0 - if cptp_penalty_factor > 0: - mdl.basis = mxBasis # set basis for jamiolkowski iso - cpPenaltyVec = _cptp_penalty(mdl, cptp_penalty_factor, mdl.basis) - ret += _np.sum(cpPenaltyVec) - - if spam_penalty_factor > 0: - mdl.basis = mxBasis - spamPenaltyVec = _spam_penalty(mdl, spam_penalty_factor, mdl.basis) - ret += _np.sum(spamPenaltyVec) + iter_log = _np.zeros((elements_per_iter_log,2)) - d = 0 nSummands = 0.0 + i = 1 for opLabel, gate in mdl.operations.items(): wt = item_weights.get(opLabel, opWeight) gate_mx = gate.to_dense() other_mx = target_model.operations[opLabel].to_dense() delta = gate_mx - other_mx + iter_log[i,0] = wt * _np.linalg.norm(delta)**2 delta = delta @ P - val = _np.linalg.norm(delta.flatten()) - d += wt * val**2 nSummands += wt * (gate.dim)**2 + val = _np.linalg.norm(delta) + iter_log[i,1] = wt * val**2 + i += 1 + + if spamWeight > 0: + for lbl, rhoV in mdl.preps.items(): + wt = item_weights.get(lbl, spamWeight) + nSummands += wt * rhoV.dim + val = rhoV.frobeniusdist_squared(target_model.preps[lbl], None, None) + iter_log[i,:] = wt * val + i += 1 + + for lbl, Evec in mdl.effects.items(): + wt = item_weights.get(lbl, spamWeight) + evec = Evec.to_dense() + other = target_model.effects[lbl].to_dense() + delta = evec - other + iter_log[i,0] = wt * _np.linalg.norm(delta)**2 + delta = delta @ P + nSummands += wt * Evec.dim + val = _np.linalg.norm(delta) + iter_log[i,1] = wt * val**2 + i += 1 + + iter_log[0,:] = _np.sqrt(_np.sum(iter_log[1:,:], axis=0) / nSummands) - for lbl, rhoV in mdl.preps.items(): - wt = item_weights.get(lbl, spamWeight) - d += wt * rhoV.frobeniusdist_squared(target_model.preps[lbl], None, None) - nSummands += wt * rhoV.dim + if "squared" in metric: + iter_log[0,:] *= iter_log[0,:] - for lbl, Evec in mdl.effects.items(): - wt = item_weights.get(lbl, spamWeight) - evec = Evec.to_dense() - other = target_model.effects[lbl].to_dense() - delta = evec - other - delta = delta @ P - val = _np.linalg.norm(delta.flatten()) - d += wt * val**2 - nSummands += wt * Evec.dim + if cptp_penalty_factor > 0: + mdl.basis = mxBasis # set basis for jamiolkowski iso + cpPenaltyVec = _cptp_penalty(mdl, cptp_penalty_factor, mdl.basis) + val = _np.sum(cpPenaltyVec) + iter_log[i,:] = val + iter_log[0,:] += val + i += 1 - val = _np.sqrt(d / nSummands) - if "squared" in metric: - val = val ** 2 - ret += val + if spam_penalty_factor > 0: + mdl.basis = mxBasis + spamPenaltyVec = _spam_penalty(mdl, spam_penalty_factor, mdl.basis) + val = _np.sum(spamPenaltyVec) + iter_log[i,:] = val + iter_log[0,:] += val + i += 1 + + iteration_log.append(iter_log) + + ret = iter_log[0,0] if n_leak == 0 else iter_log[0,1] return ret # end _objective_fn _jacobian_fn = None - return _objective_fn, _jacobian_fn + return _objective_fn, _jacobian_fn, iteration_log def _cptp_penalty_size(mdl): diff --git a/pygsti/optimize/optimize.py b/pygsti/optimize/optimize.py index a6a6a982f..664d4a518 100644 --- a/pygsti/optimize/optimize.py +++ b/pygsti/optimize/optimize.py @@ -156,7 +156,7 @@ def _basin_callback(x, f, accept): #Set options for different algorithms opts = {'maxiter': maxiter, 'disp': False} if method == "BFGS": opts['gtol'] = tol # gradient norm tolerance - elif method == "L-BFGS-B": opts['gtol'] = opts['ftol'] = tol # gradient norm and fractional y-tolerance + elif method == "L-BFGS-B": opts['gtol'] = opts['ftol'] = tol/1e4 # gradient norm and fractional y-tolerance elif method == "Nelder-Mead": opts['maxfev'] = maxfev # max fn evals (note: ftol and xtol can also be set) solution = _spo.minimize(fn, x0, options=opts, method=method, tol=tol, callback=callback, jac=jac) diff --git a/pygsti/report/reportables.py b/pygsti/report/reportables.py index b4c5eeed6..6ec1cf86f 100644 --- a/pygsti/report/reportables.py +++ b/pygsti/report/reportables.py @@ -1071,6 +1071,15 @@ def leaky_entanglement_fidelity(a, b, mx_basis): Leaky_entanglement_fidelity = _modf.opsfn_factory(leaky_entanglement_fidelity) +def leaky_gate_frob_dist(a, b, mx_basis): + n_leak = 1 + norm_a = 1.0 # _tools.subspace_restricted_fro_dist(a, _np.zeros_like(b), n_leak) + norm_b = 1.0 # _tools.subspace_restricted_fro_dist(_np.zeros_like(a), b, n_leak) + return (_tools.subspace_restricted_fro_dist(a, b, mx_basis, n_leak) / (0.5*(norm_a + norm_b)))**2 + +Leaky_gate_frob_dist = _modf.opsfn_factory(leaky_gate_frob_dist) + + def entanglement_infidelity(a, b, mx_basis): """ Entanglement infidelity between a and b @@ -1160,7 +1169,7 @@ def frobenius_diff(a, b, mx_basis): # assume vary model1, model2 fixed ------- float """ - return _tools.frobeniusdist(a, b) + return _tools.frobeniusdist(a, b)**2 Fro_diff = _modf.opsfn_factory(frobenius_diff) @@ -2455,6 +2464,8 @@ def info_of_opfn_by_name(name): "where (a_i,b_i) are corresponding eigenvalues of A and B."), "frob": ("Frobenius|Distance", "sqrt( sum( (A_ij - B_ij)^2 ) )"), + "la-frob": ("Frobenius|Distance (subspace)", + "TO WRITE"), "unmodeled": ("Un-modeled|Error", "The per-operation budget used to account for un-modeled errors (model violation)") } @@ -2539,6 +2550,9 @@ def evaluate_opfn_by_name(name, model, target_model, op_label_or_string, elif name == "evnudiamond": fn = Eigenvalue_nonunitary_diamondnorm if b else \ Circuit_eigenvalue_nonunitary_diamondnorm + elif name == "la-frob": + assert b + fn = Leaky_gate_frob_dist elif name == "frob": fn = Fro_diff if b else \ Circuit_fro_diff diff --git a/pygsti/report/section/gauge.py b/pygsti/report/section/gauge.py index c9c642e81..5927f3d2f 100644 --- a/pygsti/report/section/gauge.py +++ b/pygsti/report/section/gauge.py @@ -171,7 +171,7 @@ def final_gates_vs_target_table_gauge_var(workspace, switchboard=None, confidenc if kwargs.get('n_leak', 0) == 0: display = ('inf', 'agi', 'trace', 'diamond', 'nuinf', 'nuagi') else: - display = ('inf', 'la-inf', 'agi', 'trace', 'la-trace', 'diamond', 'nuinf', 'nuagi') + display = ('la-frob', 'la-inf', 'la-trace', 'frob', 'inf', 'trace', 'agi', 'diamond', 'nuinf', 'nuagi') return workspace.GatesVsTargetTable( switchboard.mdl_final, switchboard.mdl_target, _cri(1, switchboard, confidence_level, ci_brevity), display=display diff --git a/pygsti/tools/optools.py b/pygsti/tools/optools.py index 6a151848c..c7b11b6ea 100644 --- a/pygsti/tools/optools.py +++ b/pygsti/tools/optools.py @@ -618,6 +618,19 @@ def leading_dxd_submatrix_basis_vectors(d: int, n: int, current_basis, return_la return submatrix_basis_vectors +def subspace_restricted_fro_dist(a, b, mx_basis, n_leak=0): + diff = a - b + if n_leak == 0: + return _np.linalg.norm(diff, 'fro') + if n_leak == 1: + d = int(_np.sqrt(a.shape[0])) + assert a.shape == b.shape == (d**2, d**2) + B = leading_dxd_submatrix_basis_vectors(d-n_leak, d, mx_basis) + P = B @ B.T.conj() + return _np.linalg.norm(diff @ P) + raise ValueError() + + def average_gate_fidelity(a, b, mx_basis='pp', is_tp=None, is_unitary=None): """ Computes the average gate fidelity (AGF) between two gates. From 79d22817b7c350f93eb070a54155a01dfeaa994e Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Fri, 22 Nov 2024 08:30:45 -0500 Subject: [PATCH 45/51] bugfix in leading_dxd_submatrix_basis_vectors helper function --- pygsti/tools/optools.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/pygsti/tools/optools.py b/pygsti/tools/optools.py index c7b11b6ea..00820ec90 100644 --- a/pygsti/tools/optools.py +++ b/pygsti/tools/optools.py @@ -597,9 +597,10 @@ def leading_dxd_submatrix_basis_vectors(d: int, n: int, current_basis, return_la """ assert d <= n current_basis = _pgb.Basis.cast(current_basis, dim=n**2) - std_to_current = current_basis.create_transform_matrix('std') + X = current_basis.create_transform_matrix('std') + X = X.T.conj() if d == n: - return std_to_current + return X # we have to select a proper subset of columns in current_basis std_basis = _pgb.BuiltinBasis(name='std', dim_or_statespace=n**2) label2ind = {ell: idx for idx,ell in enumerate(std_basis.labels)} @@ -611,7 +612,7 @@ def leading_dxd_submatrix_basis_vectors(d: int, n: int, current_basis, return_la basis_ind.append(label2ind[ell]) basis_labels.append(ell) basis_ind = _np.array(basis_ind) - submatrix_basis_vectors = std_to_current[:, basis_ind] + submatrix_basis_vectors = X[:, basis_ind] if return_labels: return submatrix_basis_vectors, basis_labels if not return_labels: From 2561b7c588ac64b15b6c4ce750ef9c6d45153bf5 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Fri, 22 Nov 2024 08:56:49 -0500 Subject: [PATCH 46/51] Revert "status for notebook 2024-11-21a.ipynb" This reverts commit dfb1118230565837fcafb252035a5f9c4a7af9aa. --- pygsti/algorithms/gaugeopt.py | 159 +++++++++------------------------ pygsti/optimize/optimize.py | 2 +- pygsti/report/reportables.py | 16 +--- pygsti/report/section/gauge.py | 2 +- pygsti/tools/optools.py | 13 --- 5 files changed, 45 insertions(+), 147 deletions(-) diff --git a/pygsti/algorithms/gaugeopt.py b/pygsti/algorithms/gaugeopt.py index fc5982057..1b9e9dda3 100644 --- a/pygsti/algorithms/gaugeopt.py +++ b/pygsti/algorithms/gaugeopt.py @@ -23,39 +23,6 @@ from pygsti.models.gaugegroup import TrivialGaugeGroupElement as _TrivialGaugeGroupElement -class LoggingFriendlyGaugeoptObjective: - - def __init__(self, callme, loglist): - self.callme = callme - self.loglist = loglist - self.calls_since_callabackcall = 0 - self.log = None - self.num_callbackcalls = 0 - self.first_call = None - - - def __call__(self, gaugeGroupEl, oob_check): - self.calls_since_callabackcall += 1 - val = self.callme(gaugeGroupEl, oob_check) - if self.first_call is None: - self.first_call = self.loglist[0] - return val - - def callback(self, arg): - keep = self.loglist[-1] - self.num_callbackcalls += 1 - for _ in range(self.calls_since_callabackcall): - self.loglist.pop() - self.loglist.append(keep) - self.calls_since_callabackcall = 0 - - def finalize(self): - self.loglist.insert(0, self.first_call) - temp = _np.array(self.loglist) - temp = _np.swapaxes(_np.swapaxes(temp,0,2),1,2) - self.log = temp - - def gaugeopt_to_target(model, target_model, item_weights=None, cptp_penalty_factor=0, spam_penalty_factor=0, gates_metric="frobenius", spam_metric="frobenius", @@ -200,30 +167,20 @@ def gaugeopt_to_target(model, target_model, item_weights=None, else: raise ValueError("Invalid `convert_model_to` arguments: %s" % str(args)) - objective_fn, jacobian_fn, iter_logs = _create_objective_fn( + objective_fn, jacobian_fn = _create_objective_fn( model, target_model, item_weights, cptp_penalty_factor, spam_penalty_factor, gates_metric, spam_metric, method, comm, check_jac, n_leak) - objective_fn = LoggingFriendlyGaugeoptObjective(objective_fn, iter_logs) - result = gaugeopt_custom(model, objective_fn, gauge_group, method, maxiter, maxfev, tol, oob_check_interval, return_all, jacobian_fn, comm, verbosity) - objective_fn.finalize() - _np.set_printoptions(linewidth=200, edgeitems=20) - print(objective_fn.log) - # if len(iter_logs) > 0: - # iter_logs = _np.array(iter_logs) - # _np.set_printoptions(linewidth=200, edgeitems=20) - # print(iter_logs) - #If we've gauge optimized to a target model, declare that the # resulting model is now in the same basis as the target. - # if target_model is not None: - # newModel = result[-1] if return_all else result - # newModel.basis = target_model.basis.copy() + if target_model is not None: + newModel = result[-1] if return_all else result + newModel.basis = target_model.basis.copy() return result @@ -374,7 +331,7 @@ def _call_jacobian_fn(gauge_group_el_vec): minSol = _opt.minimize(_call_objective_fn, x0, method=method, maxiter=maxiter, maxfev=maxfev, tol=tol, jac=_call_jacobian_fn, - callback=objective_fn.callback) + callback=print_obj_func) solnX = minSol.x solnF = minSol.fun @@ -398,8 +355,7 @@ def _create_objective_fn(model, target_model, item_weights=None, Creates the objective function and jacobian (if available) for gaugeopt_to_target """ - if item_weights is None: - item_weights = {} + if item_weights is None: item_weights = {} opWeight = item_weights.get('gates', 1.0) spamWeight = item_weights.get('spam', 1.0) mxBasis = model.basis @@ -421,8 +377,6 @@ def _transform_with_oob_check(mdl, gauge_group_el, oob_check): mdl.transform_inplace(gauge_group_el) return mdl - iteration_log = [] - if method == "ls": # least-squares case where objective function returns an array of # the before-they're-squared difference terms and there's an analytic jacobian @@ -645,94 +599,65 @@ def _mock_objective_fn(v): metric = spam_metric dim = int(_np.sqrt(mxBasis.dim)) - if dim % 3 == 0: - B = _tools.leading_dxd_submatrix_basis_vectors(dim - 1, dim, mxBasis) - P = B @ B.T.conj() - assert _np.linalg.norm(P.imag) <= 1e-12 - P = P.real - else: - P = _np.eye(dim) - + B = _tools.leading_dxd_submatrix_basis_vectors(dim - n_leak, dim, mxBasis) + P = B @ B.T.conj() + assert _np.linalg.norm(P.imag) <= 1e-12 + P = P.real - elements_per_iter_log = 1 + len(model.operations) - if spamWeight > 0: - elements_per_iter_log += len(model.preps) - elements_per_iter_log += len(model.effects) - if cptp_penalty_factor > 0: - elements_per_iter_log += 1 - if spam_penalty_factor > 0: - elements_per_iter_log += 1 + log = [] def _objective_fn(gauge_group_el, oob_check): mdl = _transform_with_oob_check(model, gauge_group_el, oob_check) + ret = 0 - iter_log = _np.zeros((elements_per_iter_log,2)) + if cptp_penalty_factor > 0: + mdl.basis = mxBasis # set basis for jamiolkowski iso + cpPenaltyVec = _cptp_penalty(mdl, cptp_penalty_factor, mdl.basis) + ret += _np.sum(cpPenaltyVec) + + if spam_penalty_factor > 0: + mdl.basis = mxBasis + spamPenaltyVec = _spam_penalty(mdl, spam_penalty_factor, mdl.basis) + ret += _np.sum(spamPenaltyVec) + d = 0 nSummands = 0.0 - i = 1 for opLabel, gate in mdl.operations.items(): wt = item_weights.get(opLabel, opWeight) gate_mx = gate.to_dense() other_mx = target_model.operations[opLabel].to_dense() delta = gate_mx - other_mx - iter_log[i,0] = wt * _np.linalg.norm(delta)**2 delta = delta @ P + val = _np.linalg.norm(delta.flatten()) + d += wt * val**2 nSummands += wt * (gate.dim)**2 - val = _np.linalg.norm(delta) - iter_log[i,1] = wt * val**2 - i += 1 - - if spamWeight > 0: - for lbl, rhoV in mdl.preps.items(): - wt = item_weights.get(lbl, spamWeight) - nSummands += wt * rhoV.dim - val = rhoV.frobeniusdist_squared(target_model.preps[lbl], None, None) - iter_log[i,:] = wt * val - i += 1 - - for lbl, Evec in mdl.effects.items(): - wt = item_weights.get(lbl, spamWeight) - evec = Evec.to_dense() - other = target_model.effects[lbl].to_dense() - delta = evec - other - iter_log[i,0] = wt * _np.linalg.norm(delta)**2 - delta = delta @ P - nSummands += wt * Evec.dim - val = _np.linalg.norm(delta) - iter_log[i,1] = wt * val**2 - i += 1 - - iter_log[0,:] = _np.sqrt(_np.sum(iter_log[1:,:], axis=0) / nSummands) - - if "squared" in metric: - iter_log[0,:] *= iter_log[0,:] - - if cptp_penalty_factor > 0: - mdl.basis = mxBasis # set basis for jamiolkowski iso - cpPenaltyVec = _cptp_penalty(mdl, cptp_penalty_factor, mdl.basis) - val = _np.sum(cpPenaltyVec) - iter_log[i,:] = val - iter_log[0,:] += val - i += 1 - if spam_penalty_factor > 0: - mdl.basis = mxBasis - spamPenaltyVec = _spam_penalty(mdl, spam_penalty_factor, mdl.basis) - val = _np.sum(spamPenaltyVec) - iter_log[i,:] = val - iter_log[0,:] += val - i += 1 + for lbl, rhoV in mdl.preps.items(): + wt = item_weights.get(lbl, spamWeight) + d += wt * rhoV.frobeniusdist_squared(target_model.preps[lbl], None, None) + nSummands += wt * rhoV.dim - iteration_log.append(iter_log) + for lbl, Evec in mdl.effects.items(): + wt = item_weights.get(lbl, spamWeight) + evec = Evec.to_dense() + other = target_model.effects[lbl].to_dense() + delta = evec - other + delta = delta @ P + val = _np.linalg.norm(delta.flatten()) + d += wt * val**2 + nSummands += wt * Evec.dim - ret = iter_log[0,0] if n_leak == 0 else iter_log[0,1] + val = _np.sqrt(d / nSummands) + if "squared" in metric: + val = val ** 2 + ret += val return ret # end _objective_fn _jacobian_fn = None - return _objective_fn, _jacobian_fn, iteration_log + return _objective_fn, _jacobian_fn def _cptp_penalty_size(mdl): diff --git a/pygsti/optimize/optimize.py b/pygsti/optimize/optimize.py index 664d4a518..a6a6a982f 100644 --- a/pygsti/optimize/optimize.py +++ b/pygsti/optimize/optimize.py @@ -156,7 +156,7 @@ def _basin_callback(x, f, accept): #Set options for different algorithms opts = {'maxiter': maxiter, 'disp': False} if method == "BFGS": opts['gtol'] = tol # gradient norm tolerance - elif method == "L-BFGS-B": opts['gtol'] = opts['ftol'] = tol/1e4 # gradient norm and fractional y-tolerance + elif method == "L-BFGS-B": opts['gtol'] = opts['ftol'] = tol # gradient norm and fractional y-tolerance elif method == "Nelder-Mead": opts['maxfev'] = maxfev # max fn evals (note: ftol and xtol can also be set) solution = _spo.minimize(fn, x0, options=opts, method=method, tol=tol, callback=callback, jac=jac) diff --git a/pygsti/report/reportables.py b/pygsti/report/reportables.py index 6ec1cf86f..b4c5eeed6 100644 --- a/pygsti/report/reportables.py +++ b/pygsti/report/reportables.py @@ -1071,15 +1071,6 @@ def leaky_entanglement_fidelity(a, b, mx_basis): Leaky_entanglement_fidelity = _modf.opsfn_factory(leaky_entanglement_fidelity) -def leaky_gate_frob_dist(a, b, mx_basis): - n_leak = 1 - norm_a = 1.0 # _tools.subspace_restricted_fro_dist(a, _np.zeros_like(b), n_leak) - norm_b = 1.0 # _tools.subspace_restricted_fro_dist(_np.zeros_like(a), b, n_leak) - return (_tools.subspace_restricted_fro_dist(a, b, mx_basis, n_leak) / (0.5*(norm_a + norm_b)))**2 - -Leaky_gate_frob_dist = _modf.opsfn_factory(leaky_gate_frob_dist) - - def entanglement_infidelity(a, b, mx_basis): """ Entanglement infidelity between a and b @@ -1169,7 +1160,7 @@ def frobenius_diff(a, b, mx_basis): # assume vary model1, model2 fixed ------- float """ - return _tools.frobeniusdist(a, b)**2 + return _tools.frobeniusdist(a, b) Fro_diff = _modf.opsfn_factory(frobenius_diff) @@ -2464,8 +2455,6 @@ def info_of_opfn_by_name(name): "where (a_i,b_i) are corresponding eigenvalues of A and B."), "frob": ("Frobenius|Distance", "sqrt( sum( (A_ij - B_ij)^2 ) )"), - "la-frob": ("Frobenius|Distance (subspace)", - "TO WRITE"), "unmodeled": ("Un-modeled|Error", "The per-operation budget used to account for un-modeled errors (model violation)") } @@ -2550,9 +2539,6 @@ def evaluate_opfn_by_name(name, model, target_model, op_label_or_string, elif name == "evnudiamond": fn = Eigenvalue_nonunitary_diamondnorm if b else \ Circuit_eigenvalue_nonunitary_diamondnorm - elif name == "la-frob": - assert b - fn = Leaky_gate_frob_dist elif name == "frob": fn = Fro_diff if b else \ Circuit_fro_diff diff --git a/pygsti/report/section/gauge.py b/pygsti/report/section/gauge.py index 5927f3d2f..c9c642e81 100644 --- a/pygsti/report/section/gauge.py +++ b/pygsti/report/section/gauge.py @@ -171,7 +171,7 @@ def final_gates_vs_target_table_gauge_var(workspace, switchboard=None, confidenc if kwargs.get('n_leak', 0) == 0: display = ('inf', 'agi', 'trace', 'diamond', 'nuinf', 'nuagi') else: - display = ('la-frob', 'la-inf', 'la-trace', 'frob', 'inf', 'trace', 'agi', 'diamond', 'nuinf', 'nuagi') + display = ('inf', 'la-inf', 'agi', 'trace', 'la-trace', 'diamond', 'nuinf', 'nuagi') return workspace.GatesVsTargetTable( switchboard.mdl_final, switchboard.mdl_target, _cri(1, switchboard, confidence_level, ci_brevity), display=display diff --git a/pygsti/tools/optools.py b/pygsti/tools/optools.py index 00820ec90..9d956f18f 100644 --- a/pygsti/tools/optools.py +++ b/pygsti/tools/optools.py @@ -619,19 +619,6 @@ def leading_dxd_submatrix_basis_vectors(d: int, n: int, current_basis, return_la return submatrix_basis_vectors -def subspace_restricted_fro_dist(a, b, mx_basis, n_leak=0): - diff = a - b - if n_leak == 0: - return _np.linalg.norm(diff, 'fro') - if n_leak == 1: - d = int(_np.sqrt(a.shape[0])) - assert a.shape == b.shape == (d**2, d**2) - B = leading_dxd_submatrix_basis_vectors(d-n_leak, d, mx_basis) - P = B @ B.T.conj() - return _np.linalg.norm(diff @ P) - raise ValueError() - - def average_gate_fidelity(a, b, mx_basis='pp', is_tp=None, is_unitary=None): """ Computes the average gate fidelity (AGF) between two gates. From 35bc772a7e06dad33205ba1c8d0d99954fd6895d Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Fri, 22 Nov 2024 08:56:57 -0500 Subject: [PATCH 47/51] Revert "remove non-Frobenius objectives from gauge optimization" This reverts commit 578326394f226a87928823009c42a8cf5f82a360. --- pygsti/algorithms/gaugeopt.py | 106 ++++++++++++++++++++++++---------- 1 file changed, 74 insertions(+), 32 deletions(-) diff --git a/pygsti/algorithms/gaugeopt.py b/pygsti/algorithms/gaugeopt.py index 1b9e9dda3..51225a83a 100644 --- a/pygsti/algorithms/gaugeopt.py +++ b/pygsti/algorithms/gaugeopt.py @@ -619,40 +619,82 @@ def _objective_fn(gauge_group_el, oob_check): mdl.basis = mxBasis spamPenaltyVec = _spam_penalty(mdl, spam_penalty_factor, mdl.basis) ret += _np.sum(spamPenaltyVec) - - d = 0 - nSummands = 0.0 - for opLabel, gate in mdl.operations.items(): - wt = item_weights.get(opLabel, opWeight) - gate_mx = gate.to_dense() - other_mx = target_model.operations[opLabel].to_dense() - delta = gate_mx - other_mx - delta = delta @ P - val = _np.linalg.norm(delta.flatten()) - d += wt * val**2 - nSummands += wt * (gate.dim)**2 - - for lbl, rhoV in mdl.preps.items(): - wt = item_weights.get(lbl, spamWeight) - d += wt * rhoV.frobeniusdist_squared(target_model.preps[lbl], None, None) - nSummands += wt * rhoV.dim - for lbl, Evec in mdl.effects.items(): - wt = item_weights.get(lbl, spamWeight) - evec = Evec.to_dense() - other = target_model.effects[lbl].to_dense() - delta = evec - other - delta = delta @ P - val = _np.linalg.norm(delta.flatten()) - d += wt * val**2 - nSummands += wt * Evec.dim - - val = _np.sqrt(d / nSummands) - if "squared" in metric: - val = val ** 2 - ret += val - return ret + if "frobenius" in metric: + d = 0 + nSummands = 0.0 + for opLabel, gate in mdl.operations.items(): + wt = item_weights.get(opLabel, opWeight) + gate_mx = gate.to_dense() + other_mx = target_model.operations[opLabel].to_dense() + delta = gate_mx - other_mx + delta = delta @ P + val = _np.linalg.norm(delta.flatten()) + d += wt * val**2 + nSummands += wt * (gate.dim)**2 + + for lbl, rhoV in mdl.preps.items(): + wt = item_weights.get(lbl, spamWeight) + d += wt * rhoV.frobeniusdist_squared(target_model.preps[lbl], None, None) + nSummands += wt * rhoV.dim + + for lbl, Evec in mdl.effects.items(): + wt = item_weights.get(lbl, spamWeight) + evec = Evec.to_dense() + other = target_model.effects[lbl].to_dense() + delta = evec - other + delta = delta @ P + val = _np.linalg.norm(delta.flatten()) + d += wt * val**2 + nSummands += wt * Evec.dim + + val = _np.sqrt(d / nSummands) + if "squared" in metric: + val = val ** 2 + ret += val + return ret + + elif metric == "fidelity": + # Leakage-aware metric supported for gates, using leaky_entanglement_fidelity. + for opLbl in mdl.operations: + wt = item_weights.get(opLbl, opWeight) + top = target_model.operations[opLbl].to_dense() + mop = mdl.operations[opLbl].to_dense() + ret += wt * (1.0 - _tools.leaky_entanglement_fidelity(top, mop, mxBasis, n_leak))**2 + # Leakage-aware metrics NOT available for SPAM + for preplabel, m_prep in mdl.preps.items(): + wt = item_weights.get(preplabel, spamWeight) + rhoMx1 = _tools.vec_to_stdmx(m_prep.to_dense(), mxBasis) + t_prep = target_model.preps[preplabel] + rhoMx2 = _tools.vec_to_stdmx(t_prep.to_dense(), mxBasis) + ret += wt * (1.0 - _tools.fidelity(rhoMx1, rhoMx2))**2 + for povmlabel in mdl.povms.keys(): + wt = item_weights.get(povmlabel, spamWeight) + fidelity = _tools.povm_fidelity(mdl, target_model, povmlabel) + ret += wt * (1.0 - fidelity)**2 + return ret + + elif metric == "tracedist": + # Leakage-aware metric supported, using leaky_jtracedist. + for opLbl in mdl.operations: + wt = item_weights.get(opLbl, opWeight) + top = target_model.operations[opLbl].to_dense() + mop = mdl.operations[opLbl].to_dense() + ret += wt * _tools.leaky_jtracedist(top, mop, mxBasis, n_leak) + # Leakage-aware metrics NOT available for SPAM + for preplabel, m_prep in mdl.preps.items(): + wt = item_weights.get(preplabel, spamWeight) + rhoMx1 = _tools.vec_to_stdmx(m_prep.to_dense(), mxBasis) + t_prep = target_model.preps[preplabel] + rhoMx2 = _tools.vec_to_stdmx(t_prep.to_dense(), mxBasis) + ret += wt * _tools.tracedist(rhoMx1, rhoMx2) + for povmlabel in mdl.povms.keys(): + wt = item_weights.get(povmlabel, spamWeight) + ret += wt * _tools.povm_jtracedist(mdl, target_model, povmlabel) + return ret + else: + raise ValueError("Invalid metric: %s" % metric) # end _objective_fn _jacobian_fn = None From b598e73484a35f2b833fea2243d83a357144afbe Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Fri, 22 Nov 2024 08:57:01 -0500 Subject: [PATCH 48/51] Revert "simplify gaugeopt in the non-LS case" This reverts commit 6e9ccb10096e36c6a8a846b1d87952337363f030. --- pygsti/algorithms/gaugeopt.py | 91 ++++++++++++++++++++--------------- 1 file changed, 51 insertions(+), 40 deletions(-) diff --git a/pygsti/algorithms/gaugeopt.py b/pygsti/algorithms/gaugeopt.py index 51225a83a..823d98c99 100644 --- a/pygsti/algorithms/gaugeopt.py +++ b/pygsti/algorithms/gaugeopt.py @@ -587,24 +587,45 @@ def _mock_objective_fn(v): else: # non-least-squares case where objective function returns a single float # and (currently) there's no analytic jacobian + dim = int(_np.sqrt(mxBasis.dim)) + if n_leak > 0: + B = _tools.leading_dxd_submatrix_basis_vectors(dim - n_leak, dim, mxBasis) + """ + ^ Need to do something else. + + In a leakage-friendly basis the operation matrices can be partitioned as + [comp , semileak1] + [semileak2, fulleak ]. + Instead of projecting only on to comp, we want to keep everything + except fullleak. ... But I don't think we can implement that just by + pre-multiplying by a projector. + + I think this distinction might not be significant under certain idealized + assumptions, but small deviations from those conditions (which are certain + to happen in practice) might make it matter. + """ + P = B @ B.T.conj() + if _np.linalg.norm(P.imag) > 1e-12: + raise ValueError() + else: + P = P.real + transform_mx_arg = (P, None) + # ^ The semantics of this tuple are defined by the frobeniusdist function + # in the ExplicitOpModelCalc class. There are intended semantics for + # the second element of the tuple, but those aren't implemented yet so + # for now I'm setting the second entry to None. -- Riley + else: + transform_mx_arg = (_np.eye(mxBasis.dim), None) - assert target_model is not None assert gates_metric != "frobeniustt" assert spam_metric != "frobeniustt" + assert spam_metric == gates_metric + metric = spam_metric + # assert spam_metric == gates_metric # ^ Erik and Corey said these are rarely used. I've removed support for # them in this codepath (non-LS optimizer) in order to make it easier to # read my updated code for leakage-aware metrics. It wouldn't be hard to # add support back, but I just want to keep things simple. -- Riley - assert spam_metric == gates_metric - metric = spam_metric - - dim = int(_np.sqrt(mxBasis.dim)) - B = _tools.leading_dxd_submatrix_basis_vectors(dim - n_leak, dim, mxBasis) - P = B @ B.T.conj() - assert _np.linalg.norm(P.imag) <= 1e-12 - P = P.real - - log = [] def _objective_fn(gauge_group_el, oob_check): mdl = _transform_with_oob_check(model, gauge_group_el, oob_check) @@ -620,35 +641,25 @@ def _objective_fn(gauge_group_el, oob_check): spamPenaltyVec = _spam_penalty(mdl, spam_penalty_factor, mdl.basis) ret += _np.sum(spamPenaltyVec) + assert target_model is not None + """ + Leakage-aware metric supported, per implementation in mdl.frobeniusdist. + Refer to how mdl.frobeniusdist handles the case when transform_mx_arg + is a tuple in order to understand how the leakage-aware metric is defined. + + Idea: raise an error if mxBasis isn't leakage-friendly. + PROBLEM: if we're deep in the code then we end up needing to be + working in a basis that has the identity matrix as its + first element. The leakage-friendly basis doesn't even + have the identity matrix as an element, let alone the first element + + TODO: dig into function calls below and see where we can + access the mxBasis object. (I'm sure we can.) + Looks like it isn't accessible within ExplicitOpModelCalc objects. + It's most definitely available in ExplicitOpModel. + """ if "frobenius" in metric: - d = 0 - nSummands = 0.0 - for opLabel, gate in mdl.operations.items(): - wt = item_weights.get(opLabel, opWeight) - gate_mx = gate.to_dense() - other_mx = target_model.operations[opLabel].to_dense() - delta = gate_mx - other_mx - delta = delta @ P - val = _np.linalg.norm(delta.flatten()) - d += wt * val**2 - nSummands += wt * (gate.dim)**2 - - for lbl, rhoV in mdl.preps.items(): - wt = item_weights.get(lbl, spamWeight) - d += wt * rhoV.frobeniusdist_squared(target_model.preps[lbl], None, None) - nSummands += wt * rhoV.dim - - for lbl, Evec in mdl.effects.items(): - wt = item_weights.get(lbl, spamWeight) - evec = Evec.to_dense() - other = target_model.effects[lbl].to_dense() - delta = evec - other - delta = delta @ P - val = _np.linalg.norm(delta.flatten()) - d += wt * val**2 - nSummands += wt * Evec.dim - - val = _np.sqrt(d / nSummands) + val = mdl.frobeniusdist(target_model, transform_mx_arg, item_weights) if "squared" in metric: val = val ** 2 ret += val @@ -696,7 +707,7 @@ def _objective_fn(gauge_group_el, oob_check): else: raise ValueError("Invalid metric: %s" % metric) # end _objective_fn - + _jacobian_fn = None return _objective_fn, _jacobian_fn From 1b85e82acfd2516232813a15798e3c3fdaada86e Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Fri, 22 Nov 2024 08:57:04 -0500 Subject: [PATCH 49/51] Revert "last commit actually just included comments. Material changes are in this commit" This reverts commit 1f49e43afe0d727ba66ae5b9003a3bd48ff3d7a6. --- pygsti/algorithms/gaugeopt.py | 164 ++++++++++++++++++++-------------- 1 file changed, 96 insertions(+), 68 deletions(-) diff --git a/pygsti/algorithms/gaugeopt.py b/pygsti/algorithms/gaugeopt.py index 823d98c99..5993c152a 100644 --- a/pygsti/algorithms/gaugeopt.py +++ b/pygsti/algorithms/gaugeopt.py @@ -619,8 +619,6 @@ def _mock_objective_fn(v): assert gates_metric != "frobeniustt" assert spam_metric != "frobeniustt" - assert spam_metric == gates_metric - metric = spam_metric # assert spam_metric == gates_metric # ^ Erik and Corey said these are rarely used. I've removed support for # them in this codepath (non-LS optimizer) in order to make it easier to @@ -641,73 +639,103 @@ def _objective_fn(gauge_group_el, oob_check): spamPenaltyVec = _spam_penalty(mdl, spam_penalty_factor, mdl.basis) ret += _np.sum(spamPenaltyVec) - assert target_model is not None - """ - Leakage-aware metric supported, per implementation in mdl.frobeniusdist. - Refer to how mdl.frobeniusdist handles the case when transform_mx_arg - is a tuple in order to understand how the leakage-aware metric is defined. - - Idea: raise an error if mxBasis isn't leakage-friendly. - PROBLEM: if we're deep in the code then we end up needing to be - working in a basis that has the identity matrix as its - first element. The leakage-friendly basis doesn't even - have the identity matrix as an element, let alone the first element - - TODO: dig into function calls below and see where we can - access the mxBasis object. (I'm sure we can.) - Looks like it isn't accessible within ExplicitOpModelCalc objects. - It's most definitely available in ExplicitOpModel. - """ - if "frobenius" in metric: - val = mdl.frobeniusdist(target_model, transform_mx_arg, item_weights) - if "squared" in metric: - val = val ** 2 - ret += val - return ret - - elif metric == "fidelity": - # Leakage-aware metric supported for gates, using leaky_entanglement_fidelity. - for opLbl in mdl.operations: - wt = item_weights.get(opLbl, opWeight) - top = target_model.operations[opLbl].to_dense() - mop = mdl.operations[opLbl].to_dense() - ret += wt * (1.0 - _tools.leaky_entanglement_fidelity(top, mop, mxBasis, n_leak))**2 - # Leakage-aware metrics NOT available for SPAM - for preplabel, m_prep in mdl.preps.items(): - wt = item_weights.get(preplabel, spamWeight) - rhoMx1 = _tools.vec_to_stdmx(m_prep.to_dense(), mxBasis) - t_prep = target_model.preps[preplabel] - rhoMx2 = _tools.vec_to_stdmx(t_prep.to_dense(), mxBasis) - ret += wt * (1.0 - _tools.fidelity(rhoMx1, rhoMx2))**2 - for povmlabel in mdl.povms.keys(): - wt = item_weights.get(povmlabel, spamWeight) - fidelity = _tools.povm_fidelity(mdl, target_model, povmlabel) - ret += wt * (1.0 - fidelity)**2 - return ret - - elif metric == "tracedist": - # Leakage-aware metric supported, using leaky_jtracedist. - for opLbl in mdl.operations: - wt = item_weights.get(opLbl, opWeight) - top = target_model.operations[opLbl].to_dense() - mop = mdl.operations[opLbl].to_dense() - ret += wt * _tools.leaky_jtracedist(top, mop, mxBasis, n_leak) - # Leakage-aware metrics NOT available for SPAM - for preplabel, m_prep in mdl.preps.items(): - wt = item_weights.get(preplabel, spamWeight) - rhoMx1 = _tools.vec_to_stdmx(m_prep.to_dense(), mxBasis) - t_prep = target_model.preps[preplabel] - rhoMx2 = _tools.vec_to_stdmx(t_prep.to_dense(), mxBasis) - ret += wt * _tools.tracedist(rhoMx1, rhoMx2) - for povmlabel in mdl.povms.keys(): - wt = item_weights.get(povmlabel, spamWeight) - ret += wt * _tools.povm_jtracedist(mdl, target_model, povmlabel) - return ret + if target_model is not None: + """ + Leakage-aware metric supported, per implementation in mdl.frobeniusdist. + Refer to how mdl.frobeniusdist handles the case when transform_mx_arg + is a tuple in order to understand how the leakage-aware metric is defined. + + Idea: raise an error if mxBasis isn't leakage-friendly. + PROBLEM: if we're deep in the code then we end up needing to be + working in a basis that has the identity matrix as its + first element. The leakage-friendly basis doesn't even + have the identity matrix as an element, let alone the first element + + TODO: dig into function calls below and see where we can + access the mxBasis object. (I'm sure we can.) + Looks like it isn't accessible within ExplicitOpModelCalc objects. + It's most definitely available in ExplicitOpModel. + """ + if "frobenius" in gates_metric: + if spam_metric == gates_metric: + val = mdl.frobeniusdist(target_model, transform_mx_arg, item_weights) + else: + wts = item_weights.copy() + wts['spam'] = 0.0 + for k in wts: + if k in mdl.preps or k in mdl.povms: + wts[k] = 0.0 + val = mdl.frobeniusdist(target_model, transform_mx_arg, wts, n_leak) + if "squared" in gates_metric: + val = val ** 2 + ret += val + + elif gates_metric == "fidelity": + # Leakage-aware metric supported, using leaky_entanglement_fidelity. + for opLbl in mdl.operations: + wt = item_weights.get(opLbl, opWeight) + top = target_model.operations[opLbl].to_dense() + mop = mdl.operations[opLbl].to_dense() + ret += wt * (1.0 - _tools.leaky_entanglement_fidelity(top, mop, mxBasis, n_leak))**2 + + elif gates_metric == "tracedist": + # Leakage-aware metric supported, using leaky_jtracedist. + for opLbl in mdl.operations: + wt = item_weights.get(opLbl, opWeight) + top = target_model.operations[opLbl].to_dense() + mop = mdl.operations[opLbl].to_dense() + ret += wt * _tools.leaky_jtracedist(top, mop, mxBasis, n_leak) + + else: raise ValueError("Invalid gates_metric: %s" % gates_metric) + + if "frobenius" in spam_metric and gates_metric == spam_metric: + # We already handled SPAM error in this case. Just return. + return ret + + elif "frobenius" in spam_metric: + # Leakage-aware metric supported in principle via implementation in + # mdl.frobeniusdist (check implementation to see how it handles the + # case when transform_mx_arg is a tuple). + wts = item_weights.copy(); wts['gates'] = 0.0 + for k in wts: + if k in mdl.operations or \ + k in mdl.instruments: wts[k] = 0.0 + val = mdl.frobeniusdist(target_model, transform_mx_arg, wts) + if "squared" in spam_metric: + val = val ** 2 + ret += val + + elif spam_metric == "fidelity": + # Leakage-aware metrics NOT available + for preplabel, m_prep in mdl.preps.items(): + wt = item_weights.get(preplabel, spamWeight) + rhoMx1 = _tools.vec_to_stdmx(m_prep.to_dense(), mxBasis) + t_prep = target_model.preps[preplabel] + rhoMx2 = _tools.vec_to_stdmx(t_prep.to_dense(), mxBasis) + ret += wt * (1.0 - _tools.fidelity(rhoMx1, rhoMx2))**2 + + for povmlabel in mdl.povms.keys(): + wt = item_weights.get(povmlabel, spamWeight) + fidelity = _tools.povm_fidelity(mdl, target_model, povmlabel) + ret += wt * (1.0 - fidelity)**2 + + elif spam_metric == "tracedist": + # Leakage-aware metrics NOT available. + for preplabel, m_prep in mdl.preps.items(): + wt = item_weights.get(preplabel, spamWeight) + rhoMx1 = _tools.vec_to_stdmx(m_prep.to_dense(), mxBasis) + t_prep = target_model.preps[preplabel] + rhoMx2 = _tools.vec_to_stdmx(t_prep.to_dense(), mxBasis) + ret += wt * _tools.tracedist(rhoMx1, rhoMx2) + + for povmlabel in mdl.povms.keys(): + wt = item_weights.get(povmlabel, spamWeight) + ret += wt * _tools.povm_jtracedist(mdl, target_model, povmlabel) + + else: raise ValueError("Invalid spam_metric: %s" % spam_metric) + + return ret - else: - raise ValueError("Invalid metric: %s" % metric) - # end _objective_fn - _jacobian_fn = None return _objective_fn, _jacobian_fn From a351c3124ab8f5e7dc2ccc8e9b8e5a1bcf8be300 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Fri, 22 Nov 2024 09:14:49 -0500 Subject: [PATCH 50/51] improved HTML reporting when n_leak=1 --- pygsti/report/reportables.py | 13 +++++++++++++ pygsti/report/section/gauge.py | 2 +- pygsti/tools/optools.py | 13 +++++++++++++ 3 files changed, 27 insertions(+), 1 deletion(-) diff --git a/pygsti/report/reportables.py b/pygsti/report/reportables.py index b4c5eeed6..89c0d8f9a 100644 --- a/pygsti/report/reportables.py +++ b/pygsti/report/reportables.py @@ -1167,6 +1167,14 @@ def frobenius_diff(a, b, mx_basis): # assume vary model1, model2 fixed # init args == (model1, model2, op_label) +def leaky_gate_frob_dist(a, b, mx_basis): + n_leak = 1 + return _tools.subspace_restricted_fro_dist(a, b, mx_basis, n_leak) + + +Leaky_gate_frob_dist = _modf.opsfn_factory(leaky_gate_frob_dist) + + def jtrace_diff(a, b, mx_basis): # assume vary model1, model2 fixed """ Jamiolkowski trace distance between a and b @@ -2455,6 +2463,8 @@ def info_of_opfn_by_name(name): "where (a_i,b_i) are corresponding eigenvalues of A and B."), "frob": ("Frobenius|Distance", "sqrt( sum( (A_ij - B_ij)^2 ) )"), + "la-frob": ("Frobenius|Distance (subspace)", + "TO WRITE"), "unmodeled": ("Un-modeled|Error", "The per-operation budget used to account for un-modeled errors (model violation)") } @@ -2539,6 +2549,9 @@ def evaluate_opfn_by_name(name, model, target_model, op_label_or_string, elif name == "evnudiamond": fn = Eigenvalue_nonunitary_diamondnorm if b else \ Circuit_eigenvalue_nonunitary_diamondnorm + elif name == "la-frob": + assert b + fn = Leaky_gate_frob_dist elif name == "frob": fn = Fro_diff if b else \ Circuit_fro_diff diff --git a/pygsti/report/section/gauge.py b/pygsti/report/section/gauge.py index c9c642e81..5927f3d2f 100644 --- a/pygsti/report/section/gauge.py +++ b/pygsti/report/section/gauge.py @@ -171,7 +171,7 @@ def final_gates_vs_target_table_gauge_var(workspace, switchboard=None, confidenc if kwargs.get('n_leak', 0) == 0: display = ('inf', 'agi', 'trace', 'diamond', 'nuinf', 'nuagi') else: - display = ('inf', 'la-inf', 'agi', 'trace', 'la-trace', 'diamond', 'nuinf', 'nuagi') + display = ('la-frob', 'la-inf', 'la-trace', 'frob', 'inf', 'trace', 'agi', 'diamond', 'nuinf', 'nuagi') return workspace.GatesVsTargetTable( switchboard.mdl_final, switchboard.mdl_target, _cri(1, switchboard, confidence_level, ci_brevity), display=display diff --git a/pygsti/tools/optools.py b/pygsti/tools/optools.py index 9d956f18f..00820ec90 100644 --- a/pygsti/tools/optools.py +++ b/pygsti/tools/optools.py @@ -619,6 +619,19 @@ def leading_dxd_submatrix_basis_vectors(d: int, n: int, current_basis, return_la return submatrix_basis_vectors +def subspace_restricted_fro_dist(a, b, mx_basis, n_leak=0): + diff = a - b + if n_leak == 0: + return _np.linalg.norm(diff, 'fro') + if n_leak == 1: + d = int(_np.sqrt(a.shape[0])) + assert a.shape == b.shape == (d**2, d**2) + B = leading_dxd_submatrix_basis_vectors(d-n_leak, d, mx_basis) + P = B @ B.T.conj() + return _np.linalg.norm(diff @ P) + raise ValueError() + + def average_gate_fidelity(a, b, mx_basis='pp', is_tp=None, is_unitary=None): """ Computes the average gate fidelity (AGF) between two gates. From 0346cca1ac61b304a3d7e595befa0d84aea1e961 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Fri, 17 Jan 2025 11:18:41 -0500 Subject: [PATCH 51/51] remove EmbeddedBasis class that accidentally was kept during a rebase --- pygsti/baseobjs/basis.py | 280 +-------------------------------------- 1 file changed, 1 insertion(+), 279 deletions(-) diff --git a/pygsti/baseobjs/basis.py b/pygsti/baseobjs/basis.py index 87deebcbb..6b4c13137 100644 --- a/pygsti/baseobjs/basis.py +++ b/pygsti/baseobjs/basis.py @@ -1761,6 +1761,7 @@ def create_simple_equivalent(self, builtin_basis_name=None): builtin_basis_name = first_comp_name # if all components have the same name return BuiltinBasis(builtin_basis_name, self.elsize, sparse=self.sparse) + class OuterProdBasis(TensorProdBasis): def __init__(self, component_bases, name=None, longname=None, squeeze=True): @@ -1832,282 +1833,3 @@ def create_equivalent(self, builtin_basis_name): def create_simple_equivalent(self, builtin_basis_name=None): raise NotImplementedError() - -class EmbeddedBasis(LazyBasis): - """ - A basis that embeds a basis for a smaller state space within a larger state space. - - The elements of an EmbeddedBasis are therefore just embedded versions - of the elements of the basis that is embedded. - - Parameters - ---------- - basis_to_embed : Basis - The basis being embedded. - - state_space_labels : StateSpaceLabels - An object describing the struture of the entire state space. - - target_labels : list or tuple - The labels contained in `stateSpaceLabels` which demarcate the - portions of the state space acted on by `basis_to_embed`. - - name : str, optional - The name of this basis. If `None`, the names of `basis_to_embed` - is joined with ':' characters to the elements of `target_labels`. - - longname : str, optional - A longer description of this basis. If `None`, then a long name is - automatically generated. - """ - - @classmethod - def embed_label(cls, lbl, target_labels): - """ - Gets the EmbeddedBasis label for `lbl`. - - Convenience method that gives the EmbeddedBasis label for `lbl` - without needing to construct the `EmbeddedBasis`. E.g. `"XX:1,2"`. - - Parameters - ---------- - lbl : str - Un-embedded basis element label, e.g. `"XX"`. - - target_labels : tuple - The target state space labels upon which this basis element - will be embedded, e.g. `(1,2)` - - Returns - ------- - str - The embedded-basis-element label as an EmbeddedBasis would - assign it. E.g. `"XX:1,2"`. - """ - return "%s:%s" % (lbl, ",".join(map(str, target_labels))) - - @classmethod - def unembed_label(cls, lbl, target_labels): - """ - Convenience method that performs the reverse of :meth:`embed_label` - - Parameters - ---------- - lbl : str - Embedded basis element label, e.g. `"XX:1,2"`. - - target_labels : tuple - The target state space labels upon which this basis element - will be embedded, e.g. `(1,2)` - - Returns - ------- - str - The un-embedded label, e.g. `"XX"`. - """ - suffix = ":" + ",".join(map(str, target_labels)) - if lbl.endswith(suffix): - return lbl[:-len(suffix)] - else: - raise ValueError("Cannot unembed '%s' - doesn't end in '%s'!" % (lbl, suffix)) - - def __init__(self, basis_to_embed, state_space, target_labels, name=None, longname=None): - ''' - Create a new EmbeddedBasis. - - Parameters - ---------- - basis_to_embed : Basis - The basis being embedded. - - state_space : StateSpace - An object describing the struture of the entire state space. - - target_labels : list or tuple - The labels contained in `stateSpaceLabels` which demarcate the - portions of the state space acted on by `basis_to_embed`. - - name : str, optional - The name of this basis. If `None`, the names of `basis_to_embed` - is joined with ':' characters to the elements of `target_labels`. - - longname : str, optional - A longer description of this basis. If `None`, then a long name is - automatically generated. - ''' - from pygsti.baseobjs.statespace import StateSpace as _StateSpace - self.embedded_basis = basis_to_embed - self.target_labels = target_labels - self.state_space = _StateSpace.cast(state_space) - - if name is None: - name = ':'.join((basis_to_embed.name,) + tuple(map(str, target_labels))) - if longname is None: - longname = "Embedded %s basis as %s within %s" % \ - (basis_to_embed.name, ':'.join(map(str, target_labels)), str(self.state_space)) - - real = basis_to_embed.real - sparse = basis_to_embed.sparse - - super(EmbeddedBasis, self).__init__(name, longname, real, sparse) - - def _to_nice_serialization(self): - state = super()._to_nice_serialization() - state.update({'name': self.name, - 'longname': self.longname, - 'state_space': self.state_space.to_nice_serialization(), - 'embedded_basis': self.embedded_basis.to_nice_serialization() - }) - return state - - @classmethod - def _from_nice_serialization(cls, state): - basis_to_embed = Basis.from_nice_serialization(state['embedded_basis']) - state_space = _StateSpace.from_nice_serialization(state['state_space']) - return cls(basis_to_embed, state_space, state['target_labels'], state['name'], state['longname']) - - @property - def dim(self): - """ - The dimension of the vector space this basis fully or partially - spans. Equivalently, the length of the `vector_elements` of the - basis. - """ - return self.state_space.dim - - @property - def size(self): - """ - The number of elements (or vector-elements) in the basis. - """ - return self.embedded_basis.size - - @property - def elshape(self): - """ - The shape of each element. Typically either a length-1 or length-2 - tuple, corresponding to vector or matrix elements, respectively. - Note that *vector elements* always have shape `(dim,)` (or `(dim,1)` - in the sparse case). - """ - elndim = self.embedded_basis.elndim - if elndim == 2: # a "matrix" basis - d = int(_np.sqrt(self.dim)) - assert(d**2 == self.dim), \ - "Dimension of state_space must be a perfect square when embedding a matrix basis" - elshape = (d, d) - elif elndim == 1: - elshape = (self.dim,) - else: - raise ValueError("Can only embed bases with .elndim == 1 or 2 (received %d)!" % elndim) - return elshape - - def __hash__(self): - return hash(tuple(hash(self.embedded_basis), self.target_labels, self.state_space)) - - def _lazy_build_elements(self): - """ Take a dense or sparse basis matrix and embed it. """ - #LAZY building of elements (in case we never need them) - if self.elndim == 2: # then use EmbeddedOp to do matrix - from ..modelmembers.operations import StaticArbitraryOp - from ..modelmembers.operations import EmbeddedOp - sslbls = self.state_space.copy() - sslbls.reduce_dims_densitymx_to_state_inplace() # because we're working with basis matrices not gates - - if self.sparse: - self._elements = [] - for spmx in self.embedded_basis.elements: - mxAsOp = StaticArbitraryOp(spmx.to_dense(), evotype='statevec') - self._elements.append(EmbeddedOp(sslbls, self.target_labels, - mxAsOp).to_sparse()) - else: - self._elements = _np.zeros((self.size,) + self.elshape, 'complex') - for i, mx in enumerate(self.embedded_basis.elements): - self._elements[i] = EmbeddedOp( - sslbls, self.target_labels, StaticArbitraryOp(mx, evotype='statevec') - ).to_dense(on_space='HilbertSchmidt') - else: - # we need to perform embedding using vectors rather than matrices - doable, but - # not needed yet, so defer implementation to later. - raise NotImplementedError("Embedding *vector*-type bases not implemented yet") - - def _lazy_build_labels(self): - self._labels = [EmbeddedBasis.embed_label(lbl, self.target_labels) - for lbl in self.embedded_basis.labels] - - def _copy_with_toggled_sparsity(self): - return EmbeddedBasis(self.embedded_basis._copy_with_toggled_sparsity(), - self.state_space, - self.target_labels, - self.name, self.longname) - - def is_equivalent(self, other, sparseness_must_match=True): - """ - Tests whether this basis is equal to another basis, optionally ignoring sparseness. - - Parameters - ----------- - other : Basis or str - The basis to compare with. - - sparseness_must_match : bool, optional - If `False` then comparison ignores differing sparseness, and this function - returns `True` when the two bases are equal except for their `.sparse` values. - - Returns - ------- - bool - """ - otherIsBasis = isinstance(other, EmbeddedBasis) - if not otherIsBasis: return False # can't be equal to a non-EmbeddedBasis - if self.target_labels != other.target_labels or self.state_space != other.state_space: - return False - return self.embedded_basis.is_equivalent(other.embedded_basis, sparseness_must_match) - - def create_equivalent(self, builtin_basis_name): - """ - Create an equivalent basis with components of type `builtin_basis_name`. - - Create a Basis that is equivalent in structure & dimension to this - basis but whose simple components (perhaps just this basis itself) is - of the builtin basis type given by `builtin_basis_name`. - - Parameters - ---------- - builtin_basis_name : str - The name of a builtin basis, e.g. `"pp"`, `"gm"`, or `"std"`. Used to - construct the simple components of the returned basis. - - Returns - ------- - EmbeddedBasis - """ - equiv_embedded = self.embedded_basis.create_equivalent(builtin_basis_name) - return EmbeddedBasis(equiv_embedded, self.state_space, self.target_labels) - - def create_simple_equivalent(self, builtin_basis_name=None): - """ - Create a basis of type `builtin_basis_name` whose elements are compatible with this basis. - - Create a simple basis *and* one without components (e.g. a - :class:`TensorProdBasis`, is a simple basis w/components) of the - builtin type specified whose dimension is compatible with the - *elements* of this basis. This function might also be named - "element_equivalent", as it returns the `builtin_basis_name`-analogue - of the standard basis that this basis's elements are expressed in. - - Parameters - ---------- - builtin_basis_name : str, optional - The name of the built-in basis to use. If `None`, then a - copy of this basis is returned (if it's simple) or this - basis's name is used to try to construct a simple and - component-free version of the same builtin-basis type. - - Returns - ------- - Basis - """ - if builtin_basis_name is None: - builtin_basis_name = self.embedded_basis.name # default - return BuiltinBasis(builtin_basis_name, self.elsize, sparse=self.sparse)