From c260399d37542ade576cf2c9ed8a310fa3c1fb1a Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Wed, 3 Dec 2025 07:56:10 -0800 Subject: [PATCH 01/21] documentation for DenseOperator --- pygsti/modelmembers/operations/denseop.py | 65 +++++++++-------------- 1 file changed, 24 insertions(+), 41 deletions(-) diff --git a/pygsti/modelmembers/operations/denseop.py b/pygsti/modelmembers/operations/denseop.py index e5f600182..cce6e0088 100644 --- a/pygsti/modelmembers/operations/denseop.py +++ b/pygsti/modelmembers/operations/denseop.py @@ -255,7 +255,6 @@ def __complex__(self): return complex(self._ptr) class DenseOperator(DenseOperatorInterface, _KrausOperatorInterface, _LinearOperator): """ - TODO: update docstring An operator that behaves like a dense super-operator matrix. This class is the common base class for more specific dense operators. @@ -278,10 +277,29 @@ class DenseOperator(DenseOperatorInterface, _KrausOperatorInterface, _LinearOper The state space for this operation. If `None` a default state space with the appropriate number of qubits is used. - Attributes + Properties ---------- - base : numpy.ndarray - Direct access to the underlying process matrix data. + DenseOperatorInterface is deprecated and this class' dependence on it + will be removed in the immediate future. For the time being, the + implementation of __getattr__ in DenseOperatorInterface means that + DenseOperator has every property of its attached OpRepDenseSuperop, + which is stored in DenseOperator._rep. + + Private attributes + ------------------ + _evotype (required by ModelMember; positional arg in DenseOperator.__init__) + _rep (required by LinearOperator; either Cython or Python "OpRepDenseSuperop") + _basis (indirect requirement of _rep) + _ptr (required by DenseOperatorInterface, implemented as self._rep.base) + _state_space (required by ModelMember; inferrable from None as __init__ kwarg) + + Other private attributes + ------------------------ + ModelMember._gpindices + ModelMember._paramlbls + ModelMember._param_bounds + ModelMember._dirty + ModelMember._submember_rpindices """ @classmethod @@ -393,44 +411,9 @@ def _is_similar(self, other, rtol, atol): @property def kraus_operators(self): """A list of this operation's Kraus operators as numpy arrays.""" - # Let I index a basis element, rho be a d x d matrix, and (I // d, I mod d) := (i,ii) be the "2D" - # index corresponding to I. - # op(rho) = sum_IJ choi_IJ BI rho BJ_dag = sum_IJ (evecs_IK D_KK evecs_inv_KJ) BI rho BJ_dag - # Note: evecs can be and are assumed chosen to be orthonormal, so evecs_inv = evecs^dag - # if {Bi} is the set of matrix units then ... - # = sum_IJK(i',j') (evecs_IK D_KK evecs_inv_KJ) unitI_ii' rho_i'j' unitJ_j'j - # = sum_IJK(i',j') (evecs_(Ia,Ib)K D_KK evecs_inv_K(Ja,Jb)) unitI_ii' rho_i'j' unitJ_j'j - # using fact that sum(i') unitI_ii' ==> i'=Ib and delta_i,Ia factor - # = sum_IJK (evecs_(Ia,Ib)K D_KK evecs_inv_K(Ja,Jb)) rho_IbJa * delta_i,Ia, delta_j,Jb - # = sum_K D_KK [ sum_(Ib,Ja) evector[K]_(i,Ib) rho_IbJa dual_evector[K]_(Ja,j) ] - # -> let reshaped K-th evector be called O_K and dual O_K^dag - # = sum_K D_KK O_K rho O_K^dag assert(self._basis is not None), "Kraus operator functionality requires specifying a superoperator basis" - superop_mx = self.to_dense("HilbertSchmidt"); d = int(_np.round(_np.sqrt(superop_mx.shape[0]))) - std_basis = _Basis.cast('std', superop_mx.shape[0]) - choi_mx = _jt.jamiolkowski_iso(superop_mx, self._basis, std_basis) * d # see NOTE below - # NOTE: multiply by `d` (density mx dimension) to un-normalize choi_mx as given by - # jamiolkowski_iso. Here we *want* the trace of choi_mx to be `d`, not 1, so that - # op(rho) = sum_IJ choi_IJ BI rho BJ_dag is true. - - #CHECK 1 (to unit test?) REMOVE - #tmp_std = _bt.change_basis(superop_mx, self._basis, 'std') - #B = _bt.basis_matrices('std', superop_mx.shape[0]) - #check_superop = sum([ choi_mx[i,j] * _np.kron(B[i], B[j].conjugate()) for i in range(d*d) for j in range(d*d)]) - #assert(_np.allclose(check_superop, tmp_std)) - - evals, evecs = _np.linalg.eigh(choi_mx) - assert(_np.allclose(evecs @ _np.diag(evals) @ (evecs.conjugate().T), choi_mx)) - TOL = 1e-7 # consider lowering this tolerance as it leads to errors of this order in the Kraus decomp - if any([ev <= -TOL for ev in evals]): - raise ValueError("Cannot compute Kraus decomposition of non-positive-definite superoperator!") - kraus_ops = [evecs[:, i].reshape(d, d) * _np.sqrt(ev) for i, ev in enumerate(evals) if abs(ev) > TOL] - - #CHECK 2 (to unit test?) REMOVE - #std_superop = sum([_ot.unitary_to_std_process_mx(kop) for kop in kraus_ops]) - #assert(_np.allclose(std_superop, tmp_std)) - - return kraus_ops + kops = self._kraus_operators() + return kops def set_kraus_operators(self, kraus_operators): """ From e9bc2a39578fcca8e9a6c3bd7eecbfc1c0717600 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Wed, 3 Dec 2025 07:56:57 -0800 Subject: [PATCH 02/21] add the ability to specify custom germs in GSTModelPack.create_gst_experiment_design --- pygsti/modelpacks/_modelpack.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/pygsti/modelpacks/_modelpack.py b/pygsti/modelpacks/_modelpack.py index 358b7bf96..36bb98709 100644 --- a/pygsti/modelpacks/_modelpack.py +++ b/pygsti/modelpacks/_modelpack.py @@ -360,7 +360,8 @@ def create_gst_experiment_design(self, max_max_length, qubit_labels=None, fpr=Fa """ for k in kwargs.keys(): if k not in ('germ_length_limits', 'keep_fraction', 'keep_seed', 'include_lgst', 'nest', 'circuit_rules', - 'op_label_aliases', 'dscheck', 'action_if_missing', 'verbosity', 'add_default_protocol'): + 'op_label_aliases', 'dscheck', 'action_if_missing', 'verbosity', 'add_default_protocol', + 'germs'): raise ValueError("Invalid argument '%s' to StandardGSTDesign constructor" % k) if qubit_labels is None: qubit_labels = self._sslbls @@ -380,11 +381,13 @@ def create_gst_experiment_design(self, max_max_length, qubit_labels=None, fpr=Fa else: max_lengths_list = list(_gen_max_length(max_max_length)) + germs = kwargs.get('germs', self.germs(qubit_labels, lite)) + return _gst.StandardGSTDesign( self.processor_spec(qubit_labels), self.prep_fiducials(qubit_labels), self.meas_fiducials(qubit_labels), - self.germs(qubit_labels, lite), + germs, max_lengths_list, kwargs.get('germ_length_limits', None), fidpairs, From c5a288ffe194f4331ac744a9b0b05480f55eca93 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Wed, 3 Dec 2025 07:59:07 -0800 Subject: [PATCH 03/21] update matrixtools.eigenvalues and matrixtools.eigendecomposition to allow ndarrays where flags.writeable is False --- pygsti/tools/matrixtools.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/pygsti/tools/matrixtools.py b/pygsti/tools/matrixtools.py index a23025347..d43ac70f8 100644 --- a/pygsti/tools/matrixtools.py +++ b/pygsti/tools/matrixtools.py @@ -813,6 +813,8 @@ def eigenvalues(m: _np.ndarray, *, assume_hermitian: Optional[bool] = None, assu if assume_hermitian: # Make sure it's Hermtian in exact arithmetic. This helps with # reproducibility across different implementations of LAPACK. + if not m.flags.writeable: + m = m.copy() m += m.T.conj() m /= 2 return _np.linalg.eigvalsh(m) @@ -830,6 +832,8 @@ def eigendecomposition(m: _np.ndarray, *, assume_hermitian: Optional[bool] = Non if assume_hermitian: # Make sure it's Hermtian in exact arithmetic. This helps with # reproducibility across different implementations of LAPACK. + if not m.flags.writeable: + m = m.copy() m += m.T.conj() m /= 2 evals, evecs = _np.linalg.eigh(m) From 9d8e89395ff6288a46655e643bc1a0f3f2666750 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Wed, 3 Dec 2025 08:00:12 -0800 Subject: [PATCH 04/21] fix whitespace --- pygsti/models/model.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/pygsti/models/model.py b/pygsti/models/model.py index 562118a50..f3af1278a 100644 --- a/pygsti/models/model.py +++ b/pygsti/models/model.py @@ -1107,7 +1107,6 @@ def _get_shift(j): return _bisect.bisect_left(indices_to_remove, j) #rebuild the model index to model member map if needed. self._build_index_mm_map() - def _init_virtual_obj(self, obj): """ Initializes a "virtual object" - an object (e.g. LinearOperator) that *could* be a @@ -1232,8 +1231,6 @@ def set_parameter_value(self, index, val, close=False): """ self.set_parameter_values([index], [val], close) - - def set_parameter_values(self, indices, values, close=False): """ From a501a15a9c207be3c066ed070a2e987781c3d3b4 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Wed, 3 Dec 2025 08:03:45 -0800 Subject: [PATCH 05/21] optools.py: add minimal_kraus_decomposition function (important for this PR), rootconj_superop (important for this PR), plus partial_trace and trace_effect (nice to have but not strictly necessary for this PR). --- pygsti/tools/optools.py | 111 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 111 insertions(+) diff --git a/pygsti/tools/optools.py b/pygsti/tools/optools.py index 487f0e661..b77c81c8c 100644 --- a/pygsti/tools/optools.py +++ b/pygsti/tools/optools.py @@ -575,6 +575,117 @@ def tensorized_with_eye(op: _np.ndarray, op_basis: _Basis, ten_basis: Optional[_ return ten_op, ten_basis +def rootconj_superop(effect_superket: _np.ndarray, basis: _Basis, abstol_warn: float=1e-10, abstol_error: float=1e-8) -> _np.ndarray: + """ + Let E denote the Hermitian matrix representation effect_superket, where 0 ≤ E ≤ 1. + + This function returns the array representation (in `basis`) of the map that takes + a Hermitian matrix ρ to the Hermitian matrix E½ ρ E½. + """ + effect_mat = _bt.vec_to_stdmx(effect_superket, basis, keep_complex=True) + vecs, vals, inv_vecs = _mt.eigendecomposition(effect_mat, assume_hermitian=True) + + msg = f'Eigenvalues {vals} fall outside [0.0, 1.0], up to tolerance %s.' + if _np.any(vals < 0.0 - abstol_error) or _np.any(vals > 1.0 + abstol_error): + raise ValueError( msg % abstol_error ) + if _np.any(vals < 0.0 - abstol_warn) or _np.any(vals > 1.0 + abstol_warn ): + _warnings.warn( msg % abstol_warn ) + vals = _np.clip(vals, a_min=0.0, a_max=1.0) + + rooteffect_mat = (vecs * _np.sqrt(vals)[_np.newaxis, :]) @ inv_vecs + mx_std = _np.kron(rooteffect_mat, rooteffect_mat.T) + mx : _np.ndarray = _bt.change_basis(mx_std, 'std', basis, expect_real=True) # type: ignore + return mx + + +def partial_trace(mx: _np.ndarray, dims: tuple[int,...], axis: int) -> _np.ndarray: + """ + + Notes + ----- + This implementation is stolen from CVXPY, which was stolen from some Julia library. + """ + if mx.ndim < 2 or mx.shape[0] != mx.shape[1]: + raise ValueError("partial_trace only supports 2-d square arrays.") + if axis < 0 or axis >= len(dims): + msg = f"Invalid axis argument, should be between 0 and {len(dims)}, got {axis}." + raise ValueError(msg) + if mx.shape[0] != _np.prod(dims): + raise ValueError("Dimension of system doesn't correspond to dimension of subsystems.") + + def _term(j: int) -> _np.ndarray: + a = _sps.coo_matrix(([1.0], ([0], [0]))) + b = _sps.coo_matrix(([1.0], ([0], [0]))) + for (i_axis, dim) in enumerate(dims): + if i_axis == axis: + v = _sps.coo_matrix(([1], ([j], [0])), shape=(dim, 1)) + a = _sps.kron(a, v.T) + b = _sps.kron(b, v) + else: + eye_mat = _sps.eye_array(dim) + a = _sps.kron(a, eye_mat) + b = _sps.kron(b, eye_mat) + return a @ mx @ b + + return _np.sum([_term(j) for j in range(dims[axis])]) + + +def trace_effect(op: _np.ndarray, op_basis: BasisLike, on_space: SpaceT = 'HilbertSchmidt'): + """ + Let `op` be the array representation of a superoperator G in `op_basis`, + where G maps from and to the space of order-d Hermitian operators. + + The trace effect of G is the Heritian operator E that satifies + + trace(G(ρ)) = trace(E ρ) for all order-d Hermitian matrices ρ. + + If on_space='HilbertSchmidt', then this function returns a superket representation + of E in `op_basis`. If on_space='Hilbert', then we return E itself. + """ + d = int(_np.round(op.size ** 0.25)) + assert op.shape == (d**2, d**2) + basis = op_basis if isinstance(op_basis, _Basis) else _Basis.cast(op_basis, dim=d**2) + vecI = _bt.stdmx_to_vec(_np.eye(d), basis) + vecE = op.T.conj() @ vecI + if on_space == 'HilbertSchmidt': + return vecE + else: + E = _bt.vec_to_stdmx(vecE, op_basis) + return E + + +def minimal_kraus_decomposition(op_x: _np.ndarray, op_basis: _Basis, error_tol:float=1e-6, trunc_tol:float=1e-7) -> list[_np.ndarray]: + """ + The array `op_x` represents a completely positive superoperator X on + Hilbert-Schmidt space, using `op_basis` as the basis for that space. + + A Kraus decomposition of X is a set of square matrices, {K_i}_i, that satisfy + + X(ρ) = \\sum_i K_i ρ K_i^\\dagger. + + The matrices appearing in any such set are called Kraus operators of X. + + This function returns a minimal-length list of Kraus operators of X. + """ + d = int(_np.round(op_x.size ** 0.25)) + d2 = d**2 + assert op_x.shape == (d2, d2) + from pygsti.tools.jamiolkowski import jamiolkowski_iso + choi_mx : _np.ndarray = jamiolkowski_iso(op_x, op_basis, 'std', normalized=True) * d # type: ignore + + evecs, evals, _ = _mt.eigendecomposition(choi_mx, assume_hermitian=True) + if any([ev < -error_tol for ev in evals]): + raise ValueError("Cannot compute Kraus decomposition of non-positive-definite superoperator!") + keep = evals >= trunc_tol + evals = evals[keep] + evecs = evecs[:, keep] + out = [] + for i, ev in enumerate(evals): + temp = _np.sqrt(ev) * evecs[:, i].reshape((d, d), order='F') + out.append(temp) + return out + + 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 c1d2be0a532419538a93232202f25e0ea1995b64 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Wed, 3 Dec 2025 08:04:35 -0800 Subject: [PATCH 06/21] add LinearOperator._kraus_operators function with default implementation --- pygsti/modelmembers/operations/linearop.py | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/pygsti/modelmembers/operations/linearop.py b/pygsti/modelmembers/operations/linearop.py index 6edd84459..a1bc6695b 100644 --- a/pygsti/modelmembers/operations/linearop.py +++ b/pygsti/modelmembers/operations/linearop.py @@ -14,6 +14,7 @@ import numpy as _np +from pygsti.baseobjs.basis import Basis as _Basis from pygsti.baseobjs.opcalc import bulk_eval_compact_polynomials_complex as _bulk_eval_compact_polynomials_complex from pygsti.modelmembers import modelmember as _modelmember from pygsti.tools import optools as _ot @@ -103,6 +104,16 @@ def shape(self) -> tuple[int, int]: # are given broad freedom to define semantics of self._rep. return (self.dim, self.dim) + def _kraus_operators(self): + """A list of this operation's Kraus operators as numpy arrays.""" + basis = getattr(self, '_basis', None) + if not isinstance(basis, _Basis): + msg = "Kraus operator functionality requires specifying a superoperator basis" + raise NotImplementedError(msg) + mx = self.to_dense('HilbertSchmidt') + kops = _ot.minimal_kraus_decomposition(mx, basis, 1e-7, 1e-7) + return kops + def set_dense(self, m): """ Set the dense-matrix value of this operation. From 172a35e67b94031c543eff987f0abee51f266006 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Wed, 3 Dec 2025 08:09:01 -0800 Subject: [PATCH 07/21] address warning from not specifying default kwarg in np.linalg.lstsq --- pygsti/modelmembers/povms/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pygsti/modelmembers/povms/__init__.py b/pygsti/modelmembers/povms/__init__.py index d798b67a2..7dc36ebf6 100644 --- a/pygsti/modelmembers/povms/__init__.py +++ b/pygsti/modelmembers/povms/__init__.py @@ -612,7 +612,7 @@ def _objfn(v): tol=1e-13) if not soln.success and soln.fun > 1e-6: # not "or" because success is often not set correctly raise ValueError("Failed to find an errorgen such that Date: Wed, 3 Dec 2025 08:15:52 -0800 Subject: [PATCH 08/21] change Instrument.acton so it no longer requires a TP basis; add superket_trace to optools.py --- pygsti/modelmembers/instruments/instrument.py | 14 +++++++++----- pygsti/tools/optools.py | 9 +++++++++ 2 files changed, 18 insertions(+), 5 deletions(-) diff --git a/pygsti/modelmembers/instruments/instrument.py b/pygsti/modelmembers/instruments/instrument.py index 2e0a5a3a0..38e01ba5c 100644 --- a/pygsti/modelmembers/instruments/instrument.py +++ b/pygsti/modelmembers/instruments/instrument.py @@ -19,6 +19,9 @@ from pygsti.baseobjs import statespace as _statespace from pygsti.tools import matrixtools as _mt from pygsti.tools import slicetools as _slct +from pygsti.tools import basistools as _bt +from pygsti.tools import optools as _ot +from pygsti.baseobjs.basis import Basis as _Basis from pygsti.baseobjs.label import Label as _Label from pygsti.baseobjs.statespace import StateSpace as _StateSpace @@ -378,13 +381,14 @@ def acton(self, state): staterep = state._rep outcome_probs_and_states = _collections.OrderedDict() + for lbl, element in self.items(): output_rep = element._rep.acton(staterep) - output_unnormalized_state = output_rep.to_dense() - prob = output_unnormalized_state[0] * state.dim**0.25 - output_normalized_state = output_unnormalized_state / prob # so [0]th == 1/state_dim**0.25 - outcome_probs_and_states[lbl] = (prob, _state.StaticState(output_normalized_state, self.evotype, - self.state_space)) + unnormalized_state_array = output_rep.to_dense() + prob = _ot.superket_trace(unnormalized_state_array, element.basis) + output_state_array = unnormalized_state_array / prob + output_state = _state.StaticState(output_state_array, self.evotype, self.state_space) + outcome_probs_and_states[lbl] = (prob, output_state) return outcome_probs_and_states diff --git a/pygsti/tools/optools.py b/pygsti/tools/optools.py index b77c81c8c..84a8dcd2f 100644 --- a/pygsti/tools/optools.py +++ b/pygsti/tools/optools.py @@ -462,6 +462,15 @@ def is_trace_preserving(a: _np.ndarray, mx_basis: BasisLike='pp', tol=1e-8) -> b return check +def superket_trace(superket: _np.ndarray, basis: _Basis): + if basis.first_element_is_identity: + t = superket.ravel()[0] + else: + mx = _bt.vec_to_stdmx(superket, basis, keep_complex=True) + t = _np.real(_np.trace(mx)) + return t + + def entanglement_fidelity(a, b, mx_basis: BasisLike='pp', is_tp=None, is_unitary=None): """ Returns the "entanglement" process fidelity between gate matrices. From fe63e73948ab7c493293379511dfc96a3c164c1a Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Wed, 3 Dec 2025 08:17:33 -0800 Subject: [PATCH 09/21] CPTPLND instruments --- pygsti/modelmembers/instruments/__init__.py | 139 ++++++++++++---- .../instruments/cptpinstrument.py | 150 ++++++++++++++++++ 2 files changed, 260 insertions(+), 29 deletions(-) create mode 100644 pygsti/modelmembers/instruments/cptpinstrument.py diff --git a/pygsti/modelmembers/instruments/__init__.py b/pygsti/modelmembers/instruments/__init__.py index 8ef2c0ed2..c2f454a47 100644 --- a/pygsti/modelmembers/instruments/__init__.py +++ b/pygsti/modelmembers/instruments/__init__.py @@ -13,8 +13,17 @@ from .instrument import Instrument from .tpinstrument import TPInstrument from .tpinstrumentop import TPInstrumentOp - -from pygsti.tools import optools as _ot +from .cptpinstrument import RootConjOperator, SummedOperator + +import scipy.linalg as _la +import numpy as _np +from pygsti.tools import optools as _ot, basistools as _bt +from types import NoneType +from pygsti.baseobjs.label import Label +from pygsti.baseobjs.basis import Basis +from pygsti.modelmembers import operations as _op +from pygsti.modelmembers import povms as _pv +from pygsti.modelmembers.povms.basepovm import _BasePOVM # Avoid circular import import pygsti.modelmembers as _mm @@ -37,7 +46,6 @@ def instrument_type_from_op_type(op_type): # Limited set (only matching what is in convert) instr_conversion = { - 'auto': 'full', 'static unitary': 'static unitary', 'static clifford': 'static clifford', 'static': 'static', @@ -45,16 +53,12 @@ def instrument_type_from_op_type(op_type): 'full TP': 'full TP', 'full CPTP': 'full CPTP', 'full unitary': 'full unitary', + 'CPTPLND': 'CPTPLND' } instr_type_preferences = [] for typ in op_type_preferences: - instr_type = None - if _ot.is_valid_lindblad_paramtype(typ): - # Lindblad types are passed through as TP only (matching current convert logic) - instr_type = "full TP" - else: - instr_type = instr_conversion.get(typ, None) + instr_type = instr_conversion.get(typ, None) if instr_type is None: continue @@ -110,26 +114,103 @@ def convert(instrument, to_type, basis, ideal_instrument=None, flatten_structure The converted instrument, usually a distinct object from the object passed as input. """ - to_types = to_type if isinstance(to_type, (tuple, list)) else (to_type,) # HACK to support multiple to_type values + if not isinstance(to_type, str): + if len(to_type) > 1: + raise ValueError(f"Expected to_type to be a string, but got {to_type}") + to_type = to_type[0] + assert isinstance(to_type, str) + destination_types = {'full TP': TPInstrument} - NoneType = type(None) - - for to_type in to_types: - try: - if isinstance(instrument, destination_types.get(to_type, NoneType)): - return instrument - - if to_type == "full TP": - return TPInstrument(list(instrument.items()), instrument.evotype, instrument.state_space) - elif to_type in ("full", "static", "static unitary"): - from ..operations import convert as _op_convert - ideal_items = dict(ideal_instrument.items()) if (ideal_instrument is not None) else {} - members = [(k, _op_convert(g, to_type, basis, ideal_items.get(k, None), flatten_structure)) - for k, g in instrument.items()] - return Instrument(members, instrument.evotype, instrument.state_space) + + if isinstance(instrument, destination_types.get( to_type, NoneType ) ): + return instrument + + if to_type == "full TP": + inst_arrays = dict() + for k, v in instrument.items(): + if hasattr(v, 'to_dense'): + inst_arrays[k] = v.to_dense('HilbertSchmidt') else: - raise ValueError("Cannot convert an instrument to type %s" % to_type) - except: - pass # try next to_type + inst_arrays[k] = v + return TPInstrument(list(inst_arrays.items()), instrument.evotype, instrument.state_space) + + if to_type in ("full", "static", "static unitary"): + from ..operations import convert as _op_convert + ideal_items = dict(ideal_instrument.items()) if (ideal_instrument is not None) else {} + members = [(k, _op_convert(g, to_type, basis, ideal_items.get(k, None), flatten_structure)) + for k, g in instrument.items()] + return Instrument(members, instrument.evotype, instrument.state_space) + + if to_type == 'CPTPLND': + op_arrays = {k: v.to_dense('HilbertSchmidt') for (k,v) in instrument.items()} + inst = cptp_instrument(op_arrays, basis) + return inst + + raise ValueError("Cannot convert an instrument to type %s" % to_type) + + +def cptp_instrument(op_arrays: dict[str, _np.ndarray], basis: Basis, error_tol:float=1e-6, trunc_tol:float=1e-7, povm_errormap=None) -> Instrument: + unitaries = dict() + effects = dict() + + # Step 1. Build CPTPLND-parameterized unitaries and static POVM effects + # from Kraus decompositions of the provided operators. + for oplbl, op in op_arrays.items(): + krausops = _ot.minimal_kraus_decomposition(op, basis, error_tol, trunc_tol) + cur_effects = [] + cur_unitaries = [] + for K in krausops: + # Define F(rho) := K rho K^\\dagger. If we compute a polar decomposition + # K = u p, then we can express F(rho) = u (p rho p) u^\\dagger, which is + # a composition of a unitary channel and the "RootConjOperator" of p@p. + u, p = _la.polar(K, side='right') + u_linop = _op.StaticUnitaryOp(u, basis) # type: ignore + u_cptp = _op.convert(u_linop, 'CPTPLND', basis) + cur_unitaries.append(u_cptp) + E_superket = _bt.stdmx_to_vec(p @ p, basis) + E = _pv.StaticPOVMEffect(E_superket) + cur_effects.append(E) + effects[oplbl] = cur_effects + unitaries[oplbl] = cur_unitaries + + # Step 2. Build a CPTPLND-parameterized POVM from the static POVM effects. + # + # We can't use povms.convert(...) because we might have more + # effects than the Hilbert space dimension. + # + base_povm_effects = {} + for oplbl, cur_effects in effects.items(): + for i, E in enumerate(cur_effects): + elbl = Label((oplbl, i)) + base_povm_effects[elbl] = E + base_povm = _BasePOVM(base_povm_effects, preserve_sum=False) + udim = int(_np.round(basis.dim ** 0.5)) + if povm_errormap is None: + I_ideal = _op.StaticUnitaryOp(_np.eye(udim), basis) # type: ignore + I_cptplnd = _op.convert(I_ideal, 'CPTPLND', basis) + povm_errormap = I_cptplnd.factorops[1] + povm_cptp = _pv.ComposedPOVM(povm_errormap, base_povm) + inst_ops = dict() + + # Step 3. Assemble the CPTPLND-parameterized unitaries and POVM + # effects to represent each operator as a completely-positive + # trace-reducing map. + # + for lbl in op_arrays: + cur_effects = effects[lbl] + cur_unitaries = unitaries[lbl] + + op_summands = [] + for i, U_cptplnd in enumerate(cur_unitaries): + E_cptp = povm_cptp[Label((lbl, i))] + rcop = RootConjOperator(E_cptp, basis) + cptr = _op.ComposedOp([rcop, U_cptplnd]) + op_summands.append(cptr) + if len(op_summands) == 1: + op = op_summands[0] + else: + op = SummedOperator(op_summands, basis) - raise ValueError("Could not convert instrument to to type(s): %s" % str(to_types)) + inst_ops[lbl] = op + inst = Instrument(inst_ops) + return inst diff --git a/pygsti/modelmembers/instruments/cptpinstrument.py b/pygsti/modelmembers/instruments/cptpinstrument.py new file mode 100644 index 000000000..ba6ea5d7c --- /dev/null +++ b/pygsti/modelmembers/instruments/cptpinstrument.py @@ -0,0 +1,150 @@ +""" +Defines the CPTPInstrument class +""" +#*************************************************************************************************** +# Copyright 2015, 2019, 2025 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. +#*************************************************************************************************** + +import numpy as _np +from pygsti.pgtypes import SpaceT +from pygsti.baseobjs.basis import Basis +from pygsti.modelmembers.povms.effect import POVMEffect +from pygsti.modelmembers.operations import LinearOperator +from pygsti.tools import optools as _ot + +from typing import Union +BasisLike = Union[Basis, str] + + +class RootConjOperator(LinearOperator): + """ + A linear operator parameterized by a matrix E where 0 ≤ E ≤ 1, whose action takes ρ to E½ ρ E½ . + + Every CPTR map can be obtained by pre-composing a CPTP map with this kind of linear operator. + """ + + EIGTOL_WARNING = 1e-10 + EIGTOL_ERROR = 1e-8 + + def __init__(self, effect: POVMEffect, basis: BasisLike): + self._basis = Basis.cast(basis, effect.dim) + self._effect = effect + self._state_space = effect.state_space + self._evotype = effect.evotype + + dim = self._state_space.dim + self._rep = self._evotype.create_dense_superop_rep( _np.zeros((dim, dim)), self._basis, self._state_space ) + self._update_rep_base() + LinearOperator.__init__(self, self._rep, self._evotype) + self.init_gpindices() + + def submembers(self): + return [self._effect] + + def _update_rep_base(self): + # This function is directly analogous to TPInstrumentOp._construct_matrix. + self._rep.base.flags.writeable = True + assert(self._rep.base.shape == (self.dim, self.dim)) + effect_superket = self._effect.to_dense() + mx = _ot.rootconj_superop(effect_superket, self._basis) + self._rep.base[:] = mx + self._rep.base.flags.writeable = False + self._rep.base_has_changed() + return + + def deriv_wrt_params(self, wrt_filter=None): + raise NotImplementedError() + + def has_nonzero_hessian(self): + # This is not affine in its parameters. + return True + + def from_vector(self, v, close=False, dirty_value=True): + for sm, local_inds in zip(self.submembers(), self._submember_rpindices): + sm.from_vector(v[local_inds], close, dirty_value) + self._update_rep_base() + return + + @property + def num_params(self): + return len(self.gpindices_as_array()) + + def to_vector(self): + v = _np.empty(self.num_params, 'd') + for param_op, local_inds in zip(self.submembers(), self._submember_rpindices): + v[local_inds] = param_op.to_vector() + return v + + def to_dense(self, on_space: SpaceT = 'HilbertSchmidt') -> _np.ndarray: + assert on_space in ('HilbertSchmidt', 'minimal') + out = self._rep.base.copy() + out.flags.writeable = True + return out + + +class SummedOperator(LinearOperator): + + + def __init__(self, operators, basis: BasisLike): + op = operators[0] + self._basis = Basis.cast(basis, op.dim) + self._operators = operators + self._state_space = op.state_space + self._evotype = op.evotype + self._subreps = [op._rep for op in self._operators] + self._rep = self._evotype.create_sum_rep( self._subreps, self._state_space ) + LinearOperator.__init__(self, self._rep, self._evotype) + self.init_gpindices() + # NOTE: This class doesn't have a function analogous to _update_rep_base + # that we use in RootConjOperator. We can get away with not having such + # a function because it's the responsibility of op.from_vector(...) + # to update op's attached OpRep. + return + + def submembers(self): + out = [] + hit = set() + for op in self._operators: + temp = op.submembers() + for sm in temp: + if id(temp) not in hit: + hit.add(id(temp)) + out.append(sm) + return out + + def deriv_wrt_params(self, wrt_filter=None): + raise NotImplementedError() + + def has_nonzero_hessian(self): + return any(op.has_nonzero_hession() for op in self._operators) + + def from_vector(self, v, close=False, dirty_value=True): + for sm, local_inds in zip(self.submembers(), self._submember_rpindices): + sm.from_vector(v[local_inds], close, dirty_value) + return + + @property + def num_params(self): + return len(self.gpindices_as_array()) + + def to_vector(self): + v = _np.empty(self.num_params, 'd') + for param_op, local_inds in zip(self.submembers(), self._submember_rpindices): + v[local_inds] = param_op.to_vector() + return v + + def to_dense(self, on_space: SpaceT = 'HilbertSchmidt') -> _np.ndarray: + assert on_space in ('HilbertSchmidt', 'minimal') + on_space = 'HilbertSchmidt' + out = self._operators[0].to_dense(on_space) + if not out.flags.writeable: + out = out.copy() + for op in self._operators[1:]: + temp = op.to_dense(on_space) + out += temp + return out From ab3d7f3a8598bd9c1bf9373f8db2f1f2b5097110 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Wed, 3 Dec 2025 08:49:37 -0800 Subject: [PATCH 10/21] update Instruments.ipynb advanced example notebook --- .../objects/advanced/Instruments.ipynb | 335 ++++++++++-------- pygsti/modelmembers/instruments/__init__.py | 2 +- 2 files changed, 197 insertions(+), 140 deletions(-) diff --git a/jupyter_notebooks/Tutorials/objects/advanced/Instruments.ipynb b/jupyter_notebooks/Tutorials/objects/advanced/Instruments.ipynb index 3e4093b39..989bf4acb 100644 --- a/jupyter_notebooks/Tutorials/objects/advanced/Instruments.ipynb +++ b/jupyter_notebooks/Tutorials/objects/advanced/Instruments.ipynb @@ -12,13 +12,15 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "import pygsti\n", "from pygsti.modelpacks import smq1Q_XYI as std\n", - "import numpy as np" + "import numpy as np\n", + "from pygsti.modelmembers.instruments import Instrument\n", + "from pprint import pprint" ] }, { @@ -31,58 +33,40 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "#Make a copy so we don't modify the original\n", - "target_model = std.target_model()\n", + "mdl_ideal = std.target_model()\n", "\n", - "#Create and add the ideal instrument\n", - "E0 = target_model.effects['0']\n", - "E1 = target_model.effects['1']\n", - " # Alternate indexing that uses POVM label explicitly\n", - " # E0 = target_model['Mdefault']['0'] # 'Mdefault' = POVM label, '0' = effect label\n", - " # E1 = target_model['Mdefault']['1']\n", - "Gmz_plus = np.dot(E0,E0.T) #note effect vectors are stored as column vectors\n", + "# Create and add the ideal instrument\n", + "E0 = mdl_ideal.effects['0']\n", + "E1 = mdl_ideal.effects['1']\n", + "Gmz_plus = np.dot(E0,E0.T) # note effect vectors are stored as column vectors\n", "Gmz_minus = np.dot(E1,E1.T)\n", - "target_model[('Iz',0)] = pygsti.modelmembers.instruments.Instrument({'p0': Gmz_plus, 'p1': Gmz_minus})\n", - "\n", - "#For later use, record the identity POVM vector\n", - "povm_ident = target_model.effects['0'] + target_model.effects['1'] " + "inst = Instrument({'p0': Gmz_plus, 'p1': Gmz_minus})\n", + "mdl_ideal[('Iz',0)] = inst" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "In order to generate some simulated data later on, we'll now create a noisy version of `target_model` by depolarizing the state preparation, gates, and POVM, and also rotating the basis that is measured by the instrument and POVM." + "In order to generate some simulated data later on, we'll now create a noisy version of `mdl_ideal` by depolarizing the state preparation, gates, and POVM, and also rotating the basis that is measured by the instrument and POVM." ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 14, "metadata": {}, "outputs": [], "source": [ - "mdl_noisy = target_model.depolarize(op_noise=0.01, spam_noise=0.01)\n", - "mdl_noisy.effects.depolarize(0.01) #because above call only depolarizes the state prep, not the POVMs\n", - "\n", - "# add a rotation error to the POVM\n", - "Uerr = pygsti.rotation_gate_mx([0,0.02,0])\n", - "E0 = np.dot(mdl_noisy['Mdefault']['0'].T,Uerr).T\n", - "E1 = povm_ident - E0\n", - "mdl_noisy.povms['Mdefault'] = pygsti.modelmembers.povms.UnconstrainedPOVM({'0': E0, '1': E1},\n", - " evotype='default')\n", - "\n", - "# Use the same rotated effect vectors to \"rotate\" the instrument Iz too\n", - "E0 = mdl_noisy.effects['0']\n", - "E1 = mdl_noisy.effects['1']\n", - "Gmz_plus = np.dot(E0,E0.T)\n", - "Gmz_minus = np.dot(E1,E1.T)\n", - "mdl_noisy[('Iz',0)] = pygsti.modelmembers.instruments.Instrument({'p0': Gmz_plus, 'p1': Gmz_minus})\n", - "\n", - "#print(mdl_noisy) #print the model" + "mdl_noisy = mdl_ideal.depolarize(op_noise=0.005, spam_noise=0.01)\n", + "mdl_noisy = mdl_noisy.rotate(max_rotate=0.025, seed=2048)\n", + "mdl_noisy.effects.depolarize(0.01)\n", + "mdl_ideal.convert_members_inplace('CPTPLND')\n", + "mdl_noisy.convert_members_inplace('CPTPLND')" ] }, { @@ -95,20 +79,66 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 15, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "{('p0', '0'): 0.5000000000000003,\n", + " ('p0', '1'): 1.4597476004180728e-17,\n", + " ('p1', '0'): 5.073268178784763e-18,\n", + " ('p1', '1'): 0.5}\n", + "\n", + "{('p0', '0'): 0.481367923790963,\n", + " ('p0', '1'): 0.0024189342904068554,\n", + " ('p1', '0'): 0.0025810657095932145,\n", + " ('p1', '1'): 0.5136320762090365}\n", + "\n" + ] + } + ], "source": [ - "dict(target_model.probabilities( pygsti.circuits.Circuit(( ('Gxpi2',0) , ('Iz',0) )) ))" + "for mdl in [mdl_ideal, mdl_noisy]:\n", + " pprint(dict(mdl.probabilities( pygsti.circuits.Circuit(( ('Gxpi2',0) , ('Iz',0) )) )))\n", + " print()" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 16, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "{('p0', 'p0', '0'): 0.5000000000000002,\n", + " ('p0', 'p0', '1'): -1.6613637599827203e-18,\n", + " ('p0', 'p1', '0'): -1.1185571585378685e-17,\n", + " ('p0', 'p1', '1'): 0.4999999999999999,\n", + " ('p1', 'p0', '0'): 1.625883976416346e-17,\n", + " ('p1', 'p0', '1'): 1.8939874217871704e-34,\n", + " ('p1', 'p1', '0'): 4.5374861265545935e-34,\n", + " ('p1', 'p1', '1'): 1.625883976416347e-17}\n", + "\n", + "{('p0', 'p0', '0'): 0.4775917799207901,\n", + " ('p0', 'p0', '1'): 0.002399958693069308,\n", + " ('p0', 'p1', '0'): 0.0025624981205110728,\n", + " ('p0', 'p1', '1'): 0.5099371259816915,\n", + " ('p1', 'p0', '0'): 0.0038579004920827257,\n", + " ('p1', 'p0', '1'): 1.938643463358172e-05,\n", + " ('p1', 'p1', '0'): 1.8156751786106926e-05,\n", + " ('p1', 'p1', '1'): 0.003613193605435256}\n", + "\n" + ] + } + ], "source": [ - "dict(target_model.probabilities( pygsti.circuits.Circuit(( ('Iz',0), ('Gxpi2',0) , ('Iz',0) )) ))" + "for mdl in [mdl_ideal, mdl_noisy]:\n", + " pprint(dict(mdl.probabilities( pygsti.circuits.Circuit(( ('Iz',0), ('Gxpi2',0) , ('Iz',0) )) )))\n", + " print()" ] }, { @@ -120,11 +150,21 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 17, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "probs = {('0',): 0.5000000000000004, ('1',): 0.5}\n", + "probs['0'] = 0.5000000000000004\n", + "probs[('0',)] = 0.5000000000000004\n" + ] + } + ], "source": [ - "probs = target_model.probabilities( pygsti.circuits.Circuit([('Gxpi2',0)]) )\n", + "probs = mdl_ideal.probabilities( pygsti.circuits.Circuit([('Gxpi2',0)]) )\n", "print(\"probs = \",dict(probs))\n", "print(\"probs['0'] = \", probs['0']) #This works...\n", "print(\"probs[('0',)] = \", probs[('0',)]) # and so does this." @@ -135,158 +175,175 @@ "metadata": {}, "source": [ "## Performing tomography\n", - "\n", - "### Simulated data generation\n", - "Now let's perform tomography on a model that includes instruments. First, we'll generate some data using `mdl_noisy` in exactly the same way as we would for any other model:" + "Now let's perform tomography on a model that includes instruments. \n", + "First, we'll build an experiment design. This notebook's experiment design makes a minor modification to the default design for our working modelpack (smq1Q_XYI). " ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 18, "metadata": {}, "outputs": [], "source": [ "germs = std.germs()\n", - "germs += [pygsti.circuits.Circuit([('Iz', 0)])] # add the instrument as a germ.\n", - "\n", - "prep_fiducials = std.prep_fiducials()\n", - "meas_fiducials = std.meas_fiducials()\n", - "max_lengths = [1] # keep it simple & fast\n", - "\n", - "lsgst_list = pygsti.circuits.create_lsgst_circuits(\n", - " mdl_noisy,prep_fiducials,meas_fiducials,germs,max_lengths)\n", - "\n", - "#print(\"LinearOperator sequences:\")\n", - "#print(lsgst_list) #note that this contains LGST strings with \"Iz\"\n", - "\n", - "#Create the DataSet\n", - "ds = pygsti.data.simulate_data(mdl_noisy,lsgst_list,1000,'multinomial',seed=2018)\n", - "\n", - "#Write it to a text file to demonstrate the format:\n", - "pygsti.io.write_dataset(\"../../tutorial_files/intermediate_meas_dataset.txt\",ds)" + "germs += [pygsti.circuits.Circuit([('Iz', 0)])]\n", + "ed = std.create_gst_experiment_design(max_max_length=8, germs=germs)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Notice the format of [intermediate_meas_dataset.txt](../../tutorial_files/intermediate_meas_dataset.txt), which includes a column for each distinct outcome tuple. Since not all experiments contain data for all outcome tuples, the `\"--\"` is used as a placeholder. Now that the data is generated, we run LGST or LSGST just like we would for any other model:\n", - "\n", - "### LGST" + "### Simulated data generation\n", + "Next, we generate data using `mdl_noisy` in exactly the same way as we would for any other model. We write the data to disk so you can see how datasets look when they include measurement data." ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 19, "metadata": {}, "outputs": [], "source": [ - "#Run LGST\n", - "mdl_lgst = pygsti.run_lgst(ds, prep_fiducials, meas_fiducials, target_model)\n", - "#print(mdl_lgst)\n", - "\n", - "#Gauge optimize the result to the true data-generating model (mdl_noisy),\n", - "# and compare. Mismatch is due to finite sample noise.\n", - "mdl_lgst_opt = pygsti.gaugeopt_to_target(mdl_lgst,mdl_noisy)\n", - "print(mdl_noisy.strdiff(mdl_lgst_opt))\n", - "print(\"Frobdiff after GOpt = \",mdl_noisy.frobeniusdist(mdl_lgst_opt))" + "ds = pygsti.data.simulate_data(mdl_noisy, ed.all_circuits_needing_data, 10_000, 'multinomial', seed=2018)\n", + "pygsti.io.write_dataset(\"../../tutorial_files/intermediate_meas_dataset.txt\", ds)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### Long-sequence GST\n", - "Instruments just add parameters to a `Model` like gates, state preparations, and POVMs do. The total number of parameters in our model is \n", + "Notice the format of [intermediate_meas_dataset.txt](../../tutorial_files/intermediate_meas_dataset.txt), which includes a column for each distinct outcome tuple. Since not all experiments contain data for all outcome tuples, the `\"--\"` is used as a placeholder. Now that the data is generated, we run LGST or LSGST just like we would for any other model:\n", "\n", - "$4$ (prep) + $2\\times 4$ (2 effects) + $5\\times 16$ (3 gates and 2 instrument members) $ = 92$." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "target_model.num_params" + "### Running the GST fit and generating a report" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 20, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "-- Std Practice: Iter 1 of 2 (full TP) --: \n", + " --- Iterative GST: [##################################################] 100.0% 592 circuits ---\n", + "-- Std Practice: Iter 2 of 2 (CPTPLND) --: \n", + " --- Iterative GST: [##################################################] 100.0% 592 circuits ---\n" + ] + } + ], "source": [ - "#Run long sequence GST\n", - "results = pygsti.run_long_sequence_gst(ds,target_model,prep_fiducials,meas_fiducials,germs,max_lengths)" + "from pygsti.protocols import StandardGST, ProtocolData\n", + "gst = StandardGST(\n", + " modes=('full TP', 'CPTPLND'), target_model=mdl_ideal, verbosity=2\n", + ")\n", + "pd = ProtocolData(ed, ds)\n", + "res = gst.run(pd)" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 21, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Running idle tomography\n", + "Computing switchable properties\n", + "Found standard clifford compilation from smq1Q_XYI\n", + "Found standard clifford compilation from smq1Q_XYI\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/rjmurr/Documents/pygsti-general/pyGSTi/pygsti/models/model.py:145: UserWarning:\n", + "\n", + "Model.num_modeltest_params could not obtain number of *non-gauge* parameters - using total instead\n", + "\n", + "/Users/rjmurr/Documents/pygsti-general/pyGSTi/pygsti/tools/optools.py:152: UserWarning:\n", + "\n", + "The input matrix a is not trace-1 up to tolerance 1.4901161193847656e-08. Beware result!\n", + "\n", + "/Users/rjmurr/Documents/pygsti-general/pyGSTi/pygsti/tools/optools.py:154: UserWarning:\n", + "\n", + "The input matrix b is not trace-1 up to tolerance 1.4901161193847656e-08. Beware result!\n", + "\n", + "/Users/rjmurr/Documents/pygsti-general/pyGSTi/pygsti/tools/optools.py:179: UserWarning:\n", + "\n", + "\n", + " Input matrix is not PSD up to tolerance 1.4901161193847656e-08.\n", + " We'll project out the bad eigenspaces to only work with the PSD part.\n", + " \n", + "\n", + "/Users/rjmurr/Documents/pygsti-general/pyGSTi/pygsti/report/workspacetables.py:2741: RuntimeWarning:\n", + "\n", + "divide by zero encountered in log\n", + "\n", + "/Users/rjmurr/Documents/pygsti-general/pyGSTi/pygsti/forwardsims/mapforwardsim.py:732: UserWarning:\n", + "\n", + "Generating dense process matrix representations of circuits or gates \n", + "can be inefficient and should be avoided for the purposes of forward \n", + "simulation/calculation of circuit outcome probability distributions \n", + "when using the MapForwardSimulator.\n", + "\n" + ] + } + ], "source": [ - "#Compare estimated model (after gauge opt) to data-generating one\n", - "mdl_est = results.estimates['GateSetTomography'].models['go0']\n", - "mdl_est_opt = pygsti.gaugeopt_to_target(mdl_est,mdl_noisy)\n", - "print(\"Frobdiff after GOpt = \", mdl_noisy.frobeniusdist(mdl_est_opt))" + "from pygsti.report import construct_standard_report\n", + "report_dir = '../../tutorial_files/cptp-instrument-report'\n", + "report_object = construct_standard_report(res, title='CPTP intrument GST')\n", + "report_object.write_html(report_dir)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "The same analysis can be done for a trace-preserving model, whose instruments are constrained to *add* to a perfectly trace-preserving map. The number of parameters in the model are now: \n", + "**Thats it!** You've done tomography on a model with intermediate measurments (instruments).\n", "\n", - "$3$ (prep) + $1\\times 4$ (effect and complement) + $3\\times 12$ (3 gates) + $(2\\times 16 - 3)$ (TP instrument) $ = 71$" + "As a bonus, the code below checks for violation of complete positivity of the instrument operations from the two fits. We see that the CPTP fit has no violation, while the full TP fit has small violations." ] }, { "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "mdl_targetTP = target_model.copy()\n", - "mdl_targetTP.set_all_parameterizations(\"full TP\")\n", - "print(\"POVM type = \",type(mdl_targetTP[\"Mdefault\"]),\" Np=\",mdl_targetTP[\"Mdefault\"].num_params)\n", - "print(\"Instrument type = \",type(mdl_targetTP[(\"Iz\",0)]),\" Np=\",mdl_targetTP[(\"Iz\",0)].num_params)\n", - "print(\"Number of model parameters = \", mdl_targetTP.num_params)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "resultsTP = pygsti.run_long_sequence_gst(ds,mdl_targetTP,prep_fiducials,meas_fiducials,germs,max_lengths)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "#Again compare estimated model (after gauge opt) to data-generating one\n", - "mdl_est = resultsTP.estimates['GateSetTomography'].models['go0']\n", - "mdl_est_opt = pygsti.gaugeopt_to_target(mdl_est,mdl_noisy)\n", - "print(\"Frobdiff after GOpt = \", mdl_noisy.frobeniusdist(mdl_est_opt))" - ] - }, - { - "cell_type": "markdown", + "execution_count": 22, "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Iz:0\n", + "p0: 0.0009697360237423556\n", + "p1: 2.6230206864367674e-05\n", + "\n", + "Iz:0\n", + "p0: 0\n", + "p1: 0\n", + "\n" + ] + } + ], "source": [ - "**Thats it!** You've done tomography on a model with intermediate measurments (instruments)." + "for mdl in [res.estimates['full TP'].models['stdgaugeopt'], res.estimates['CPTPLND'].models['stdgaugeopt']]:\n", + " for instlbl, inst in mdl.instruments.items():\n", + " print(instlbl)\n", + " for ioplbl, iop in inst.items():\n", + " violation = max(0, pygsti.tools.sum_of_negative_choi_eigenvalues_gate(iop.to_dense(), 'pp'))\n", + " print(ioplbl + ': ' + str(violation) )\n", + " print()" ] } ], "metadata": { "kernelspec": { - "display_name": "Python 3 (ipykernel)", + "display_name": "pgdev311", "language": "python", "name": "python3" }, @@ -300,7 +357,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.9.13" + "version": "3.11.13" } }, "nbformat": 4, diff --git a/pygsti/modelmembers/instruments/__init__.py b/pygsti/modelmembers/instruments/__init__.py index c2f454a47..9864f0987 100644 --- a/pygsti/modelmembers/instruments/__init__.py +++ b/pygsti/modelmembers/instruments/__init__.py @@ -190,12 +190,12 @@ def cptp_instrument(op_arrays: dict[str, _np.ndarray], basis: Basis, error_tol:f I_cptplnd = _op.convert(I_ideal, 'CPTPLND', basis) povm_errormap = I_cptplnd.factorops[1] povm_cptp = _pv.ComposedPOVM(povm_errormap, base_povm) - inst_ops = dict() # Step 3. Assemble the CPTPLND-parameterized unitaries and POVM # effects to represent each operator as a completely-positive # trace-reducing map. # + inst_ops = dict() for lbl in op_arrays: cur_effects = effects[lbl] cur_unitaries = unitaries[lbl] From 14bf71a759405b24a6ea696d555b3a51b79e4dc4 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Wed, 17 Dec 2025 08:26:39 -0800 Subject: [PATCH 11/21] drop unnecessary NotImplementedError in favor of acceptable base class implementation --- pygsti/modelmembers/instruments/cptpinstrument.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pygsti/modelmembers/instruments/cptpinstrument.py b/pygsti/modelmembers/instruments/cptpinstrument.py index ba6ea5d7c..7f7988c54 100644 --- a/pygsti/modelmembers/instruments/cptpinstrument.py +++ b/pygsti/modelmembers/instruments/cptpinstrument.py @@ -58,7 +58,7 @@ def _update_rep_base(self): return def deriv_wrt_params(self, wrt_filter=None): - raise NotImplementedError() + return LinearOperator.deriv_wrt_params(self, wrt_filter) def has_nonzero_hessian(self): # This is not affine in its parameters. From b5177bb5d31ddebffa2d9ba76703a79c5db873a3 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Thu, 29 Jan 2026 09:48:45 -0800 Subject: [PATCH 12/21] change to allow any non-CPTPLND Lindblad parameterization. If anything other than CPTPLND or GLND we raise a warning. --- pygsti/modelmembers/instruments/__init__.py | 19 +++++++++++++++++-- 1 file changed, 17 insertions(+), 2 deletions(-) diff --git a/pygsti/modelmembers/instruments/__init__.py b/pygsti/modelmembers/instruments/__init__.py index 9864f0987..edbe6f6ff 100644 --- a/pygsti/modelmembers/instruments/__init__.py +++ b/pygsti/modelmembers/instruments/__init__.py @@ -15,6 +15,7 @@ from .tpinstrumentop import TPInstrumentOp from .cptpinstrument import RootConjOperator, SummedOperator +import warnings as _warnings import scipy.linalg as _la import numpy as _np from pygsti.tools import optools as _ot, basistools as _bt @@ -53,6 +54,9 @@ def instrument_type_from_op_type(op_type): 'full TP': 'full TP', 'full CPTP': 'full CPTP', 'full unitary': 'full unitary', + 'GLND': 'full TP', + # ^ It's pretty harmless to associate GLND operations with "full TP" + # instruments. In both cases we're relaxing positivity constraints. 'CPTPLND': 'CPTPLND' } @@ -60,8 +64,19 @@ def instrument_type_from_op_type(op_type): for typ in op_type_preferences: instr_type = instr_conversion.get(typ, None) - if instr_type is None: - continue + if instr_type is None and _ot.is_valid_lindblad_paramtype(typ): + # NOTE: need to update the message below if more lindblad + # types are added as keys to the instr_conversion dict. + msg = \ + f""" + Operation type {typ} is a Lindblad parameterization, but + is neither 'GLND' or 'CPTPLND'. That means you might be + asking for a reduced-order model that's a subset of all + CPTP models. We don't support that parameterization at this + time. We're falling back to a 'full TP' parameterization! + """ + _warnings.warn(msg) + instr_type = 'full TP' # non-CPTPLND falls back to full TP. if instr_type not in instr_type_preferences: instr_type_preferences.append(instr_type) From ceba6ed648953efaf8978e36f8c82bb23c96f008 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Thu, 29 Jan 2026 09:50:33 -0800 Subject: [PATCH 13/21] add back in explicit auto handling --- pygsti/modelmembers/instruments/__init__.py | 1 + 1 file changed, 1 insertion(+) diff --git a/pygsti/modelmembers/instruments/__init__.py b/pygsti/modelmembers/instruments/__init__.py index edbe6f6ff..ff7b2ae44 100644 --- a/pygsti/modelmembers/instruments/__init__.py +++ b/pygsti/modelmembers/instruments/__init__.py @@ -47,6 +47,7 @@ def instrument_type_from_op_type(op_type): # Limited set (only matching what is in convert) instr_conversion = { + 'auto': 'full TP', 'static unitary': 'static unitary', 'static clifford': 'static clifford', 'static': 'static', From d5140cb8b9b74d531b12c6855cd5b90b9dfba57f Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Thu, 29 Jan 2026 18:31:16 -0800 Subject: [PATCH 14/21] remove redundant function that raised NotImplementedError. Formatting. --- pygsti/modelmembers/states/cptpstate.py | 33 +++---------------------- 1 file changed, 4 insertions(+), 29 deletions(-) diff --git a/pygsti/modelmembers/states/cptpstate.py b/pygsti/modelmembers/states/cptpstate.py index 79357efe5..9ee09f372 100644 --- a/pygsti/modelmembers/states/cptpstate.py +++ b/pygsti/modelmembers/states/cptpstate.py @@ -315,10 +315,10 @@ def deriv_wrt_params(self, wrt_filter=None): # Derivative of vector wrt params; shape == [vecLen,dmDim,dmDim] *not dealing with TP condition yet* # (first get derivative assuming last diagonal el of Lmx *is* a parameter, then use chain rule) - dVdp = _np.einsum('aml,mb,ab->lab', conj_basis_mxs, Lbar, F1) # only a >= b nonzero (F1) - dVdp += _np.einsum('mal,mb,ab->lab', conj_basis_mxs, L, F1) # ditto - dVdp += _np.einsum('bml,ma,ab->lab', conj_basis_mxs, Lbar, F2) # only b > a nonzero (F2) - dVdp += _np.einsum('mbl,ma,ab->lab', conj_basis_mxs, L, F2.conjugate()) # ditto + dVdp = _np.einsum('aml,mb,ab->lab', conj_basis_mxs, Lbar, F1) # only a >= b nonzero (F1) + dVdp += _np.einsum('mal,mb,ab->lab', conj_basis_mxs, L, F1) # ditto + dVdp += _np.einsum('bml,ma,ab->lab', conj_basis_mxs, Lbar, F2) # only b > a nonzero (F2) + dVdp += _np.einsum('mbl,ma,ab->lab', conj_basis_mxs, L, F2.conj()) # ditto dVdp.shape = [dVdp.shape[0], nP] # jacobian with respect to "p" params, # which don't include normalization for TP-constraint @@ -363,28 +363,3 @@ def has_nonzero_hessian(self): bool """ return True - - def hessian_wrt_params(self, wrt_filter1=None, wrt_filter2=None): - """ - Construct the Hessian of this state vector with respect to its parameters. - - This function returns a tensor whose first axis corresponds to the - flattened operation matrix and whose 2nd and 3rd axes correspond to the - parameters that are differentiated with respect to. - - Parameters - ---------- - wrt_filter1 : list or numpy.ndarray - List of parameter indices to take 1st derivatives with respect to. - (None means to use all the this operation's parameters.) - - wrt_filter2 : list or numpy.ndarray - List of parameter indices to take 2nd derivatives with respect to. - (None means to use all the this operation's parameters.) - - Returns - ------- - numpy array - Hessian with shape (dimension, num_params1, num_params2) - """ - raise NotImplementedError("TODO: add hessian computation for CPTPState") From 9116bddefc510f6d90e24f526ddda0ce36618e08 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Mon, 2 Feb 2026 12:48:45 -0800 Subject: [PATCH 15/21] check-in new ProjectedState class --- pygsti/modelmembers/states/projectedstate.py | 345 +++++++++++++++++++ pygsti/tools/matrixtools.py | 50 ++- pygsti/tools/optools.py | 34 +- pygsti/tools/sdptools.py | 66 +++- 4 files changed, 471 insertions(+), 24 deletions(-) create mode 100644 pygsti/modelmembers/states/projectedstate.py diff --git a/pygsti/modelmembers/states/projectedstate.py b/pygsti/modelmembers/states/projectedstate.py new file mode 100644 index 000000000..ebc8bdaee --- /dev/null +++ b/pygsti/modelmembers/states/projectedstate.py @@ -0,0 +1,345 @@ +""" +A State that holds an internal TPState, and wraps that by projecting +as needed in from_vector(). We implement to_vector() by extracting +the projected density matrix and then returning +""" +from __future__ import annotations +from typing import TYPE_CHECKING, Literal, Any, Tuple + +if TYPE_CHECKING: + import cvxpy as cp + +import numpy as np +import jax +import jax.numpy as jnp +from pygsti.baseobjs.basis import ExplicitBasis, Basis + + +def strict_triu_op_data(N: int): + sT = N*(N-1)//2 + row_idx, col_idx = jnp.triu_indices(N, k=1) + P_rows = N * row_idx + col_idx + P_cols = np.arange(sT) + return (P_rows, P_cols), (N**2, sT) + + +def jax_vec_to_strict_upper_tri(vec): + ell = vec.shape[0] + n = round( ((8 * ell + 1) ** 0.5 + 1) / 2 ) + if not (n * (n - 1) // 2 == ell): + raise ValueError("The size of the vector must be a triangular number.") + coords, shape = strict_triu_op_data(n) + P_dense = jnp.zeros(shape) + P_dense = P_dense.at[coords].set(1.0) + out = (P_dense @ vec).reshape((n, n), order='F').T + return out + + +def tpvec_to_herm_parts(vec: cp.Expression | jnp.ndarray): + import cvxpy as cp + N2m1 = vec.shape[0] + N = round((N2m1 + 1)**0.5) + block0 = vec[ : N*(N-1)//2 ] + block1 = vec[ N*(N-1)//2 : N*(N-1) ] + block2 = vec[ N*(N-1) : ] + if isinstance(vec, np.ndarray): + vec = jnp.array(vec) + + if isinstance(vec, cp.Expression): + tri_real = cp.vec_to_upper_tri(block0, strict=True) + tri_imag = cp.vec_to_upper_tri(block1, strict=True) + diag = cp.diag(cp.concatenate([block2, 1 - cp.sum(block2)])) + elif isinstance(vec, jnp.ndarray): + tri_real = jax_vec_to_strict_upper_tri(block0) + tri_imag = jax_vec_to_strict_upper_tri(block1) + block3 = jnp.atleast_1d( 1 - jnp.sum(block2, axis=0) ) # type: ignore + diag = jnp.diag( + jnp.concatenate([block2, block3]) # type: ignore + ) + else: + raise ValueError() + + herm_real = tri_real + tri_real.T + diag + herm_imag = tri_imag - tri_imag.T + return herm_real, herm_imag + + +def herm_parts_to_tpvec(herm_parts: tuple[jnp.ndarray, jnp.ndarray] | tuple[Any, Any], full_superket=False): + herm_real, herm_imag = herm_parts + n = herm_imag.shape[0] + triui = np.triu_indices(n, k=1) + flat = (n*(n-1)//2,) + herm_imag.shape[2:] + tri_imag = herm_imag[triui].reshape(flat) + tri_real = herm_real[triui].reshape(flat) + num_diag = n if full_superket else n - 1 + diag = herm_real[np.diag_indices(num_diag)] + if isinstance(herm_real, jnp.ndarray): + return jnp.concatenate([tri_real, tri_imag, diag]) + else: + return np.concatenate([tri_real, tri_imag, diag]) + + +def herm_parts_basis(N: int) -> ExplicitBasis: + triui_rows, triu_cols = np.triu_indices(N, k=1) + sT = N*(N-1)//2 + mats = [] + labels = [] + for k in range(N**2): + mat = np.zeros((N, N), dtype=complex) + if k < sT: + # upper triangle, real, read in row-major + r, c = triui_rows[k], triu_cols[k] + label = '[ReH_%s]' % str({r,c}) + mat[r, c] = 1.0 + mat[c, r] = 1.0 + elif k < 2*sT: + # upper triangle, skew, read in row-major + r, c = triui_rows[k-sT], triu_cols[k-sT] + label = '[ImH_%s]' % str({r,c}) + mat[r, c] = 1.0j + mat[c, r] = -1.0j + else: + i = k - 2*sT + label = '[Diag_{%s}]' % i + mat[i, i] = 1.0 + mats.append(mat) + labels.append(label) + name = f'Hermitian matrix-unit basis of {N**2}-dimensional Hilbert-Schmidt space' + B = ExplicitBasis(mats, labels, name) + return B + + +def eigh_via_symmetricification(herm_parts: tuple[jnp.ndarray, jnp.ndarray] | tuple[Any, Any]) -> tuple[ jnp.ndarray | Any , tuple[Any, Any]]: + """ + Let A be real-symmetric and B be real skew-symmetric, both of order n. + If c is a scalar and x and y are real n-vectors satisfying + + c [x] = [ A B ][x] + [y] = [ B' A ][y]. + + then x - 1j * y is an eigenvector of A + 1j * B with eigenvalue c. + + This function returns a complete set of eigenvalues and eigenvectors + for the matrix herm_parts[0] + 1j * herm_parts[1], where the + eigenvectors are returned in separate matrices containing their real + and imaginary parts. + + The implementation assumes that `eigh` returns eigenpairs ordered by + eigenvalue (either increasing or decreasing is fine). + """ + herm_real, herm_imag = herm_parts + n = herm_imag.shape[0] + if isinstance(herm_real, jnp.ndarray): + sym = jnp.block([[herm_real, herm_imag],[herm_imag.T, herm_real]]) + sevals, sevecs = jnp.linalg.eigh(sym) + else: + sym = np.block([[herm_real, herm_imag],[herm_imag.T, herm_real]]) + sevals, sevecs = np.linalg.eigh(sym) + evals = sevals[::2] + evecs_real = sevecs[ :n, ::2 ] + evecs_imag = -sevecs[ n:, ::2 ] + return evals, (evecs_real, evecs_imag) + + +def recollect_expanded_eigendecomposition(evals, evecs_real, evecs_imag): + """ + evals, evecs_real, and evecs_imag are all real-valued. + + Define D = diag(evals) and U = evecs_real + 1j * evecs_imag. + + This function returns matrices herm_real and herm_imag where + + herm_real + 1j * herm_imag = U @ D @ U.T.conj(). + """ + temp_real = evals[:, None] * evecs_real.T + temp_imag = evals[:, None] * evecs_imag.T + herm_real = evecs_real @ temp_real + evecs_imag @ temp_imag + herm_imag = evecs_imag @ temp_real - evecs_real @ temp_imag + return herm_real, herm_imag + + +def density_projection_model(N: int, param_name='parameters', var_name='projection'): + import cvxpy as cp + param = cp.Parameter(shape=(N**2 - 1,), name=param_name) + herm_real, herm_imag = tpvec_to_herm_parts(param) + herm_expr = herm_real + 1j * herm_imag + + X = cp.Variable(shape=(N, N), hermitian=True, name=var_name) + F = cp.sum_squares(herm_expr - X) + C = [ 0 << X, cp.trace(cp.real(X)) == 1 ] + problem = cp.Problem( cp.Minimize(F), C ) + + """ + To map V := X.value back to `param` space, + set triui = np.triu_indices_from(V, k=1), + then concatenate ... + V.real[triui], + V.imag[triui], and + np.diag(V.real)[:-1]. + """ + return problem + + +def project_onto_simplex(vec: jnp.ndarray) -> jnp.ndarray: + # Sort in descending order + u = jnp.sort(vec)[::-1] + css = jnp.cumsum(u, axis=0) + # Find the largest index rho such that + # u[rho] > (css[rho] - 1) / (rho + 1) + inds = jnp.arange(vec.shape[0]) + 1 + cond = u - (css - 1) / inds > 0 + if not jnp.any(cond): + # if no positive entries, project uniformly on first coordinate + theta = (css[0] - 1) / 1 + else: + rho = jnp.where(cond)[0][-1] + theta = (css[rho] - 1) / (rho + 1) + proj = jnp.maximum(vec - theta, 0) + return proj + + +def project_onto_densities(herm_mx: jnp.ndarray) -> jnp.ndarray: + v, u = jnp.linalg.eigh(herm_mx) + udag = u.T.conj() + proj_v = project_onto_simplex(v) + density_mx = u @ (proj_v[:, None] * udag) + return density_mx + + +def project_onto_densities_realification(herm_parts: jnp.ndarray) -> jnp.ndarray: + """ + Example usage: + + jac = jax.jacobian(project_onto_densities_realification) + + herm_mx = np.random.randn(3,3) + 1j * np.random.randn(3,3) + herm_mx += temp.T.conj() + + herm_mx_reif = jnp.stack((herm_mx.real, temp.imag)) + # ^ shape (2, 3, 3) + out = jac(herm_mx_reif) + # ^ shape (2, 3, 3, 2, 3, 3) + + """ + evals, evecs = eigh_via_symmetricification(herm_parts) # type: ignore + proj_evals = project_onto_simplex(evals) + density_mx = recollect_expanded_eigendecomposition(proj_evals, *evecs) + density_mx = jnp.stack(density_mx) + return density_mx + + +def project_onto_densities_tpvec( tpvec, full_superket=False ): + herm_parts = tpvec_to_herm_parts( tpvec ) + herm_mx = jnp.stack(herm_parts) # type: ignore + density_mx = project_onto_densities_realification(herm_mx) + proj_parts = (density_mx[0], density_mx[1]) + tpvec_proj = herm_parts_to_tpvec(proj_parts, full_superket) + return tpvec_proj + + + +class ProjectedState_CVXPY: + + def __init__(self, N: int) -> None: + problem = density_projection_model(N) + self._cvxpy_model = problem + self._cvxpy_param = problem.param_dict['parameters'] + return + +from pygsti.modelmembers import ModelMember +from pygsti.modelmembers.states import State, TPState +from pygsti.tools.basistools import vec_to_stdmx, change_basis + + +# TODO: need a routine for changing basis of to_vector()/from_vector() with TPState +# versus ProjectedState. + + + + + +class ProjectedState(TPState): + + def __init__(self, vec: np.ndarray, basis: Basis, evotype="default", state_space=None): + super().__init__(vec, basis, evotype, state_space) + self._working_basis = herm_parts_basis(self.state_space.udim) + + self._n2w_np = self._basis.create_transform_matrix(self._working_basis).real + self._w2n_np = self._working_basis.create_transform_matrix(self._basis).real + self._n2w = jnp.array( self._n2w_np ) + self._w2n = jnp.array( self._w2n_np ) + + self._natural_firstEl = np.array([self.basis.elsize ** -0.25]) + return + + def to_vector(self) -> np.ndarray: + return super().to_vector() + + def _natural_vec_to_working_vec(self, natvec): + if isinstance(natvec, np.ndarray): + workvec = self._n2w_np[:-1, :] @ np.concatenate([self._natural_firstEl, natvec]) + else: + workvec = self._n2w[:-1, :] @ jnp.concatenate([self._natural_firstEl, natvec]) + return workvec + + def _working_vec_to_natural_vec(self, workvec): + if isinstance(workvec, np.ndarray): + tail = np.atleast_1d(1 - np.sum(workvec[-(self.state_space.udim - 1):])) + natvec = self._w2n_np[1:, :] @ np.concatenate([workvec, tail]) + else: + tail = jnp.atleast_1d(1 - jnp.sum(workvec[-(self.state_space.udim - 1):])) + natvec = self._w2n[1:, :] @ jnp.concatenate([workvec, tail]) + return natvec + + def from_vector(self, v, close=False, dirty_value=True) -> None: + v = jnp.array(v) + tpvec_before = self._natural_vec_to_working_vec(v) + tpvec_after = project_onto_densities_tpvec( tpvec_before ) + v_projected = self._working_vec_to_natural_vec( tpvec_after ) + v_projected = np.asarray(v_projected).astype(np.double) + super().from_vector(v_projected, False, True) + return + + def set_dense(self, vec): + super().set_dense(vec) + self.from_vector(self.to_vector(), False, True) + return + + def deriv_wrt_params(self, wrt_filter=None): + assert wrt_filter is None + + def projected_tpvec_onto_densities_as_superket( tpvec ): + todense_working = project_onto_densities_tpvec( tpvec, full_superket=True ) + todense_natural = self._w2n @ todense_working + return todense_natural + + jac_fun = jax.jacobian(projected_tpvec_onto_densities_as_superket) + natural_sket = jnp.array(self.to_dense('HilbertSchmidt')) + current_tpvec = self._n2w[:-1, :] @ natural_sket + jac_val = jac_fun( current_tpvec ) + jac_val = jac_val @ self._n2w[:-1, :] + # ^ Jacobian mapping natural_sket_in to natural_sket_out. + jac_val = jac_val[:, 1:] + # ^ Jacobian mapping the last `from_vector` arg `v` to self.to_dense(). + + jac_val = np.asarray( jac_val, dtype=np.double ) + return jac_val + + def stateless_data(self) -> Tuple[int]: + raise NotImplementedError() + + @staticmethod + def torch_base(sd, t_param): + raise NotImplementedError() + + def has_nonzero_hessian(self): + return False + + +if __name__ == '__main__': + temp = np.random.randn(3,3) + 1j * np.random.randn(3,3) + temp += temp.T.conj() + tpvec = np.random.rand(3) + jac2 = jax.jacobian(project_onto_densities_tpvec) + x,y = tpvec_to_herm_parts(jnp.array(tpvec)) + print() diff --git a/pygsti/tools/matrixtools.py b/pygsti/tools/matrixtools.py index d43ac70f8..85a20e347 100644 --- a/pygsti/tools/matrixtools.py +++ b/pygsti/tools/matrixtools.py @@ -14,7 +14,7 @@ import itertools as _itertools import warnings as _warnings -from typing import Protocol, Any, runtime_checkable, TypeVar, Optional, Union +from typing import Protocol, Any, runtime_checkable, TypeVar, Optional, Union, Literal import numpy as _np import scipy.linalg as _spl @@ -150,6 +150,54 @@ def is_valid_density_mx(mx, tol=1e-9): return abs(_np.trace(mx) - 1.0) < tol and is_pos_def(mx, tol) +def project_onto_simplex(vec: _np.ndarray, ord:Literal[1,2,'inf']) -> _np.ndarray: + """ + Return the solution, `proj`, to the optimization problem + + min{ ||vec - proj||_ord : 0 <= proj, sum(proj) == 1 }, + + where ||*||_ord refers to the standard vector `ord`-norm. + """ + vec = vec.ravel() + if ord == 2: + # Sort in descending order + u = _np.sort(vec)[::-1] + css = _np.cumsum(u) + # Find the largest index rho such that + # u[rho] > (css[rho] - 1) / (rho + 1) + inds = _np.arange(vec.size) + 1 + cond = u - (css - 1) / inds > 0 + if not _np.any(cond): + # fallback: if no positive entries, project uniformly on first coordinate + theta = (css[0] - 1) / 1 + else: + rho = _np.where(cond)[0][-1] + theta = (css[rho] - 1) / (rho + 1) + # The projection is max(flat - theta, 0) + proj = _np.maximum(vec - theta, 0) + else: + raise NotImplementedError() + return proj + + +def project_onto_densities(herm_mx: _np.ndarray, ord:Literal['fro','nuc']='fro') -> _np.ndarray: + if ord != 'fro': + raise NotImplementedError() + assert_hermitian(herm_mx, 1e-8) + u, v, udag = eigendecomposition(herm_mx, assume_hermitian=True) + proj_v = project_onto_simplex(v, 2) + density_mx = u @ (proj_v[:, None] * udag) + return density_mx + + +def clip_eigenvalues(herm_mx: _np.ndarray, lower: float, upper:float) -> _np.ndarray: + assert_hermitian(herm_mx, 1e-8) + u, v, udag = eigendecomposition(herm_mx, assume_hermitian=True) + proj_v = _np.clip(v, lower, upper) + clipped_mx = u @ (proj_v[:, None] * udag) + return clipped_mx + + def nullspace(m, tol=1e-7): """ Compute the nullspace of a matrix. diff --git a/pygsti/tools/optools.py b/pygsti/tools/optools.py index 84a8dcd2f..6adfd5edd 100644 --- a/pygsti/tools/optools.py +++ b/pygsti/tools/optools.py @@ -373,7 +373,6 @@ def diamonddist(a, b, mx_basis='pp', return_x=False): `dm = trace( |(J(a)-J(b)).T * W| )`. """ assert _sdps.CVXPY_ENABLED - import cvxpy as _cp assert a.shape == b.shape mx_basis = _bt.create_basis_for_matrix(a, mx_basis) @@ -382,31 +381,22 @@ def diamonddist(a, b, mx_basis='pp', return_x=False): # It will convert c a "single-block" basis representation # when mx_basis has multiple blocks. J = _jam.fast_jamiolkowski_iso_std(c, mx_basis, normalized=False) - prob, vars = _sdps.diamond_norm_model_jamiolkowski(J) + problem, vars = _sdps.diamond_norm_model_jamiolkowski(J) + objval, varvals = _sdps.solve_sdp(problem) - objective_val = -2 - varvals = [_np.zeros_like(J), None, None] - sdp_solvers = ['MOSEK', 'CLARABEL', 'CVXOPT'] - for i, solver in enumerate(sdp_solvers): - try: - prob.solve(solver=solver) - objective_val = prob.value - varvals = [v.value for v in vars] - break - except (AssertionError, _cp.SolverError) as e: - if solver != 'MOSEK': - msg = f"Received error {e} when trying to use solver={solver}." - if i + 1 == len(sdp_solvers): - failure_msg = "Out of solvers. Returning -2 for diamonddist." - else: - failure_msg = f"Trying {sdp_solvers[i+1]} next." - msg += f'\n{failure_msg}' - _warnings.warn(msg) + if _np.isnan(objval): + objval = -2 + varvals = [_np.zeros_like(J), None, None] + else: + # the returned `varvals` is a dict keyed by + # variable name. We care about the variables + # as ordered by `vars`. + varvals = [v.value for v in vars] if return_x: - return objective_val, varvals + return objval, varvals else: - return objective_val + return objval 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 index 071195a20..02f2e6248 100644 --- a/pygsti/tools/sdptools.py +++ b/pygsti/tools/sdptools.py @@ -13,13 +13,15 @@ from __future__ import annotations import numpy as np +import warnings -from typing import Union, List, Tuple, TYPE_CHECKING +from typing import Union, List, Tuple, Sequence, TYPE_CHECKING if TYPE_CHECKING: import cvxpy as cp ExpressionLike = Union[cp.Expression, np.ndarray] from pygsti.baseobjs import Basis +from pygsti.tools.matrixtools import assert_hermitian from pygsti.tools.basistools import stdmx_to_vec from pygsti.tools.jamiolkowski import jamiolkowski_iso @@ -35,6 +37,68 @@ CVXPY_ENABLED = False +SDP_SOLVER_PRIORITY = ['MOSEK', 'CLARABEL', 'CVXOPT'] + + +def solve_sdp(prob: cp.Problem) -> tuple[np.floating, dict[str, np.ndarray]]: + + objective_val = np.array(np.NaN).item() + varvals = dict() + for i, solver in enumerate(SDP_SOLVER_PRIORITY): + try: + prob.solve(solver=solver) + objective_val = prob.value + varvals.update({k: v.value for (k, v) in prob.var_dict}) + break + except (AssertionError, cp.SolverError) as e: + if solver != 'MOSEK': + msg = f"Received error {e} when trying to use solver={solver}." + if i + 1 == len(SDP_SOLVER_PRIORITY): + failure_msg = "Out of solvers. Returning NaN." + else: + failure_msg = f"Trying {SDP_SOLVER_PRIORITY[i+1]} next." + msg += f'\n{failure_msg}' + warnings.warn(msg) + + return objective_val, varvals + + +def povm_projection_model(effects: Sequence[np.ndarray], independent_complement: bool=False) -> cp.Problem: + N = round(effects[0].size ** 0.5) + cumsum_E = np.zeros((N, N), dtype=complex) + for E in effects: + assert_hermitian(E, tol=1e-8) + cumsum_E += E + if np.real(np.trace(cumsum_E)) < (N - 0.01): + # we're going to assume the complement POVM effect was omitted. + complement_E = np.eye(N) - np.sum(effects, axis=0) + effects = [E for E in effects] # shallow-copy + effects.append( complement_E ) + + proj_effects, constraints = povm_effect_variables(N, len(effects), independent_complement) + objective_expr = cp.Constant(0.0) + for E_proj, E in zip( proj_effects, effects ): + objective_expr = objective_expr + cp.sum_squares(E_proj - E) + objective = cp.Minimize(objective_expr) + + prob = cp.Problem(objective, constraints) + return prob + + +def povm_effect_variables(N: int, num_effects: int, independent_complement: bool) -> tuple[list[cp.Variable], list[cp.Constraint]]: + num_vars = num_effects if independent_complement else num_effects - 1 + vars = [cp.Variable(shape=(N, N), hermitian=True) for _ in range(num_vars)] + cons = [v >> 0 for v in vars] + if num_vars == num_effects: + expr = cp.sum(vars) + cons.append( expr == np.eye(N) ) + else: + expr = np.eye(N) - cp.sum(vars) + cons.append( expr >> 0 ) + vars.append( expr ) + return vars, cons + + def diamond_norm_model_jamiolkowski(J: ExpressionLike) -> tuple[cp.Problem, List[cp.Variable]]: # return a model for computing the diamond norm. # From 1b195c869de47ad6f476afd61fcfa4881f1263a5 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Fri, 13 Mar 2026 13:44:52 -0700 Subject: [PATCH 16/21] QoL updates to sdptools.py --- pygsti/tools/sdptools.py | 25 ++++++++++++------------- 1 file changed, 12 insertions(+), 13 deletions(-) diff --git a/pygsti/tools/sdptools.py b/pygsti/tools/sdptools.py index 0417d28eb..704aa463b 100644 --- a/pygsti/tools/sdptools.py +++ b/pygsti/tools/sdptools.py @@ -15,7 +15,7 @@ import numpy as np import warnings -from typing import Union, List, Tuple, Sequence, TYPE_CHECKING +from typing import Any, Union, List, Tuple, Sequence, TYPE_CHECKING if TYPE_CHECKING: import cvxpy as cp ExpressionLike = Union[cp.Expression, np.ndarray] @@ -40,18 +40,17 @@ SDP_SOLVER_PRIORITY = ['MOSEK', 'CLARABEL', 'CVXOPT'] -def solve_sdp(prob: cp.Problem) -> tuple[np.floating, dict[str, np.ndarray]]: +def solve_sdp(prob: cp.Problem, **kwargs) -> tuple[np.floating, dict[str, np.ndarray]]: - objective_val = np.array(np.NaN).item() - varvals = dict() + objective_val : np.floating = np.array(np.NaN).item() + varvals : dict[str, np.ndarray] = dict() for i, solver in enumerate(SDP_SOLVER_PRIORITY): try: with warnings.catch_warnings(): warnings.filterwarnings('ignore','.*Solution may be inaccurate.*', UserWarning) - prob.solve(solver=solver) - prob.solve(solver=solver) - objective_val = prob.value - varvals.update({k: v.value for (k, v) in prob.var_dict}) + prob.solve(solver=solver, **kwargs) + objective_val = prob.value # type: ignore + varvals.update({k: v.value for (k, v) in prob.var_dict.items()}) # type: ignore break except (AssertionError, cp.SolverError) as e: if solver != 'MOSEK': @@ -88,15 +87,15 @@ def povm_projection_model(effects: Sequence[np.ndarray], independent_complement: return prob -def povm_effect_variables(N: int, num_effects: int, independent_complement: bool) -> tuple[list[cp.Variable], list[cp.Constraint]]: +def povm_effect_variables(N: int, num_effects: int, independent_complement: bool) -> tuple[list[cp.Expression], list[cp.Constraint]]: num_vars = num_effects if independent_complement else num_effects - 1 - vars = [cp.Variable(shape=(N, N), hermitian=True) for _ in range(num_vars)] - cons = [v >> 0 for v in vars] + vars : list[cp.Expression] = [cp.Variable(shape=(N, N), hermitian=True) for _ in range(num_vars)] # type: ignore + cons : list[cp.Constraint] = [v >> 0 for v in vars] if num_vars == num_effects: - expr = cp.sum(vars) + expr : cp.Expression = cp.sum(vars) # type: ignore cons.append( expr == np.eye(N) ) else: - expr = np.eye(N) - cp.sum(vars) + expr : cp.Expression = np.eye(N) - cp.sum(vars) # type: ignore cons.append( expr >> 0 ) vars.append( expr ) return vars, cons From af571515837ce4c371142d6e99793f37392a0da4 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Fri, 13 Mar 2026 13:48:28 -0700 Subject: [PATCH 17/21] add CVXPY-backed ProjectedState prototype --- pygsti/modelmembers/states/projectedstate.py | 152 +++++++++++++++++-- 1 file changed, 139 insertions(+), 13 deletions(-) diff --git a/pygsti/modelmembers/states/projectedstate.py b/pygsti/modelmembers/states/projectedstate.py index ebc8bdaee..f0d9142ac 100644 --- a/pygsti/modelmembers/states/projectedstate.py +++ b/pygsti/modelmembers/states/projectedstate.py @@ -4,16 +4,23 @@ the projected density matrix and then returning """ from __future__ import annotations -from typing import TYPE_CHECKING, Literal, Any, Tuple +from typing import TYPE_CHECKING, Literal, Any, Tuple, Optional, Callable if TYPE_CHECKING: import cvxpy as cp +import cvxpy as cp import numpy as np import jax +jax.config.update("jax_enable_x64", True) import jax.numpy as jnp from pygsti.baseobjs.basis import ExplicitBasis, Basis +from pygsti.modelmembers import ModelMember +from pygsti.modelmembers.states import State, TPState +from pygsti.tools.basistools import vec_to_stdmx, change_basis, stdmx_to_vec +from pygsti.tools import sdptools as _sdps + def strict_triu_op_data(N: int): sT = N*(N-1)//2 @@ -48,7 +55,8 @@ def tpvec_to_herm_parts(vec: cp.Expression | jnp.ndarray): if isinstance(vec, cp.Expression): tri_real = cp.vec_to_upper_tri(block0, strict=True) tri_imag = cp.vec_to_upper_tri(block1, strict=True) - diag = cp.diag(cp.concatenate([block2, 1 - cp.sum(block2)])) + tail = 1 - cp.sum(block2) + diag = cp.diag(cp.concatenate([block2, tail.reshape((1,))])) elif isinstance(vec, jnp.ndarray): tri_real = jax_vec_to_strict_upper_tri(block0) tri_imag = jax_vec_to_strict_upper_tri(block1) @@ -158,13 +166,53 @@ def recollect_expanded_eigendecomposition(evals, evecs_real, evecs_imag): return herm_real, herm_imag -def density_projection_model(N: int, param_name='parameters', var_name='projection'): +def demo_derivatives_cvxpy(): + import cvxpy as cp + + __var = cp.Variable(shape=(3,), pos=True) + x = __var[0] + y = __var[1] + z = __var[2] + + __param_ab = cp.Parameter(shape=(2,), pos=True) + a = __param_ab[0] + b = __param_ab[1] + c = 0.5 + + objective_fn = 1/(x*y*z) + objective = cp.Minimize(objective_fn) + constraints = [a*(x*y + x*z + y*z) <= b, x >= y**c] + problem = cp.Problem(objective, constraints) + + problem.is_dgp(dpp=True) + + __param_ab.value = np.array([2.0, 1.0]) + problem.solve(gp=True, requires_grad=True, solver='CLARABEL') + + d__var = [] + EPS = 1e-4 + for i in range(__param_ab.size): + __param_ab.delta = np.zeros_like(__param_ab.value) + __param_ab.delta[i] = EPS + problem.derivative() + d__var.append(__var.delta.copy() / EPS ) + + return d__var + + +def density_projection_model(N: int, firstEl: np.ndarray, basis, param_name='parameters', var_name='projection', param_transform: Optional[Callable]=None): import cvxpy as cp param = cp.Parameter(shape=(N**2 - 1,), name=param_name) + if param_transform is not None: + param = param_transform(param) herm_real, herm_imag = tpvec_to_herm_parts(param) herm_expr = herm_real + 1j * herm_imag - X = cp.Variable(shape=(N, N), hermitian=True, name=var_name) + free_params = cp.Variable(shape=(N**2 - 1,), name=var_name) + assert firstEl.size == 1 + superket = cp.concatenate([firstEl.ravel(), free_params]) + basisMx = np.column_stack([B.ravel() for B in basis.elements]) + X = cp.hermitian_wrap((basisMx @ superket).reshape((N, N))) F = cp.sum_squares(herm_expr - X) C = [ 0 << X, cp.trace(cp.real(X)) == 1 ] problem = cp.Problem( cp.Minimize(F), C ) @@ -238,23 +286,99 @@ def project_onto_densities_tpvec( tpvec, full_superket=False ): -class ProjectedState_CVXPY: +class ProjectedState_CVXPY(TPState): + + def __init__(self, vec: np.ndarray, basis: Basis, evotype="default", state_space=None): + super().__init__(vec, basis, evotype, state_space) + self._working_basis = herm_parts_basis(self.state_space.udim) - def __init__(self, N: int) -> None: - problem = density_projection_model(N) - self._cvxpy_model = problem - self._cvxpy_param = problem.param_dict['parameters'] + self._n2w_np = self._basis.create_transform_matrix(self._working_basis).real + self._w2n_np = self._working_basis.create_transform_matrix(self._basis).real + + self._natural_firstEl = np.array([self._basis.elsize ** -0.25]) + self._cvxpy_model = density_projection_model( + self.state_space.udim, self._natural_firstEl, + self._basis, param_transform=self._natural_vec_to_working_vec + ) + # ^ the param_transform means that the + self._cvxpy_param = self._cvxpy_model.param_dict['parameters'] return -from pygsti.modelmembers import ModelMember -from pygsti.modelmembers.states import State, TPState -from pygsti.tools.basistools import vec_to_stdmx, change_basis + @staticmethod + def _project_onto_densities_natural(model, param, natvec): + param.value = natvec + _, vars = _sdps.solve_sdp(model) + superket = vars['projection'] + return superket + + def to_vector(self) -> np.ndarray: + return super().to_vector() + + def _natural_vec_to_working_vec(self, natvec): + if isinstance(natvec, cp.Expression): + sket = cp.concatenate([self._natural_firstEl, natvec]) + else: + sket = np.concatenate([self._natural_firstEl, natvec]) + workvec = self._n2w_np[:-1, :] @ sket + return workvec + + def from_vector(self, v, close=False, dirty_value=True) -> None: + v_projected = ProjectedState_CVXPY._project_onto_densities_natural( self._cvxpy_model, self._cvxpy_param, v ) + super().from_vector(v_projected, False, True) + return + + def set_dense(self, vec): + super().set_dense(vec) + self.from_vector(self.to_vector(), False, True) + return + + def deriv_wrt_params(self, wrt_filter=None): + # return the Jacobian of the map from [the vector `v` most recently + # passed to `self.from_vector(v)`] to [self.to_dense()]. + assert wrt_filter is None + ALL_SOLVERS = _sdps.SDP_SOLVER_PRIORITY + _sdps.SDP_SOLVER_PRIORITY = ['CLARABEL'] + _sdps.solve_sdp(self._cvxpy_model, requires_grad=True) + var = self._cvxpy_model.var_dict['projection'] + param = self._cvxpy_param + jac_cols = [] + EPS = 1e-4 + param.value = self.to_vector() + for i in range(param.size): + param.delta = np.zeros_like(param.value) + param.delta[i] = EPS + self._cvxpy_model.derivative() + jac_cols.append(var.delta.copy() / EPS) + jac_val = np.column_stack(jac_cols) + jac_val = np.row_stack([np.zeros((1, self.state_space.dim-1)), jac_val]) + _sdps.SDP_SOLVER_PRIORITY = ALL_SOLVERS + return jac_val + + def stateless_data(self) -> Tuple[int]: + raise NotImplementedError() + + @staticmethod + def torch_base(sd, t_param): + raise NotImplementedError() + + def has_nonzero_hessian(self): + return False + # TODO: need a routine for changing basis of to_vector()/from_vector() with TPState # versus ProjectedState. +""" +TODO: modify to_vector so that it returns the raw un-projected v from the last from_vector. +Only have the projection happen in to_dense() and deriv_wrt_params(). This viewpoint +means the underlying vector in a ProjectedState is viewed through an equivalence class +based on how the resulting Hermitian matrix projects onto the set of density matrices. + +^ Would be desirable if we wanted to satisfy + assert( modelmember.to_vector(modelmember.from_vector(v)) == v) +""" @@ -269,7 +393,7 @@ def __init__(self, vec: np.ndarray, basis: Basis, evotype="default", state_space self._n2w = jnp.array( self._n2w_np ) self._w2n = jnp.array( self._w2n_np ) - self._natural_firstEl = np.array([self.basis.elsize ** -0.25]) + self._natural_firstEl = np.array([self._basis.elsize ** -0.25]) return def to_vector(self) -> np.ndarray: @@ -343,3 +467,5 @@ def has_nonzero_hessian(self): jac2 = jax.jacobian(project_onto_densities_tpvec) x,y = tpvec_to_herm_parts(jnp.array(tpvec)) print() + demo_derivatives_cvxpy() + print() \ No newline at end of file From d2cf2adb920cb4389cdd078195dcabee1bf6b92c Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Fri, 13 Mar 2026 14:46:38 -0700 Subject: [PATCH 18/21] Revert "add CVXPY-backed ProjectedState prototype" This reverts commit af571515837ce4c371142d6e99793f37392a0da4. --- pygsti/modelmembers/states/projectedstate.py | 152 ++----------------- 1 file changed, 13 insertions(+), 139 deletions(-) diff --git a/pygsti/modelmembers/states/projectedstate.py b/pygsti/modelmembers/states/projectedstate.py index f0d9142ac..ebc8bdaee 100644 --- a/pygsti/modelmembers/states/projectedstate.py +++ b/pygsti/modelmembers/states/projectedstate.py @@ -4,23 +4,16 @@ the projected density matrix and then returning """ from __future__ import annotations -from typing import TYPE_CHECKING, Literal, Any, Tuple, Optional, Callable +from typing import TYPE_CHECKING, Literal, Any, Tuple if TYPE_CHECKING: import cvxpy as cp -import cvxpy as cp import numpy as np import jax -jax.config.update("jax_enable_x64", True) import jax.numpy as jnp from pygsti.baseobjs.basis import ExplicitBasis, Basis -from pygsti.modelmembers import ModelMember -from pygsti.modelmembers.states import State, TPState -from pygsti.tools.basistools import vec_to_stdmx, change_basis, stdmx_to_vec -from pygsti.tools import sdptools as _sdps - def strict_triu_op_data(N: int): sT = N*(N-1)//2 @@ -55,8 +48,7 @@ def tpvec_to_herm_parts(vec: cp.Expression | jnp.ndarray): if isinstance(vec, cp.Expression): tri_real = cp.vec_to_upper_tri(block0, strict=True) tri_imag = cp.vec_to_upper_tri(block1, strict=True) - tail = 1 - cp.sum(block2) - diag = cp.diag(cp.concatenate([block2, tail.reshape((1,))])) + diag = cp.diag(cp.concatenate([block2, 1 - cp.sum(block2)])) elif isinstance(vec, jnp.ndarray): tri_real = jax_vec_to_strict_upper_tri(block0) tri_imag = jax_vec_to_strict_upper_tri(block1) @@ -166,53 +158,13 @@ def recollect_expanded_eigendecomposition(evals, evecs_real, evecs_imag): return herm_real, herm_imag -def demo_derivatives_cvxpy(): - import cvxpy as cp - - __var = cp.Variable(shape=(3,), pos=True) - x = __var[0] - y = __var[1] - z = __var[2] - - __param_ab = cp.Parameter(shape=(2,), pos=True) - a = __param_ab[0] - b = __param_ab[1] - c = 0.5 - - objective_fn = 1/(x*y*z) - objective = cp.Minimize(objective_fn) - constraints = [a*(x*y + x*z + y*z) <= b, x >= y**c] - problem = cp.Problem(objective, constraints) - - problem.is_dgp(dpp=True) - - __param_ab.value = np.array([2.0, 1.0]) - problem.solve(gp=True, requires_grad=True, solver='CLARABEL') - - d__var = [] - EPS = 1e-4 - for i in range(__param_ab.size): - __param_ab.delta = np.zeros_like(__param_ab.value) - __param_ab.delta[i] = EPS - problem.derivative() - d__var.append(__var.delta.copy() / EPS ) - - return d__var - - -def density_projection_model(N: int, firstEl: np.ndarray, basis, param_name='parameters', var_name='projection', param_transform: Optional[Callable]=None): +def density_projection_model(N: int, param_name='parameters', var_name='projection'): import cvxpy as cp param = cp.Parameter(shape=(N**2 - 1,), name=param_name) - if param_transform is not None: - param = param_transform(param) herm_real, herm_imag = tpvec_to_herm_parts(param) herm_expr = herm_real + 1j * herm_imag - free_params = cp.Variable(shape=(N**2 - 1,), name=var_name) - assert firstEl.size == 1 - superket = cp.concatenate([firstEl.ravel(), free_params]) - basisMx = np.column_stack([B.ravel() for B in basis.elements]) - X = cp.hermitian_wrap((basisMx @ superket).reshape((N, N))) + X = cp.Variable(shape=(N, N), hermitian=True, name=var_name) F = cp.sum_squares(herm_expr - X) C = [ 0 << X, cp.trace(cp.real(X)) == 1 ] problem = cp.Problem( cp.Minimize(F), C ) @@ -286,99 +238,23 @@ def project_onto_densities_tpvec( tpvec, full_superket=False ): -class ProjectedState_CVXPY(TPState): - - def __init__(self, vec: np.ndarray, basis: Basis, evotype="default", state_space=None): - super().__init__(vec, basis, evotype, state_space) - self._working_basis = herm_parts_basis(self.state_space.udim) - - self._n2w_np = self._basis.create_transform_matrix(self._working_basis).real - self._w2n_np = self._working_basis.create_transform_matrix(self._basis).real - - self._natural_firstEl = np.array([self._basis.elsize ** -0.25]) - self._cvxpy_model = density_projection_model( - self.state_space.udim, self._natural_firstEl, - self._basis, param_transform=self._natural_vec_to_working_vec - ) - # ^ the param_transform means that the - self._cvxpy_param = self._cvxpy_model.param_dict['parameters'] - return - - @staticmethod - def _project_onto_densities_natural(model, param, natvec): - param.value = natvec - _, vars = _sdps.solve_sdp(model) - superket = vars['projection'] - return superket +class ProjectedState_CVXPY: - def to_vector(self) -> np.ndarray: - return super().to_vector() - - def _natural_vec_to_working_vec(self, natvec): - if isinstance(natvec, cp.Expression): - sket = cp.concatenate([self._natural_firstEl, natvec]) - else: - sket = np.concatenate([self._natural_firstEl, natvec]) - workvec = self._n2w_np[:-1, :] @ sket - return workvec - - def from_vector(self, v, close=False, dirty_value=True) -> None: - v_projected = ProjectedState_CVXPY._project_onto_densities_natural( self._cvxpy_model, self._cvxpy_param, v ) - super().from_vector(v_projected, False, True) + def __init__(self, N: int) -> None: + problem = density_projection_model(N) + self._cvxpy_model = problem + self._cvxpy_param = problem.param_dict['parameters'] return - - def set_dense(self, vec): - super().set_dense(vec) - self.from_vector(self.to_vector(), False, True) - return - - def deriv_wrt_params(self, wrt_filter=None): - # return the Jacobian of the map from [the vector `v` most recently - # passed to `self.from_vector(v)`] to [self.to_dense()]. - assert wrt_filter is None - ALL_SOLVERS = _sdps.SDP_SOLVER_PRIORITY - _sdps.SDP_SOLVER_PRIORITY = ['CLARABEL'] - _sdps.solve_sdp(self._cvxpy_model, requires_grad=True) - var = self._cvxpy_model.var_dict['projection'] - param = self._cvxpy_param - jac_cols = [] - EPS = 1e-4 - param.value = self.to_vector() - for i in range(param.size): - param.delta = np.zeros_like(param.value) - param.delta[i] = EPS - self._cvxpy_model.derivative() - jac_cols.append(var.delta.copy() / EPS) - jac_val = np.column_stack(jac_cols) - jac_val = np.row_stack([np.zeros((1, self.state_space.dim-1)), jac_val]) - _sdps.SDP_SOLVER_PRIORITY = ALL_SOLVERS - return jac_val - - def stateless_data(self) -> Tuple[int]: - raise NotImplementedError() - - @staticmethod - def torch_base(sd, t_param): - raise NotImplementedError() - - def has_nonzero_hessian(self): - return False +from pygsti.modelmembers import ModelMember +from pygsti.modelmembers.states import State, TPState +from pygsti.tools.basistools import vec_to_stdmx, change_basis # TODO: need a routine for changing basis of to_vector()/from_vector() with TPState # versus ProjectedState. -""" -TODO: modify to_vector so that it returns the raw un-projected v from the last from_vector. -Only have the projection happen in to_dense() and deriv_wrt_params(). This viewpoint -means the underlying vector in a ProjectedState is viewed through an equivalence class -based on how the resulting Hermitian matrix projects onto the set of density matrices. - -^ Would be desirable if we wanted to satisfy - assert( modelmember.to_vector(modelmember.from_vector(v)) == v) -""" @@ -393,7 +269,7 @@ def __init__(self, vec: np.ndarray, basis: Basis, evotype="default", state_space self._n2w = jnp.array( self._n2w_np ) self._w2n = jnp.array( self._w2n_np ) - self._natural_firstEl = np.array([self._basis.elsize ** -0.25]) + self._natural_firstEl = np.array([self.basis.elsize ** -0.25]) return def to_vector(self) -> np.ndarray: @@ -467,5 +343,3 @@ def has_nonzero_hessian(self): jac2 = jax.jacobian(project_onto_densities_tpvec) x,y = tpvec_to_herm_parts(jnp.array(tpvec)) print() - demo_derivatives_cvxpy() - print() \ No newline at end of file From a8d9c733ca18669e589265dd30008da6b8951096 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Fri, 13 Mar 2026 14:47:38 -0700 Subject: [PATCH 19/21] remove hacky proof of concept implementation --- pygsti/modelmembers/states/projectedstate.py | 345 ------------------- 1 file changed, 345 deletions(-) delete mode 100644 pygsti/modelmembers/states/projectedstate.py diff --git a/pygsti/modelmembers/states/projectedstate.py b/pygsti/modelmembers/states/projectedstate.py deleted file mode 100644 index ebc8bdaee..000000000 --- a/pygsti/modelmembers/states/projectedstate.py +++ /dev/null @@ -1,345 +0,0 @@ -""" -A State that holds an internal TPState, and wraps that by projecting -as needed in from_vector(). We implement to_vector() by extracting -the projected density matrix and then returning -""" -from __future__ import annotations -from typing import TYPE_CHECKING, Literal, Any, Tuple - -if TYPE_CHECKING: - import cvxpy as cp - -import numpy as np -import jax -import jax.numpy as jnp -from pygsti.baseobjs.basis import ExplicitBasis, Basis - - -def strict_triu_op_data(N: int): - sT = N*(N-1)//2 - row_idx, col_idx = jnp.triu_indices(N, k=1) - P_rows = N * row_idx + col_idx - P_cols = np.arange(sT) - return (P_rows, P_cols), (N**2, sT) - - -def jax_vec_to_strict_upper_tri(vec): - ell = vec.shape[0] - n = round( ((8 * ell + 1) ** 0.5 + 1) / 2 ) - if not (n * (n - 1) // 2 == ell): - raise ValueError("The size of the vector must be a triangular number.") - coords, shape = strict_triu_op_data(n) - P_dense = jnp.zeros(shape) - P_dense = P_dense.at[coords].set(1.0) - out = (P_dense @ vec).reshape((n, n), order='F').T - return out - - -def tpvec_to_herm_parts(vec: cp.Expression | jnp.ndarray): - import cvxpy as cp - N2m1 = vec.shape[0] - N = round((N2m1 + 1)**0.5) - block0 = vec[ : N*(N-1)//2 ] - block1 = vec[ N*(N-1)//2 : N*(N-1) ] - block2 = vec[ N*(N-1) : ] - if isinstance(vec, np.ndarray): - vec = jnp.array(vec) - - if isinstance(vec, cp.Expression): - tri_real = cp.vec_to_upper_tri(block0, strict=True) - tri_imag = cp.vec_to_upper_tri(block1, strict=True) - diag = cp.diag(cp.concatenate([block2, 1 - cp.sum(block2)])) - elif isinstance(vec, jnp.ndarray): - tri_real = jax_vec_to_strict_upper_tri(block0) - tri_imag = jax_vec_to_strict_upper_tri(block1) - block3 = jnp.atleast_1d( 1 - jnp.sum(block2, axis=0) ) # type: ignore - diag = jnp.diag( - jnp.concatenate([block2, block3]) # type: ignore - ) - else: - raise ValueError() - - herm_real = tri_real + tri_real.T + diag - herm_imag = tri_imag - tri_imag.T - return herm_real, herm_imag - - -def herm_parts_to_tpvec(herm_parts: tuple[jnp.ndarray, jnp.ndarray] | tuple[Any, Any], full_superket=False): - herm_real, herm_imag = herm_parts - n = herm_imag.shape[0] - triui = np.triu_indices(n, k=1) - flat = (n*(n-1)//2,) + herm_imag.shape[2:] - tri_imag = herm_imag[triui].reshape(flat) - tri_real = herm_real[triui].reshape(flat) - num_diag = n if full_superket else n - 1 - diag = herm_real[np.diag_indices(num_diag)] - if isinstance(herm_real, jnp.ndarray): - return jnp.concatenate([tri_real, tri_imag, diag]) - else: - return np.concatenate([tri_real, tri_imag, diag]) - - -def herm_parts_basis(N: int) -> ExplicitBasis: - triui_rows, triu_cols = np.triu_indices(N, k=1) - sT = N*(N-1)//2 - mats = [] - labels = [] - for k in range(N**2): - mat = np.zeros((N, N), dtype=complex) - if k < sT: - # upper triangle, real, read in row-major - r, c = triui_rows[k], triu_cols[k] - label = '[ReH_%s]' % str({r,c}) - mat[r, c] = 1.0 - mat[c, r] = 1.0 - elif k < 2*sT: - # upper triangle, skew, read in row-major - r, c = triui_rows[k-sT], triu_cols[k-sT] - label = '[ImH_%s]' % str({r,c}) - mat[r, c] = 1.0j - mat[c, r] = -1.0j - else: - i = k - 2*sT - label = '[Diag_{%s}]' % i - mat[i, i] = 1.0 - mats.append(mat) - labels.append(label) - name = f'Hermitian matrix-unit basis of {N**2}-dimensional Hilbert-Schmidt space' - B = ExplicitBasis(mats, labels, name) - return B - - -def eigh_via_symmetricification(herm_parts: tuple[jnp.ndarray, jnp.ndarray] | tuple[Any, Any]) -> tuple[ jnp.ndarray | Any , tuple[Any, Any]]: - """ - Let A be real-symmetric and B be real skew-symmetric, both of order n. - If c is a scalar and x and y are real n-vectors satisfying - - c [x] = [ A B ][x] - [y] = [ B' A ][y]. - - then x - 1j * y is an eigenvector of A + 1j * B with eigenvalue c. - - This function returns a complete set of eigenvalues and eigenvectors - for the matrix herm_parts[0] + 1j * herm_parts[1], where the - eigenvectors are returned in separate matrices containing their real - and imaginary parts. - - The implementation assumes that `eigh` returns eigenpairs ordered by - eigenvalue (either increasing or decreasing is fine). - """ - herm_real, herm_imag = herm_parts - n = herm_imag.shape[0] - if isinstance(herm_real, jnp.ndarray): - sym = jnp.block([[herm_real, herm_imag],[herm_imag.T, herm_real]]) - sevals, sevecs = jnp.linalg.eigh(sym) - else: - sym = np.block([[herm_real, herm_imag],[herm_imag.T, herm_real]]) - sevals, sevecs = np.linalg.eigh(sym) - evals = sevals[::2] - evecs_real = sevecs[ :n, ::2 ] - evecs_imag = -sevecs[ n:, ::2 ] - return evals, (evecs_real, evecs_imag) - - -def recollect_expanded_eigendecomposition(evals, evecs_real, evecs_imag): - """ - evals, evecs_real, and evecs_imag are all real-valued. - - Define D = diag(evals) and U = evecs_real + 1j * evecs_imag. - - This function returns matrices herm_real and herm_imag where - - herm_real + 1j * herm_imag = U @ D @ U.T.conj(). - """ - temp_real = evals[:, None] * evecs_real.T - temp_imag = evals[:, None] * evecs_imag.T - herm_real = evecs_real @ temp_real + evecs_imag @ temp_imag - herm_imag = evecs_imag @ temp_real - evecs_real @ temp_imag - return herm_real, herm_imag - - -def density_projection_model(N: int, param_name='parameters', var_name='projection'): - import cvxpy as cp - param = cp.Parameter(shape=(N**2 - 1,), name=param_name) - herm_real, herm_imag = tpvec_to_herm_parts(param) - herm_expr = herm_real + 1j * herm_imag - - X = cp.Variable(shape=(N, N), hermitian=True, name=var_name) - F = cp.sum_squares(herm_expr - X) - C = [ 0 << X, cp.trace(cp.real(X)) == 1 ] - problem = cp.Problem( cp.Minimize(F), C ) - - """ - To map V := X.value back to `param` space, - set triui = np.triu_indices_from(V, k=1), - then concatenate ... - V.real[triui], - V.imag[triui], and - np.diag(V.real)[:-1]. - """ - return problem - - -def project_onto_simplex(vec: jnp.ndarray) -> jnp.ndarray: - # Sort in descending order - u = jnp.sort(vec)[::-1] - css = jnp.cumsum(u, axis=0) - # Find the largest index rho such that - # u[rho] > (css[rho] - 1) / (rho + 1) - inds = jnp.arange(vec.shape[0]) + 1 - cond = u - (css - 1) / inds > 0 - if not jnp.any(cond): - # if no positive entries, project uniformly on first coordinate - theta = (css[0] - 1) / 1 - else: - rho = jnp.where(cond)[0][-1] - theta = (css[rho] - 1) / (rho + 1) - proj = jnp.maximum(vec - theta, 0) - return proj - - -def project_onto_densities(herm_mx: jnp.ndarray) -> jnp.ndarray: - v, u = jnp.linalg.eigh(herm_mx) - udag = u.T.conj() - proj_v = project_onto_simplex(v) - density_mx = u @ (proj_v[:, None] * udag) - return density_mx - - -def project_onto_densities_realification(herm_parts: jnp.ndarray) -> jnp.ndarray: - """ - Example usage: - - jac = jax.jacobian(project_onto_densities_realification) - - herm_mx = np.random.randn(3,3) + 1j * np.random.randn(3,3) - herm_mx += temp.T.conj() - - herm_mx_reif = jnp.stack((herm_mx.real, temp.imag)) - # ^ shape (2, 3, 3) - out = jac(herm_mx_reif) - # ^ shape (2, 3, 3, 2, 3, 3) - - """ - evals, evecs = eigh_via_symmetricification(herm_parts) # type: ignore - proj_evals = project_onto_simplex(evals) - density_mx = recollect_expanded_eigendecomposition(proj_evals, *evecs) - density_mx = jnp.stack(density_mx) - return density_mx - - -def project_onto_densities_tpvec( tpvec, full_superket=False ): - herm_parts = tpvec_to_herm_parts( tpvec ) - herm_mx = jnp.stack(herm_parts) # type: ignore - density_mx = project_onto_densities_realification(herm_mx) - proj_parts = (density_mx[0], density_mx[1]) - tpvec_proj = herm_parts_to_tpvec(proj_parts, full_superket) - return tpvec_proj - - - -class ProjectedState_CVXPY: - - def __init__(self, N: int) -> None: - problem = density_projection_model(N) - self._cvxpy_model = problem - self._cvxpy_param = problem.param_dict['parameters'] - return - -from pygsti.modelmembers import ModelMember -from pygsti.modelmembers.states import State, TPState -from pygsti.tools.basistools import vec_to_stdmx, change_basis - - -# TODO: need a routine for changing basis of to_vector()/from_vector() with TPState -# versus ProjectedState. - - - - - -class ProjectedState(TPState): - - def __init__(self, vec: np.ndarray, basis: Basis, evotype="default", state_space=None): - super().__init__(vec, basis, evotype, state_space) - self._working_basis = herm_parts_basis(self.state_space.udim) - - self._n2w_np = self._basis.create_transform_matrix(self._working_basis).real - self._w2n_np = self._working_basis.create_transform_matrix(self._basis).real - self._n2w = jnp.array( self._n2w_np ) - self._w2n = jnp.array( self._w2n_np ) - - self._natural_firstEl = np.array([self.basis.elsize ** -0.25]) - return - - def to_vector(self) -> np.ndarray: - return super().to_vector() - - def _natural_vec_to_working_vec(self, natvec): - if isinstance(natvec, np.ndarray): - workvec = self._n2w_np[:-1, :] @ np.concatenate([self._natural_firstEl, natvec]) - else: - workvec = self._n2w[:-1, :] @ jnp.concatenate([self._natural_firstEl, natvec]) - return workvec - - def _working_vec_to_natural_vec(self, workvec): - if isinstance(workvec, np.ndarray): - tail = np.atleast_1d(1 - np.sum(workvec[-(self.state_space.udim - 1):])) - natvec = self._w2n_np[1:, :] @ np.concatenate([workvec, tail]) - else: - tail = jnp.atleast_1d(1 - jnp.sum(workvec[-(self.state_space.udim - 1):])) - natvec = self._w2n[1:, :] @ jnp.concatenate([workvec, tail]) - return natvec - - def from_vector(self, v, close=False, dirty_value=True) -> None: - v = jnp.array(v) - tpvec_before = self._natural_vec_to_working_vec(v) - tpvec_after = project_onto_densities_tpvec( tpvec_before ) - v_projected = self._working_vec_to_natural_vec( tpvec_after ) - v_projected = np.asarray(v_projected).astype(np.double) - super().from_vector(v_projected, False, True) - return - - def set_dense(self, vec): - super().set_dense(vec) - self.from_vector(self.to_vector(), False, True) - return - - def deriv_wrt_params(self, wrt_filter=None): - assert wrt_filter is None - - def projected_tpvec_onto_densities_as_superket( tpvec ): - todense_working = project_onto_densities_tpvec( tpvec, full_superket=True ) - todense_natural = self._w2n @ todense_working - return todense_natural - - jac_fun = jax.jacobian(projected_tpvec_onto_densities_as_superket) - natural_sket = jnp.array(self.to_dense('HilbertSchmidt')) - current_tpvec = self._n2w[:-1, :] @ natural_sket - jac_val = jac_fun( current_tpvec ) - jac_val = jac_val @ self._n2w[:-1, :] - # ^ Jacobian mapping natural_sket_in to natural_sket_out. - jac_val = jac_val[:, 1:] - # ^ Jacobian mapping the last `from_vector` arg `v` to self.to_dense(). - - jac_val = np.asarray( jac_val, dtype=np.double ) - return jac_val - - def stateless_data(self) -> Tuple[int]: - raise NotImplementedError() - - @staticmethod - def torch_base(sd, t_param): - raise NotImplementedError() - - def has_nonzero_hessian(self): - return False - - -if __name__ == '__main__': - temp = np.random.randn(3,3) + 1j * np.random.randn(3,3) - temp += temp.T.conj() - tpvec = np.random.rand(3) - jac2 = jax.jacobian(project_onto_densities_tpvec) - x,y = tpvec_to_herm_parts(jnp.array(tpvec)) - print() From 3af6652a309d86fe75c57f59f23235a04f808762 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Sat, 14 Mar 2026 08:29:58 -0700 Subject: [PATCH 20/21] move helper operation classes for CPTP instruments into modelmembers/operations/cptrop.py. Add tests for those helper classes in test/unit/modelmembers/instruments.py. Move a TPInstrumentOp test class from test/unit/modelmembers/operations.py into the new test/unit/modelmembers/instrument.py. Docs for pygsti/modelmembers/instruments::cptp_instrument(...) --- pygsti/modelmembers/instruments/__init__.py | 54 ++++- pygsti/modelmembers/operations/__init__.py | 2 + .../cptrop.py} | 116 +++++++--- pygsti/tools/optools.py | 2 +- pygsti/tools/sdptools.py | 2 +- test/unit/modelmembers/test_instrument.py | 211 ++++++++++++++++++ test/unit/modelmembers/test_operation.py | 33 +-- test/unit/objects/test_instrument.py | 74 ------ 8 files changed, 353 insertions(+), 141 deletions(-) rename pygsti/modelmembers/{instruments/cptpinstrument.py => operations/cptrop.py} (54%) create mode 100644 test/unit/modelmembers/test_instrument.py delete mode 100644 test/unit/objects/test_instrument.py diff --git a/pygsti/modelmembers/instruments/__init__.py b/pygsti/modelmembers/instruments/__init__.py index ff7b2ae44..e3331013b 100644 --- a/pygsti/modelmembers/instruments/__init__.py +++ b/pygsti/modelmembers/instruments/__init__.py @@ -13,7 +13,7 @@ from .instrument import Instrument from .tpinstrument import TPInstrument from .tpinstrumentop import TPInstrumentOp -from .cptpinstrument import RootConjOperator, SummedOperator +from ..operations.cptrop import RootConjOperator, SummedOperator import warnings as _warnings import scipy.linalg as _la @@ -165,7 +165,57 @@ def convert(instrument, to_type, basis, ideal_instrument=None, flatten_structure raise ValueError("Cannot convert an instrument to type %s" % to_type) -def cptp_instrument(op_arrays: dict[str, _np.ndarray], basis: Basis, error_tol:float=1e-6, trunc_tol:float=1e-7, povm_errormap=None) -> Instrument: +def cptp_instrument(op_arrays: dict[str, _np.ndarray], basis: Basis, + error_tol: float = 1e-6, trunc_tol: float = 1e-7, + povm_errormap=None) -> Instrument: + """ + Construct a CPTPLND-parameterized instrument from a dictionary of dense superoperator arrays. + + Each operator in `op_arrays` must be a completely positive trace-reducing (CPTR) map. + The function Kraus-decomposes each operator, polar-decomposes each Kraus operator into + a unitary and a positive semidefinite part, then assembles the result using + CPTPLND-parameterized unitaries and a shared `ComposedPOVM` that enforces the + completeness (trace-preservation) constraint across all instrument outcomes. + + The construction follows these steps: + + 1. For each outcome label, compute a minimal Kraus decomposition and polar-decompose + each Kraus operator ``K = u p`` (``p`` PSD, ``u`` unitary). + 2. Represent each Kraus term as the composition of a ``RootConjOperator`` (encoding + ``ρ ↦ p² ρ p²``, parameterized by the POVM effect for ``p²``) and a + CPTPLND-parameterized unitary channel. + 3. Collect all POVM effects into a shared ``ComposedPOVM`` whose error map is + constrained to be CPTPLND. This enforces that the effects sum to ≤ I, which + ensures the full instrument is CPTR. + 4. If an outcome requires more than one Kraus term, wrap the summands in a + ``SummedOperator``; otherwise use the single term directly. + + Parameters + ---------- + op_arrays : dict of str → numpy.ndarray + Mapping from outcome label to dense superoperator matrix (in the given basis). + Each matrix must represent a CPTR map; the function raises or warns if the + Kraus decomposition is inconsistent with this. + basis : Basis or str + The operator basis in which the superoperators are expressed (e.g. ``'pp'``). + error_tol : float, optional + Tolerance for eigenvalue errors in the minimal Kraus decomposition. Eigenvalues + below ``-error_tol`` cause an error to be raised. Default ``1e-6``. + trunc_tol : float, optional + Truncation tolerance for the minimal Kraus decomposition. Eigenvalues between + ``-error_tol`` and ``trunc_tol`` are silently set to zero. Default ``1e-7``. + povm_errormap : LinearOperator, optional + A CPTPLND-parameterized error map to use as the shared POVM error map. + If ``None`` (default), the identity channel is used as the starting point + and converted to CPTPLND parameterization. + + Returns + ------- + Instrument + An ``Instrument`` whose member operations are CPTR maps sharing a common set + of CPTPLND parameters. The instrument is not itself a ``TPInstrument`` because + the completeness constraint is enforced implicitly through the shared POVM. + """ unitaries = dict() effects = dict() diff --git a/pygsti/modelmembers/operations/__init__.py b/pygsti/modelmembers/operations/__init__.py index d0f340c4f..e585e36ea 100644 --- a/pygsti/modelmembers/operations/__init__.py +++ b/pygsti/modelmembers/operations/__init__.py @@ -39,11 +39,13 @@ from .stochasticop import StochasticNoiseOp from .lindbladcoefficients import LindbladCoefficientBlock as _LindbladCoefficientBlock from .affineshiftop import AffineShiftOp +from .cptrop import RootConjOperator, SummedOperator from pygsti.baseobjs import statespace as _statespace from pygsti.tools import basistools as _bt from pygsti.tools import optools as _ot from pygsti import SpaceT + def create_from_unitary_mx(unitary_mx, op_type, basis='pp', stdname=None, evotype='default', state_space=None): """ TODO: docstring - note that op_type can be a list/tuple of types in order of precedence """ op_type_preferences = verbose_type_from_op_type(op_type) diff --git a/pygsti/modelmembers/instruments/cptpinstrument.py b/pygsti/modelmembers/operations/cptrop.py similarity index 54% rename from pygsti/modelmembers/instruments/cptpinstrument.py rename to pygsti/modelmembers/operations/cptrop.py index 7f7988c54..886ad30e6 100644 --- a/pygsti/modelmembers/instruments/cptpinstrument.py +++ b/pygsti/modelmembers/operations/cptrop.py @@ -1,5 +1,14 @@ """ -Defines the CPTPInstrument class +Operators used in CPTR (completely positive trace-reducing) parameterizations. + +A completely positive trace-reducing (CPTR) map is a completely positive map +whose Kraus operators satisfy sum_i K_i^† K_i ≤ I (i.e., trace is non-increasing). +Every CPTR map can be decomposed as a CPTP map followed by a "root-conjugation" +by a positive semidefinite matrix E with 0 ≤ E ≤ I. + +The two classes here, `RootConjOperator` and `SummedOperator`, are `LinearOperator` +subclasses used as building blocks when constructing CPTPLND-parameterized instruments +via :func:`pygsti.modelmembers.instruments.cptp_instrument`. """ #*************************************************************************************************** # Copyright 2015, 2019, 2025 National Technology & Engineering Solutions of Sandia, LLC (NTESS). @@ -14,7 +23,7 @@ from pygsti.pgtypes import SpaceT from pygsti.baseobjs.basis import Basis from pygsti.modelmembers.povms.effect import POVMEffect -from pygsti.modelmembers.operations import LinearOperator +from pygsti.modelmembers.operations.linearop import LinearOperator from pygsti.tools import optools as _ot from typing import Union @@ -23,9 +32,31 @@ class RootConjOperator(LinearOperator): """ - A linear operator parameterized by a matrix E where 0 ≤ E ≤ 1, whose action takes ρ to E½ ρ E½ . - - Every CPTR map can be obtained by pre-composing a CPTP map with this kind of linear operator. + A linear operator parameterized by a PSD matrix E with 0 ≤ E ≤ I, acting as ρ ↦ E^½ ρ E^½. + + Every CPTR (completely positive trace-reducing) map can be expressed as the + composition of a CPTP map with a ``RootConjOperator``. This class implements + the "root-conjugation" factor: given a POVM effect E (representing the PSD matrix + in superket form), the superoperator for ρ ↦ E^½ ρ E^½ is constructed and kept + up to date as the effect's parameters change. + + The effect E encodes the constraint 0 ≤ E ≤ I; its parameters are inherited + directly as the parameters of this operator. + + Parameters + ---------- + effect : POVMEffect + A POVM effect whose superket encodes the PSD matrix E. The effect's + parameter vector is shared (via gpindices) with this operator. + basis : Basis or str + The operator basis in which the superoperator is represented. + + Attributes + ---------- + EIGTOL_WARNING : float + Eigenvalue tolerance below which a warning is issued during construction. + EIGTOL_ERROR : float + Eigenvalue tolerance below which an error is raised during construction. """ EIGTOL_WARNING = 1e-10 @@ -38,16 +69,18 @@ def __init__(self, effect: POVMEffect, basis: BasisLike): self._evotype = effect.evotype dim = self._state_space.dim - self._rep = self._evotype.create_dense_superop_rep( _np.zeros((dim, dim)), self._basis, self._state_space ) + self._rep = self._evotype.create_dense_superop_rep( + _np.zeros((dim, dim)), self._basis, self._state_space + ) self._update_rep_base() LinearOperator.__init__(self, self._rep, self._evotype) self.init_gpindices() - + def submembers(self): return [self._effect] - + def _update_rep_base(self): - # This function is directly analogous to TPInstrumentOp._construct_matrix. + # Analogous to TPInstrumentOp._construct_matrix. self._rep.base.flags.writeable = True assert(self._rep.base.shape == (self.dim, self.dim)) effect_superket = self._effect.to_dense() @@ -55,21 +88,19 @@ def _update_rep_base(self): self._rep.base[:] = mx self._rep.base.flags.writeable = False self._rep.base_has_changed() - return - + def deriv_wrt_params(self, wrt_filter=None): return LinearOperator.deriv_wrt_params(self, wrt_filter) - + def has_nonzero_hessian(self): - # This is not affine in its parameters. + # Not affine in its parameters. return True def from_vector(self, v, close=False, dirty_value=True): for sm, local_inds in zip(self.submembers(), self._submember_rpindices): sm.from_vector(v[local_inds], close, dirty_value) self._update_rep_base() - return - + @property def num_params(self): return len(self.gpindices_as_array()) @@ -88,7 +119,33 @@ def to_dense(self, on_space: SpaceT = 'HilbertSchmidt') -> _np.ndarray: class SummedOperator(LinearOperator): - + """ + A linear operator whose superoperator is the sum of a list of constituent operators. + + ``SummedOperator`` wraps a sequence of ``LinearOperator`` objects and exposes + their pointwise sum as a single ``LinearOperator``. It is used inside + :func:`~pygsti.modelmembers.instruments.cptp_instrument` when a CPTR map requires + more than one Kraus term: each term contributes a separate summand, and the full + map is their sum. + + The parameter vector is the concatenation (with shared-index deduplication via + gpindices) of the submembers' parameter vectors. The constituent operators' ``_rep`` + objects are linked directly, so calling ``from_vector`` on a submember automatically + updates this operator's representation. + + Parameters + ---------- + operators : list of LinearOperator + The operators to sum. All must share the same evotype, state space, and + superoperator dimension. + basis : Basis or str + The operator basis (used for bookkeeping; the operators themselves must + already be expressed in this basis). + + Notes + ----- + ``deriv_wrt_params`` is not implemented for this class. + """ def __init__(self, operators, basis: BasisLike): op = operators[0] @@ -97,15 +154,13 @@ def __init__(self, operators, basis: BasisLike): self._state_space = op.state_space self._evotype = op.evotype self._subreps = [op._rep for op in self._operators] - self._rep = self._evotype.create_sum_rep( self._subreps, self._state_space ) + self._rep = self._evotype.create_sum_rep(self._subreps, self._state_space) LinearOperator.__init__(self, self._rep, self._evotype) self.init_gpindices() - # NOTE: This class doesn't have a function analogous to _update_rep_base - # that we use in RootConjOperator. We can get away with not having such - # a function because it's the responsibility of op.from_vector(...) - # to update op's attached OpRep. - return - + # NOTE: No _update_rep_base analogue is needed here. Each constituent + # operator's from_vector(...) updates its own attached OpRep, and the + # sum rep reflects those changes automatically. + def submembers(self): out = [] hit = set() @@ -116,18 +171,17 @@ def submembers(self): hit.add(id(temp)) out.append(sm) return out - + def deriv_wrt_params(self, wrt_filter=None): raise NotImplementedError() - + def has_nonzero_hessian(self): - return any(op.has_nonzero_hession() for op in self._operators) + return any(op.has_nonzero_hessian() for op in self._operators) def from_vector(self, v, close=False, dirty_value=True): for sm, local_inds in zip(self.submembers(), self._submember_rpindices): sm.from_vector(v[local_inds], close, dirty_value) - return - + @property def num_params(self): return len(self.gpindices_as_array()) @@ -140,11 +194,9 @@ def to_vector(self): def to_dense(self, on_space: SpaceT = 'HilbertSchmidt') -> _np.ndarray: assert on_space in ('HilbertSchmidt', 'minimal') - on_space = 'HilbertSchmidt' - out = self._operators[0].to_dense(on_space) + out = self._operators[0].to_dense('HilbertSchmidt') if not out.flags.writeable: out = out.copy() for op in self._operators[1:]: - temp = op.to_dense(on_space) - out += temp + out += op.to_dense('HilbertSchmidt') return out diff --git a/pygsti/tools/optools.py b/pygsti/tools/optools.py index 212f1370d..12b457f44 100644 --- a/pygsti/tools/optools.py +++ b/pygsti/tools/optools.py @@ -1210,7 +1210,7 @@ def povm_diamonddist(model, target_model, povmlbl, _premultiplier=None): return diamonddist(povm_mx, target_povm_mx, target_model.basis) except AssertionError as e: assert '`dim` must be a perfect square' in str(e) - return _np.NaN + return _np.nan def instrument_infidelity(a, b, mx_basis): diff --git a/pygsti/tools/sdptools.py b/pygsti/tools/sdptools.py index 704aa463b..40d7fc4e9 100644 --- a/pygsti/tools/sdptools.py +++ b/pygsti/tools/sdptools.py @@ -42,7 +42,7 @@ def solve_sdp(prob: cp.Problem, **kwargs) -> tuple[np.floating, dict[str, np.ndarray]]: - objective_val : np.floating = np.array(np.NaN).item() + objective_val : np.floating = np.array(np.nan).item() varvals : dict[str, np.ndarray] = dict() for i, solver in enumerate(SDP_SOLVER_PRIORITY): try: diff --git a/test/unit/modelmembers/test_instrument.py b/test/unit/modelmembers/test_instrument.py new file mode 100644 index 000000000..f79bb72df --- /dev/null +++ b/test/unit/modelmembers/test_instrument.py @@ -0,0 +1,211 @@ +import numpy as np + +import pygsti.modelmembers.operations as op +from pygsti.modelmembers import instruments as inst +from pygsti.modelmembers.instruments import cptp_instrument +from pygsti.modelmembers.instruments import Instrument, TPInstrument +from pygsti.modelmembers.operations import RootConjOperator, SummedOperator +from pygsti.modelmembers import povms as pv +from pygsti.modelpacks.legacy import std1Q_XYI as std +from pygsti.models.gaugegroup import FullGaugeGroupElement +from ..util import BaseCase +from .test_operation import ImmutableDenseOpBase + + +class InstrumentMethodBase(object): + def test_num_elements(self): + self.assertEqual(self.instrument.num_elements, self.n_elements) + + def test_copy(self): + inst_copy = self.instrument.copy() + # TODO assert correctness + + def test_to_string(self): + inst_str = str(self.instrument) + # TODO assert correctness + + def test_transform(self): + T = FullGaugeGroupElement( + np.array([[1, 0, 0, 0], + [0, 1, 0, 0], + [0, 0, 0, 1], + [0, 0, 1, 0]], 'd')) + self.instrument.transform_inplace(T) + # TODO assert correctness + + def test_simplify_operations(self): + gates = self.instrument.simplify_operations(prefix="ABC") + # TODO assert correctness + + def test_constructor_raises_on_non_none_param_conflict(self): + with self.assertRaises(AssertionError): + self.constructor(["Non-none-matrices"], 'default', None, False, ["Non-none-items"]) # can't both be non-None + + def test_constructor_raises_on_bad_op_matrices_type(self): + with self.assertRaises(ValueError): + self.constructor("foobar") # op_matrices must be a list or dict + + def test_convert_raises_on_unknown_basis(self): + with self.assertRaises(ValueError): + inst.convert(self.instrument, "foobar", self.model.basis) + + +class InstrumentInstanceBase(object): + def setUp(self): + # Initialize standard target model for instruments + # XXX can instruments be tested independently of a model? EGN: yes, I was just lazy; but they should also be tested within a model. + self.n_elements = 32 + + self.model = std.target_model() + E = self.model.povms['Mdefault']['0'] + Erem = self.model.povms['Mdefault']['1'] + self.Gmz_plus = np.dot(E, E.T) + self.Gmz_minus = np.dot(Erem, Erem.T) + # XXX is this used? + self.povm_ident = self.model.povms['Mdefault']['0'] + self.model.povms['Mdefault']['1'] + self.instrument = self.constructor({'plus': self.Gmz_plus, 'minus': self.Gmz_minus}) + self.model.instruments['Iz'] = self.instrument + super(InstrumentInstanceBase, self).setUp() + + +class InstrumentInstanceTester(InstrumentMethodBase, InstrumentInstanceBase, BaseCase): + constructor = inst.Instrument + + +class TPInstrumentInstanceTester(InstrumentMethodBase, InstrumentInstanceBase, BaseCase): + constructor = inst.TPInstrument + + def test_raise_on_modify(self): + with self.assertRaises(ValueError): + self.instrument['plus'] = None # can't set value of a TP Instrument element + + +class CPTPInstrumentTester(BaseCase): + """Tests for cptp_instrument and the CPTR operation primitives.""" + + def setUp(self): + model = std.target_model() + E = model.povms['Mdefault']['0'] + Erem = model.povms['Mdefault']['1'] + self.Gmz_plus = np.dot(E, E.T) + self.Gmz_minus = np.dot(Erem, Erem.T) + self.basis = model.basis + self.op_arrays = {'plus': self.Gmz_plus, 'minus': self.Gmz_minus} + self.instrument = cptp_instrument(self.op_arrays, self.basis) + + # --- cptp_instrument factory tests --- + + def test_construction_returns_instrument(self): + self.assertIsInstance(self.instrument, Instrument) + + def test_construction_keys(self): + self.assertSetEqual(set(self.instrument.keys()), {'plus', 'minus'}) + + def test_to_from_vector_roundtrip(self): + v0 = self.instrument.to_vector() + self.assertEqual(len(v0), self.instrument.num_params) + # Perturb and restore; dense output should change then return. + v_perturbed = v0 + 1e-3 * np.random.default_rng(0).standard_normal(len(v0)) + self.instrument.from_vector(v_perturbed) + v1 = self.instrument.to_vector() + self.assertFalse(np.allclose(v0, v1)) + self.instrument.from_vector(v0) + v2 = self.instrument.to_vector() + np.testing.assert_allclose(v0, v2, atol=1e-12) + + def test_convert_to_cptplnd(self): + plain = inst.Instrument({'plus': self.Gmz_plus, 'minus': self.Gmz_minus}) + converted = inst.convert(plain, 'CPTPLND', self.basis) + self.assertIsInstance(converted, Instrument) + self.assertSetEqual(set(converted.keys()), {'plus', 'minus'}) + + def test_member_dense_shapes(self): + dim = self.basis.dim + for op in self.instrument.values(): + mx = op.to_dense('HilbertSchmidt') + self.assertEqual(mx.shape, (dim, dim)) + + def test_sum_of_members_is_cptp_like(self): + # The sum of all instrument members' superoperators should equal the + # superoperator of the total channel (which is CPTP), i.e. the first + # row of the sum must be [1, 0, 0, 0] in the pp basis. + total = sum(op.to_dense('HilbertSchmidt') for op in self.instrument.values()) + # First row encodes trace preservation: tr(rho) is preserved. + np.testing.assert_allclose(total[0, 0], 1.0, atol=1e-6) + np.testing.assert_allclose(total[0, 1:], 0.0, atol=1e-6) + + # --- RootConjOperator unit tests --- + + def test_root_conj_operator_to_dense_shape(self): + # Build a real StaticPOVMEffect from the plus POVM effect + effect = pv.StaticPOVMEffect(self.Gmz_plus[:, 0]) + rcop = RootConjOperator(effect, self.basis) + mx = rcop.to_dense('HilbertSchmidt') + dim = self.basis.dim + self.assertEqual(mx.shape, (dim, dim)) + + def test_root_conj_operator_num_params(self): + effect = pv.StaticPOVMEffect(self.Gmz_plus[:, 0]) + rcop = RootConjOperator(effect, self.basis) + # StaticPOVMEffect has zero parameters; RootConjOperator inherits them. + self.assertEqual(rcop.num_params, 0) + + def test_root_conj_operator_has_nonzero_hessian(self): + effect = pv.StaticPOVMEffect(self.Gmz_plus[:, 0]) + rcop = RootConjOperator(effect, self.basis) + self.assertTrue(rcop.has_nonzero_hessian()) + + # --- SummedOperator unit tests --- + + def test_summed_operator_to_dense_equals_sum(self): + e1 = pv.StaticPOVMEffect(self.Gmz_plus[:, 0]) + e2 = pv.StaticPOVMEffect(self.Gmz_minus[:, 0]) + r1 = RootConjOperator(e1, self.basis) + r2 = RootConjOperator(e2, self.basis) + sop = SummedOperator([r1, r2], self.basis) + expected = r1.to_dense() + r2.to_dense() + np.testing.assert_allclose(sop.to_dense(), expected, atol=1e-12) + + def test_summed_operator_num_params(self): + e1 = pv.StaticPOVMEffect(self.Gmz_plus[:, 0]) + e2 = pv.StaticPOVMEffect(self.Gmz_minus[:, 0]) + r1 = RootConjOperator(e1, self.basis) + r2 = RootConjOperator(e2, self.basis) + sop = SummedOperator([r1, r2], self.basis) + self.assertEqual(sop.num_params, 0) + + def test_summed_operator_deriv_raises(self): + e1 = pv.StaticPOVMEffect(self.Gmz_plus[:, 0]) + r1 = RootConjOperator(e1, self.basis) + sop = SummedOperator([r1], self.basis) + with self.assertRaises(NotImplementedError): + sop.deriv_wrt_params() + + +class TPInstrumentOpTester(ImmutableDenseOpBase, BaseCase): + n_params = 28 + + @staticmethod + def build_gate(): + Gmz_plus = np.array([[0.5, 0, 0, 0.5], + [0, 0, 0, 0], + [0, 0, 0, 0], + [0.5, 0, 0, 0.5]]) + Gmz_minus = np.array([[0.5, 0, 0, -0.5], + [0, 0, 0, 0], + [0, 0, 0, 0], + [-0.5, 0, 0, 0.5]]) + evotype = 'default' + instrument = TPInstrument({ + 'plus': op.FullArbitraryOp(Gmz_plus, 'pp', evotype), + 'minus': op.FullArbitraryOp(Gmz_minus, 'pp', evotype) + }) + return instrument['plus'] + + def test_vector_conversion(self): + self.gate.to_vector() # now to_vector is allowed + + def test_deriv_wrt_params_shape(self): + super(TPInstrumentOpTester, self).test_deriv_wrt_params() + deriv = self.gate.deriv_wrt_params([0]) + self.assertEqual(deriv.shape[1], 1) diff --git a/test/unit/modelmembers/test_operation.py b/test/unit/modelmembers/test_operation.py index a9861d951..a0990eedb 100644 --- a/test/unit/modelmembers/test_operation.py +++ b/test/unit/modelmembers/test_operation.py @@ -10,7 +10,6 @@ import pygsti.tools.optools as gt from pygsti.models.modelconstruction import create_spam_vector, create_operation from pygsti.evotypes import Evotype -from pygsti.modelmembers.instruments import TPInstrument from pygsti.modelmembers.states import FullState from pygsti.models import ExplicitOpModel from pygsti.baseobjs import statespace, basisconstructors as bc @@ -219,6 +218,7 @@ def test_rotate(self): self.gate.rotate([0.01, 0.02, 0.03], 'gm') # TODO assert correctness + class ImmutableDenseOpBase(DenseOpBase): def test_raises_on_set_value(self): M = np.asarray(self.gate) # gate as a matrix @@ -430,6 +430,7 @@ def test_first_row_read_only(self): with self.assertRaises(ValueError): self.gate[0][1:2] = [0] + class AffineShiftOpTester(DenseOpBase, BaseCase): n_params = 3 @@ -773,36 +774,6 @@ def test_constructor_raises_on_state_space_label_mismatch(self): op.EmbeddedOp(state_space, ['Q0', 'Q1'], op.FullArbitraryOp(mx, evotype=evotype, state_space=None)) -class TPInstrumentOpTester(ImmutableDenseOpBase, BaseCase): - n_params = 28 - - @staticmethod - def build_gate(): - # XXX can this be constructed directly? EGN: what do you mean? - Gmz_plus = np.array([[0.5, 0, 0, 0.5], - [0, 0, 0, 0], - [0, 0, 0, 0], - [0.5, 0, 0, 0.5]]) - Gmz_minus = np.array([[0.5, 0, 0, -0.5], - [0, 0, 0, 0], - [0, 0, 0, 0], - [-0.5, 0, 0, 0.5]]) - evotype = 'default' - inst = TPInstrument({'plus': op.FullArbitraryOp(Gmz_plus, evotype=evotype), 'minus': op.FullArbitraryOp( - Gmz_minus, evotype=evotype)}) - return inst['plus'] - - def test_vector_conversion(self): - self.gate.to_vector() # now to_vector is allowed - - def test_deriv_wrt_params(self): - super(TPInstrumentOpTester, self).test_deriv_wrt_params() - - # XXX does this check anything meaningful? EGN: yes, this checks that when I give deriv_wrt_params a length-1 list it's return value has the right shape. - deriv = self.gate.deriv_wrt_params([0]) - self.assertEqual(deriv.shape[1], 1) - - class StochasticNoiseOpTester(BaseCase): def test_instance(self): state_space = statespace.default_space_for_dim(4) diff --git a/test/unit/objects/test_instrument.py b/test/unit/objects/test_instrument.py deleted file mode 100644 index b11850a5c..000000000 --- a/test/unit/objects/test_instrument.py +++ /dev/null @@ -1,74 +0,0 @@ -import numpy as np - -from pygsti.modelmembers import instruments as inst -from pygsti.modelpacks.legacy import std1Q_XYI as std -from pygsti.models.gaugegroup import FullGaugeGroupElement -from ..util import BaseCase - - -class InstrumentMethodBase(object): - def test_num_elements(self): - self.assertEqual(self.instrument.num_elements, self.n_elements) - - def test_copy(self): - inst_copy = self.instrument.copy() - # TODO assert correctness - - def test_to_string(self): - inst_str = str(self.instrument) - # TODO assert correctness - - def test_transform(self): - T = FullGaugeGroupElement( - np.array([[1, 0, 0, 0], - [0, 1, 0, 0], - [0, 0, 0, 1], - [0, 0, 1, 0]], 'd')) - self.instrument.transform_inplace(T) - # TODO assert correctness - - def test_simplify_operations(self): - gates = self.instrument.simplify_operations(prefix="ABC") - # TODO assert correctness - - def test_constructor_raises_on_non_none_param_conflict(self): - with self.assertRaises(AssertionError): - self.constructor(["Non-none-matrices"], 'default', None, False, ["Non-none-items"]) # can't both be non-None - - def test_constructor_raises_on_bad_op_matrices_type(self): - with self.assertRaises(ValueError): - self.constructor("foobar") # op_matrices must be a list or dict - - def test_convert_raises_on_unknown_basis(self): - with self.assertRaises(ValueError): - inst.convert(self.instrument, "foobar", self.model.basis) - - -class InstrumentInstanceBase(object): - def setUp(self): - # Initialize standard target model for instruments - # XXX can instruments be tested independently of a model? EGN: yes, I was just lazy; but they should also be tested within a model. - self.n_elements = 32 - - self.model = std.target_model() - E = self.model.povms['Mdefault']['0'] - Erem = self.model.povms['Mdefault']['1'] - self.Gmz_plus = np.dot(E, E.T) - self.Gmz_minus = np.dot(Erem, Erem.T) - # XXX is this used? - self.povm_ident = self.model.povms['Mdefault']['0'] + self.model.povms['Mdefault']['1'] - self.instrument = self.constructor({'plus': self.Gmz_plus, 'minus': self.Gmz_minus}) - self.model.instruments['Iz'] = self.instrument - super(InstrumentInstanceBase, self).setUp() - - -class InstrumentInstanceTester(InstrumentMethodBase, InstrumentInstanceBase, BaseCase): - constructor = inst.Instrument - - -class TPInstrumentInstanceTester(InstrumentMethodBase, InstrumentInstanceBase, BaseCase): - constructor = inst.TPInstrument - - def test_raise_on_modify(self): - with self.assertRaises(ValueError): - self.instrument['plus'] = None # can't set value of a TP Instrument element From 424081ff5765a0d0c58e34be1c2e7eb5e58ec9df Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Sat, 14 Mar 2026 09:12:43 -0700 Subject: [PATCH 21/21] improve tests --- test/unit/modelmembers/test_instrument.py | 83 +++++++++++++---------- 1 file changed, 47 insertions(+), 36 deletions(-) diff --git a/test/unit/modelmembers/test_instrument.py b/test/unit/modelmembers/test_instrument.py index f79bb72df..2ae9f4329 100644 --- a/test/unit/modelmembers/test_instrument.py +++ b/test/unit/modelmembers/test_instrument.py @@ -12,17 +12,43 @@ from .test_operation import ImmutableDenseOpBase -class InstrumentMethodBase(object): +class InstrumentTestBase(BaseCase): + """Shared setUp and test methods for Instrument and TPInstrument. + + Subclasses must define a class attribute ``constructor`` pointing to the + instrument class under test. + """ + __test__ = False + constructor: type + + def setUp(self): + self.n_elements = 32 + self.model = std.target_model() + E = self.model.povms['Mdefault']['0'] + Erem = self.model.povms['Mdefault']['1'] + self.Gmz_plus = np.dot(E, E.T) + self.Gmz_minus = np.dot(Erem, Erem.T) + self.instrument: Instrument = self.constructor({'plus': self.Gmz_plus, 'minus': self.Gmz_minus}) + self.model.instruments['Iz'] = self.instrument + def test_num_elements(self): self.assertEqual(self.instrument.num_elements, self.n_elements) def test_copy(self): inst_copy = self.instrument.copy() - # TODO assert correctness + self.assertIsInstance(inst_copy, type(self.instrument)) + self.assertEqual(list(inst_copy.keys()), list(self.instrument.keys())) + for key in self.instrument.keys(): + actual = inst_copy[key].to_dense() + expected = self.instrument[key].to_dense() + self.assertArraysEqual(actual, expected) def test_to_string(self): inst_str = str(self.instrument) - # TODO assert correctness + self.assertIsInstance(inst_str, str) + self.assertIn("Instrument with elements:", inst_str) + for key in self.instrument.keys(): + self.assertIn(str(key), inst_str) def test_transform(self): T = FullGaugeGroupElement( @@ -30,54 +56,39 @@ def test_transform(self): [0, 1, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0]], 'd')) + T_mx = T.transform_matrix + T_inv = T.transform_matrix_inverse + originals = {k: v.to_dense().copy() for k, v in self.instrument.items()} self.instrument.transform_inplace(T) - # TODO assert correctness + for key, E_orig in originals.items(): + expected = T_mx @ E_orig @ T_inv + self.assertArraysAlmostEqual(self.instrument[key].to_dense(), expected) def test_simplify_operations(self): gates = self.instrument.simplify_operations(prefix="ABC") - # TODO assert correctness - - def test_constructor_raises_on_non_none_param_conflict(self): - with self.assertRaises(AssertionError): - self.constructor(["Non-none-matrices"], 'default', None, False, ["Non-none-items"]) # can't both be non-None + self.assertEqual(len(gates), len(self.instrument)) + expected_keys = ["ABC_" + k for k in self.instrument.keys()] + self.assertEqual(list(gates.keys()), expected_keys) + for gate, orig in zip(gates.values(), self.instrument.values()): + self.assertIs(gate, orig) - def test_constructor_raises_on_bad_op_matrices_type(self): + def test_convert_raises_on_unknown_parameterization(self): with self.assertRaises(ValueError): - self.constructor("foobar") # op_matrices must be a list or dict - - def test_convert_raises_on_unknown_basis(self): - with self.assertRaises(ValueError): - inst.convert(self.instrument, "foobar", self.model.basis) - - -class InstrumentInstanceBase(object): - def setUp(self): - # Initialize standard target model for instruments - # XXX can instruments be tested independently of a model? EGN: yes, I was just lazy; but they should also be tested within a model. - self.n_elements = 32 - - self.model = std.target_model() - E = self.model.povms['Mdefault']['0'] - Erem = self.model.povms['Mdefault']['1'] - self.Gmz_plus = np.dot(E, E.T) - self.Gmz_minus = np.dot(Erem, Erem.T) - # XXX is this used? - self.povm_ident = self.model.povms['Mdefault']['0'] + self.model.povms['Mdefault']['1'] - self.instrument = self.constructor({'plus': self.Gmz_plus, 'minus': self.Gmz_minus}) - self.model.instruments['Iz'] = self.instrument - super(InstrumentInstanceBase, self).setUp() + inst.convert(self.instrument, "H+S", self.model.basis) -class InstrumentInstanceTester(InstrumentMethodBase, InstrumentInstanceBase, BaseCase): +class InstrumentInstanceTester(InstrumentTestBase): + __test__ = True constructor = inst.Instrument -class TPInstrumentInstanceTester(InstrumentMethodBase, InstrumentInstanceBase, BaseCase): +class TPInstrumentInstanceTester(InstrumentTestBase): + __test__ = True constructor = inst.TPInstrument def test_raise_on_modify(self): with self.assertRaises(ValueError): - self.instrument['plus'] = None # can't set value of a TP Instrument element + self.instrument['plus'] = None class CPTPInstrumentTester(BaseCase):