diff --git a/jupyter_notebooks/Examples/Leakage-automagic.ipynb b/jupyter_notebooks/Examples/Leakage-automagic.ipynb new file mode 100644 index 000000000..c586ff789 --- /dev/null +++ b/jupyter_notebooks/Examples/Leakage-automagic.ipynb @@ -0,0 +1,87 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from pygsti.modelpacks import smq1Q_XY, smq1Q_ZN\n", + "from pygsti.tools.leakage import leaky_qubit_model_from_pspec, construct_leakage_report\n", + "from pygsti.data import simulate_data\n", + "from pygsti.protocols import StandardGST, ProtocolData" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## GST: modeling a leaky qubit as a qutrit\n", + "\n", + "This short notebook shows how (data from) an experiment design for a two-level system can be used to fit a three-level sytem model, and how to generate a special report to provide insights for these models. The report includes special gate error metrics that reflect the distinguished role of the first two levels in the three-level system." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "mp = smq1Q_XY\n", + "ed = mp.create_gst_experiment_design(max_max_length=32)\n", + "tm3 = leaky_qubit_model_from_pspec(mp.processor_spec(), mx_basis='l2p1')\n", + "# ^ We could use basis = 'gm' instead of 'l2p1'. We prefer 'l2p1'\n", + "# because it makes process matrices easier to interpret in leakage\n", + "# modeling.\n", + "ds = simulate_data(tm3, ed.all_circuits_needing_data, num_samples=1000, seed=1997)\n", + "gst = StandardGST( modes=('CPTPLND',), target_model=tm3, verbosity=2)\n", + "pd = ProtocolData(ed, ds)\n", + "res = gst.run(pd)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "report_dir = 'example_files/leakage-report-automagic'\n", + "report_object, updated_res = construct_leakage_report(res, title='easy leakage analysis!')\n", + "# ^ Each estimate in updated_res has a new gauge-optimized model.\n", + "# The gauge optimization was done to reflect how our target gates\n", + "# are only _really_ defined on the first two levels of our\n", + "# three-level system.\n", + "# \n", + "report_object.write_html(report_dir)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "leak311", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.13" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/jupyter_notebooks/Examples/Leakage.ipynb b/jupyter_notebooks/Examples/Leakage-manual.ipynb similarity index 100% rename from jupyter_notebooks/Examples/Leakage.ipynb rename to jupyter_notebooks/Examples/Leakage-manual.ipynb diff --git a/pygsti/algorithms/gaugeopt.py b/pygsti/algorithms/gaugeopt.py index 941e72845..c332553fa 100644 --- a/pygsti/algorithms/gaugeopt.py +++ b/pygsti/algorithms/gaugeopt.py @@ -20,7 +20,12 @@ from pygsti import tools as _tools from pygsti.tools import mpitools as _mpit from pygsti.tools import slicetools as _slct -from pygsti.models.gaugegroup import TrivialGaugeGroupElement as _TrivialGaugeGroupElement +from pygsti.models.gaugegroup import ( + TrivialGaugeGroupElement as _TrivialGaugeGroupElement, + GaugeGroupElement as _GaugeGroupElement +) + +from typing import Callable, Union, Optional def gaugeopt_to_target(model, target_model, item_weights=None, @@ -29,7 +34,7 @@ def gaugeopt_to_target(model, target_model, item_weights=None, gauge_group=None, method='auto', maxiter=100000, maxfev=None, tol=1e-8, oob_check_interval=0, convert_model_to=None, return_all=False, comm=None, - verbosity=0, check_jac=False): + verbosity=0, check_jac=False, n_leak=0): """ Optimize the gauge degrees of freedom of a model to that of a target. @@ -170,7 +175,7 @@ def gaugeopt_to_target(model, target_model, item_weights=None, objective_fn, jacobian_fn = _create_objective_fn( model, target_model, item_weights, cptp_penalty_factor, spam_penalty_factor, - gates_metric, spam_metric, method, comm, check_jac) + gates_metric, spam_metric, method, comm, check_jac, n_leak) result = gaugeopt_custom(model, objective_fn, gauge_group, method, maxiter, maxfev, tol, oob_check_interval, @@ -307,9 +312,6 @@ def _call_jacobian_fn(gauge_group_el_vec): printer.log("--- Gauge Optimization (%s method, %s) ---" % (method, str(type(gauge_group))), 2) if method == 'ls': - #minSol = _opt.least_squares(_call_objective_fn, x0, #jac=_call_jacobian_fn, - # max_nfev=maxfev, ftol=tol) - #solnX = minSol.x assert(_call_jacobian_fn is not None), "Cannot use 'ls' method unless jacobian is available" ralloc = _baseobjs.ResourceAllocation(comm) # FUTURE: plumb up a resource alloc object? test_f = _call_objective_fn(x0) @@ -354,10 +356,15 @@ def _call_jacobian_fn(gauge_group_el_vec): return newModel -def _create_objective_fn(model, target_model, item_weights=None, - cptp_penalty_factor=0, spam_penalty_factor=0, +GGElObjective = Callable[[_GaugeGroupElement,bool], Union[float, _np.ndarray]] + +GGElJacobian = Union[None, Callable[[_GaugeGroupElement], _np.ndarray]] + + +def _create_objective_fn(model, target_model, item_weights: Optional[dict[str,float]]=None, + cptp_penalty_factor: float=0.0, spam_penalty_factor: float=0.0, gates_metric="frobenius", spam_metric="frobenius", - method=None, comm=None, check_jac=False): + method=None, comm=None, check_jac=False, n_leak=0) -> tuple[GGElObjective, GGElJacobian]: """ Creates the objective function and jacobian (if available) for gaugeopt_to_target @@ -595,17 +602,32 @@ def _mock_objective_fn(v): # non-least-squares case where objective function returns a single float # and (currently) there's no analytic jacobian + assert gates_metric != "frobeniustt" + assert spam_metric != "frobeniustt" + # ^ PR #410 removed support for Frobenius transform-target metrics in this codepath. + + dim = int(_np.sqrt(mxBasis.dim)) + if n_leak > 0: + B = _tools.leading_dxd_submatrix_basis_vectors(dim - n_leak, dim, mxBasis) + P = B @ B.T.conj() + if _np.linalg.norm(P.imag) > 1e-12: + msg = f"Attempting to run leakage-aware gauge optimization with basis {mxBasis}\n" + msg += "is resulting an orthogonal projector onto the computational subspace that\n" + msg += "is not real-valued. Try again with a different basis, like 'l2p1' or 'gm'." + raise ValueError(msg) + else: + P = P.real + transform_mx_arg = (P, _tools.matrixtools.IdentityOperator()) + # ^ The semantics of this tuple are defined by the frobeniusdist function + # in the ExplicitOpModelCalc class. + else: + transform_mx_arg = None + # ^ It would be equivalent to set this to a pair of IdentityOperator objects. + def _objective_fn(gauge_group_el, oob_check): mdl = _transform_with_oob_check(model, gauge_group_el, oob_check) ret = 0 - if gates_metric == "frobeniustt" or spam_metric == "frobeniustt": - full_target_model = target_model.copy() - full_target_model.convert_members_inplace("full") # so we can gauge-transform the target model. - transformed_target = _transform_with_oob_check(full_target_model, gauge_group_el.inverse(), oob_check) - else: - transformed_target = None - if cptp_penalty_factor > 0: mdl.basis = mxBasis # set basis for jamiolkowski iso cpPenaltyVec = _cptp_penalty(mdl, cptp_penalty_factor, mdl.basis) @@ -616,84 +638,86 @@ def _objective_fn(gauge_group_el, oob_check): spamPenaltyVec = _spam_penalty(mdl, spam_penalty_factor, mdl.basis) ret += _np.sum(spamPenaltyVec) - if target_model is not None: - if gates_metric == "frobenius": - if spam_metric == "frobenius": - ret += mdl.frobeniusdist(target_model, None, item_weights) - else: - wts = item_weights.copy(); wts['spam'] = 0.0 - for k in wts: - if k in mdl.preps or \ - k in mdl.povms: wts[k] = 0.0 - ret += mdl.frobeniusdist(target_model, None, wts) - - elif gates_metric == "frobeniustt": - if spam_metric == "frobeniustt": - ret += transformed_target.frobeniusdist(model, None, item_weights) - else: - wts = item_weights.copy(); wts['spam'] = 0.0 - for k in wts: - if k in mdl.preps or \ - k in mdl.povms: wts[k] = 0.0 - ret += transformed_target.frobeniusdist(model, None, wts) - - elif gates_metric == "fidelity": - for opLbl in mdl.operations: - wt = item_weights.get(opLbl, opWeight) - ret += wt * (1.0 - _tools.entanglement_fidelity( - target_model.operations[opLbl], mdl.operations[opLbl]))**2 - - elif gates_metric == "tracedist": - for opLbl in mdl.operations: - wt = item_weights.get(opLbl, opWeight) - ret += opWeight * _tools.jtracedist( - target_model.operations[opLbl], mdl.operations[opLbl]) - - else: raise ValueError("Invalid gates_metric: %s" % gates_metric) - - if spam_metric == "frobenius": - if gates_metric != "frobenius": # otherwise handled above to match normalization in frobeniusdist - wts = item_weights.copy(); wts['gates'] = 0.0 - for k in wts: - if k in mdl.operations or \ - k in mdl.instruments: wts[k] = 0.0 - ret += mdl.frobeniusdist(target_model, None, wts) - - elif spam_metric == "frobeniustt": - if gates_metric != "frobeniustt": # otherwise handled above to match normalization in frobeniusdist - wts = item_weights.copy(); wts['gates'] = 0.0 - for k in wts: - if k in mdl.operations or \ - k in mdl.instruments: wts[k] = 0.0 - ret += transformed_target.frobeniusdist(model, None, wts) - - elif spam_metric == "fidelity": - for preplabel, prep in mdl.preps.items(): - wt = item_weights.get(preplabel, spamWeight) - rhoMx1 = _tools.vec_to_stdmx(prep, mxBasis) - rhoMx2 = _tools.vec_to_stdmx( - target_model.preps[preplabel], mxBasis) - ret += wt * (1.0 - _tools.fidelity(rhoMx1, rhoMx2))**2 - - for povmlabel, povm in mdl.povms.items(): - wt = item_weights.get(povmlabel, spamWeight) - ret += wt * (1.0 - _tools.povm_fidelity( - mdl, target_model, povmlabel))**2 - - elif spam_metric == "tracedist": - for preplabel, prep in mdl.preps.items(): - wt = item_weights.get(preplabel, spamWeight) - rhoMx1 = _tools.vec_to_stdmx(prep, mxBasis) - rhoMx2 = _tools.vec_to_stdmx( - target_model.preps[preplabel], mxBasis) - ret += wt * _tools.tracedist(rhoMx1, rhoMx2) - - for povmlabel, povm in mdl.povms.items(): - wt = item_weights.get(povmlabel, spamWeight) - ret += wt * (1.0 - _tools.povm_jtracedist( - mdl, target_model, povmlabel))**2 - - else: raise ValueError("Invalid spam_metric: %s" % spam_metric) + if target_model is None: + return ret + + if "frobenius" in gates_metric: + if spam_metric == gates_metric: + val = mdl.frobeniusdist(target_model, transform_mx_arg, item_weights) + else: + wts = item_weights.copy() + wts['spam'] = 0.0 + for k in wts: + if k in mdl.preps or k in mdl.povms: + wts[k] = 0.0 + val = mdl.frobeniusdist(target_model, transform_mx_arg, wts, n_leak) + if "squared" in gates_metric: + val = val ** 2 + ret += val + + elif gates_metric == "fidelity": + # If n_leak==0, then subspace_entanglement_fidelity is just entanglement_fidelity + for opLbl in mdl.operations: + wt = item_weights.get(opLbl, opWeight) + top = target_model.operations[opLbl].to_dense() + mop = mdl.operations[opLbl].to_dense() + ret += wt * (1.0 - _tools.subspace_entanglement_fidelity(top, mop, mxBasis, n_leak))**2 + + elif gates_metric == "tracedist": + # If n_leak==0, then subspace_jtracedist is just jtracedist. + for opLbl in mdl.operations: + wt = item_weights.get(opLbl, opWeight) + top = target_model.operations[opLbl].to_dense() + mop = mdl.operations[opLbl].to_dense() + ret += wt * _tools.subspace_jtracedist(top, mop, mxBasis, n_leak) + + else: + raise ValueError("Invalid gates_metric: %s" % gates_metric) + + if "frobenius" in spam_metric and gates_metric == spam_metric: + # We already handled SPAM error in this case. Just return. + return ret + + if "frobenius" in spam_metric: + # SPAM and gates can have different choices for squared vs non-squared. + wts = item_weights.copy(); wts['gates'] = 0.0 + for k in wts: + if k in mdl.operations or k in mdl.instruments: + wts[k] = 0.0 + val = mdl.frobeniusdist(target_model, transform_mx_arg, wts) + if "squared" in spam_metric: + val = val ** 2 + ret += val + + elif spam_metric == "fidelity": + # Leakage-aware metrics NOT available + for preplabel, m_prep in mdl.preps.items(): + wt = item_weights.get(preplabel, spamWeight) + rhoMx1 = _tools.vec_to_stdmx(m_prep.to_dense(), mxBasis) + t_prep = target_model.preps[preplabel] + rhoMx2 = _tools.vec_to_stdmx(t_prep.to_dense(), mxBasis) + ret += wt * (1.0 - _tools.fidelity(rhoMx1, rhoMx2))**2 + + for povmlabel in mdl.povms.keys(): + wt = item_weights.get(povmlabel, spamWeight) + fidelity = _tools.povm_fidelity(mdl, target_model, povmlabel) + ret += wt * (1.0 - fidelity)**2 + + elif spam_metric == "tracedist": + # Leakage-aware metrics NOT available. + for preplabel, m_prep in mdl.preps.items(): + wt = item_weights.get(preplabel, spamWeight) + rhoMx1 = _tools.vec_to_stdmx(m_prep.to_dense(), mxBasis) + t_prep = target_model.preps[preplabel] + rhoMx2 = _tools.vec_to_stdmx(t_prep.to_dense(), mxBasis) + ret += wt * _tools.tracedist(rhoMx1, rhoMx2) + + for povmlabel in mdl.povms.keys(): + wt = item_weights.get(povmlabel, spamWeight) + ret += wt * _tools.povm_jtracedist(mdl, target_model, povmlabel) + + else: + raise ValueError("Invalid spam_metric: %s" % spam_metric) return ret diff --git a/pygsti/baseobjs/basis.py b/pygsti/baseobjs/basis.py index d84f203eb..9e5734a1e 100644 --- a/pygsti/baseobjs/basis.py +++ b/pygsti/baseobjs/basis.py @@ -794,6 +794,22 @@ def ellookup(self): self._lazy_build_labels() self._ellookup = {lbl: el for lbl, el in zip(self._labels, self._elements)} return self._ellookup + + @property + def elindlookup(self) -> dict: + """ + A dictionary mapping labels to the index in self.elements that holds + that label's basis element. + + Returns + ------- + dict + """ + if self._ellookup is None: + if self._labels is None: + self._lazy_build_labels() + self._elindlookup = {lbl: ind for ind, lbl in enumerate(self._labels)} + return self._elindlookup @property def elements(self): diff --git a/pygsti/baseobjs/basisconstructors.py b/pygsti/baseobjs/basisconstructors.py index 5edd81f47..5ba4e40bd 100644 --- a/pygsti/baseobjs/basisconstructors.py +++ b/pygsti/baseobjs/basisconstructors.py @@ -697,6 +697,52 @@ def gm_labels(matrix_dim): lblList.extend(["Z_{%d}" % (k) for k in range(1, d)]) return lblList +def lf_labels(matrix_dim: int) -> tuple[str,...]: + if matrix_dim != 3: + raise NotImplementedError() + # Divides the underlying Hilbert space into levels (0, 1, L), where (0, 1) + # is the computational subspace and (L,) is leakage space. + lbls = ( + # The first element is proportional to the identity on the computational subspace. + "I", + # The next seven elements also appear in the Gell-Mann basis when matrix_dim == 3. + "X", # --> X_{0,1} + "Y", # --> Y_{0,1} + "Z", # --> Z_{1} + "LX0", # --> X_{0,2} + "LX1", # --> X_{1,2} + "LY0", # --> Y_{0,2} + "LY1", # --> Y_{1,2} + # The ninth and final element is proportional to the identity on leakage space. + "L" + ) + return lbls + +def lf_matrices(matrix_dim: int) -> list[_np.ndarray]: + """ + This basis is used to isolate the parts of Hilbert-Schmidt space that act on + the computational subspace induced from a partition of 3-dimensional complex + Hilbert space into a 2-dimensional computational subspace and a 1-dimensional + leakage space. + """ + if matrix_dim != 3: + raise NotImplementedError() + + gm_basis = gm_matrices(3) + leakage_basis_mxs = [ + _np.sqrt(2) / 3 * (_np.sqrt(3) * gm_basis[0] + 0.5 * _np.sqrt(6) * gm_basis[8]), + gm_basis[1], + gm_basis[4], + gm_basis[7], + gm_basis[2], + gm_basis[3], + gm_basis[5], + gm_basis[6], + 1 / 3 * (_np.sqrt(3) * gm_basis[0] - _np.sqrt(6) * gm_basis[8]), + ] + return leakage_basis_mxs + + def qsim_matrices(matrix_dim): """ @@ -1269,6 +1315,7 @@ def unknown_labels(dim): _basis_constructor_dict['PP'] = MatrixBasisConstructor('Pauli-Product basis', PP_matrices, pp_labels, True, True) _basis_constructor_dict['qsim'] = MatrixBasisConstructor('QuantumSim basis', qsim_matrices, qsim_labels, True, False) _basis_constructor_dict['qt'] = MatrixBasisConstructor('Qutrit basis', qt_matrices, qt_labels, True, True) +_basis_constructor_dict['l2p1'] = MatrixBasisConstructor('Leakage basis for a 2+1 level system', lf_matrices, lf_labels, True, False) _basis_constructor_dict['id'] = SingleElementMatrixBasisConstructor('Identity-only subbasis', identity_matrices, identity_labels, True, True) _basis_constructor_dict['clmx'] = DiagonalMatrixBasisConstructor('Diagonal Matrix-unit basis', cl_vectors, diff --git a/pygsti/circuits/circuit.py b/pygsti/circuits/circuit.py index af57c3d18..78e077021 100644 --- a/pygsti/circuits/circuit.py +++ b/pygsti/circuits/circuit.py @@ -977,11 +977,10 @@ def __eq__(self, x): if isinstance(x, Circuit): if len(self) != len(x): return False + elif self._static and x._static and self._hash != x._hash: + return False else: - if self._static and x._static: - return self._hash == x._hash - else: - return self.tup == x.tup + return self.tup == x.tup elif x is None: return False else: diff --git a/pygsti/data/dataset.py b/pygsti/data/dataset.py index 6b109259e..135062d8c 100644 --- a/pygsti/data/dataset.py +++ b/pygsti/data/dataset.py @@ -26,6 +26,7 @@ from pygsti.circuits import circuit as _cir from pygsti.baseobjs import outcomelabeldict as _ld, _compatibility as _compat from pygsti.baseobjs.mongoserializable import MongoSerializable as _MongoSerializable +from pygsti.baseobjs.label import Label as _Label from pygsti.tools import NamedDict as _NamedDict from pygsti.tools import listtools as _lt from pygsti.tools.legacytools import deprecate as _deprecated_fn @@ -1176,7 +1177,6 @@ def _get_row(self, circuit): # needed because name-only Labels don't hash the same as strings # so key lookups need to be done at least with tuples of Labels. circuit = _cir.Circuit.cast(circuit) - #Note: cirIndex value is either an int (non-static) or a slice (static) cirIndex = self.cirIndex[circuit] repData = self.repData[cirIndex] if (self.repData is not None) else None @@ -1292,11 +1292,28 @@ def gate_labels(self, prefix='G'): list of strings A list where each element is a operation label. """ - opLabels = [] - for opLabelString in self: - for opLabel in opLabelString: - if not prefix or opLabel.name.startswith(prefix): - if opLabel not in opLabels: opLabels.append(opLabel) + opLabels = set() + empty_label = _Label(()) + for circuit in self: + for layer in circuit.layertup: + if isinstance(layer, list): + for lbl in layer: + if lbl.name.startswith(prefix): + opLabels.add(lbl.name) + elif lbl == empty_label: + opLabels.add('[]') + elif isinstance(layer, _Label): + if layer.name.startswith(prefix): + opLabels.add(layer.name) + else: + raise ValueError() + + bracket_idle = '[]' in opLabels + opLabels = [op for op in opLabels if op != '[]'] + opLabels.sort() + if bracket_idle: + opLabels = ['[]'] + opLabels + return opLabels def degrees_of_freedom(self, circuits=None, method="present_outcomes-1", diff --git a/pygsti/evotypes/basereps.py b/pygsti/evotypes/basereps.py index 01eb8246d..449c0e2b5 100644 --- a/pygsti/evotypes/basereps.py +++ b/pygsti/evotypes/basereps.py @@ -78,7 +78,7 @@ def __init__(self, int_coeff_dict, max_num_vars, vindices_per_int): def reinit(self, int_coeff_dict): """ - Reinitialize this polynomial using new coefficents. + Reinitialize this polynomial using new coefficients. Parameters ---------- @@ -131,7 +131,7 @@ def copy(self): def abs(self): """ - Return a polynomial whose coefficents are the absolute values of this PolynomialRep's coefficients. + Return a polynomial whose coefficients are the absolute values of this PolynomialRep's coefficients. Returns ------- diff --git a/pygsti/extras/idletomography/idtcore.py b/pygsti/extras/idletomography/idtcore.py index 71dc3edf5..2ff10de74 100644 --- a/pygsti/extras/idletomography/idtcore.py +++ b/pygsti/extras/idletomography/idtcore.py @@ -601,7 +601,8 @@ def extract_action(g, cur_sslbls, ql): # StaticArbitraryOp, LindbladDenseOp, other gates... if len(cur_sslbls) == 1 and cur_sslbls[0] == ql: mx = g.to_dense() - assert(mx.shape == (4, 4)) + if mx.shape != (4,4): + raise ValueError(f'Expected matrix of shape (4,4), got {mx.shape}.') return mx else: mx = g.to_dense() diff --git a/pygsti/modelmembers/modelmembergraph.py b/pygsti/modelmembers/modelmembergraph.py index 014bd9a98..5b2b2e743 100644 --- a/pygsti/modelmembers/modelmembergraph.py +++ b/pygsti/modelmembers/modelmembergraph.py @@ -19,8 +19,8 @@ class ModelMemberGraph(object): """A directed acyclic graph of dependencies of ModelMembers""" - @classmethod - def load_modelmembers_from_serialization_dict(cls, sdict, parent_model): + @staticmethod + def load_modelmembers_from_serialization_dict(sdict, parent_model): """Create a nested dictionary of model members from a previously serialized graph. Parameters diff --git a/pygsti/modelmembers/operations/linearop.py b/pygsti/modelmembers/operations/linearop.py index a2ab0c0ea..7a323145c 100644 --- a/pygsti/modelmembers/operations/linearop.py +++ b/pygsti/modelmembers/operations/linearop.py @@ -17,6 +17,8 @@ from pygsti.tools import optools as _ot from pygsti import SpaceT +from typing import Any + #Note on initialization sequence of Operations within a Model: # 1) a Model is constructed (empty) # 2) a LinearOperator is constructed - apart from a Model if it's locally parameterized, @@ -27,7 +29,7 @@ # 3) the LinearOperator is assigned/added to a dict within the Model. As a part of this # process, the LinearOperator's 'gpindices' member is set, if it isn't already, and the # Model's "global" parameter vector (and number of params) is updated as -# needed to accomodate new parameters. +# needed to accommodate new parameters. # # Note: gpindices may be None (before initialization) or any valid index # into a 1D numpy array (e.g. a slice or integer array). It may NOT have @@ -44,7 +46,7 @@ class LinearOperator(_modelmember.ModelMember): """ - Base class for all operation representations + Base class for all *square* operation representations Parameters ---------- @@ -92,6 +94,14 @@ def hilbert_schmidt_size(self): """ return (self.dim)**2 + @property + def shape(self) -> tuple[int, int]: + # Provide this function to mimic numpy array semantics. + # + # We can't rely on self._rep.shape since superclasses + # are given broad freedom to define semantics of self._rep. + return (self.dim, self.dim) + def set_dense(self, m): """ Set the dense-matrix value of this operation. @@ -381,13 +391,13 @@ def taylor_order_terms_above_mag(self, order, max_polynomial_vars, min_term_mag) return [t for t in terms_at_order if t.magnitude >= min_term_mag] - def frobeniusdist_squared(self, other_op, transform=None, inv_transform=None): + def frobeniusdist_squared(self, other_op, transform=None, inv_transform=None) -> _np.floating[Any]: """ Return the squared frobenius difference between this operation and `other_op` Optionally transforms this operation first using matrices `transform` and `inv_transform`. Specifically, this operation gets - transfomed as: `O => inv_transform * O * transform` before comparison with + transformed as: `O => inv_transform * O * transform` before comparison with `other_op`. Parameters @@ -405,12 +415,13 @@ def frobeniusdist_squared(self, other_op, transform=None, inv_transform=None): ------- float """ - if transform is None and inv_transform is None: - return _ot.frobeniusdist_squared(self.to_dense("minimal"), other_op.to_dense("minimal")) - else: - return _ot.frobeniusdist_squared(_np.dot( - inv_transform, _np.dot(self.to_dense("minimal"), transform)), - other_op.to_dense("minimal")) + self_mx = self.to_dense("minimal") + if transform is not None: + self_mx = self_mx @ transform + if inv_transform is not None: + self_mx = inv_transform @ self_mx + return _ot.frobeniusdist_squared(self_mx, other_op.to_dense("minimal")) + def frobeniusdist(self, other_op, transform=None, inv_transform=None): """ @@ -418,7 +429,7 @@ def frobeniusdist(self, other_op, transform=None, inv_transform=None): Optionally transforms this operation first using matrices `transform` and `inv_transform`. Specifically, this operation gets - transfomed as: `O => inv_transform * O * transform` before comparison with + transformed as: `O => inv_transform * O * transform` before comparison with `other_op`. Parameters diff --git a/pygsti/modelmembers/povms/__init__.py b/pygsti/modelmembers/povms/__init__.py index 332102871..f05a636e2 100644 --- a/pygsti/modelmembers/povms/__init__.py +++ b/pygsti/modelmembers/povms/__init__.py @@ -547,10 +547,14 @@ def convert(povm, to_type, basis, ideal_povm=None, flatten_structure=False, cp_p #It is often the case that there are more error generators than physical degrees of freedom in the POVM #We define a function which finds linear comb. of errgens that span these degrees of freedom. #This has been called "the trivial gauge", and this function is meant to avoid it + eegb = 'PP' if basis.name == 'pp' else basis def calc_physical_subspace(dense_ideal_povm, epsilon = 1e-4): degrees_of_freedom = (dense_ideal_povm.shape[0] - 1) * dense_ideal_povm.shape[1] - errgen = _LindbladErrorgen.from_error_generator(povm.state_space.dim, parameterization=to_type) + errgen = _LindbladErrorgen.from_error_generator( + povm.state_space.dim, parameterization=to_type, + elementary_errorgen_basis=eegb, mx_basis=basis + ) if degrees_of_freedom > errgen.num_params: warnings.warn("POVM has more degrees of freedom than the available number of parameters, representation in this parameterization is not guaranteed") diff --git a/pygsti/modelmembers/povms/composedpovm.py b/pygsti/modelmembers/povms/composedpovm.py index 345ccbfbe..9f2f93bcd 100644 --- a/pygsti/modelmembers/povms/composedpovm.py +++ b/pygsti/modelmembers/povms/composedpovm.py @@ -338,7 +338,7 @@ def transform_inplace(self, s): """ self.error_map.spam_transform_inplace(s, 'effect') # only do this *once* for lbl, effect in self.items(): - #effect._update_rep() # these two lines mimic the bookeeping in + #effect._update_rep() # these two lines mimic the bookkeeping in effect.dirty = True # a "effect.transform_inplace(s, 'effect')" call. self.dirty = True diff --git a/pygsti/modelmembers/povms/effect.py b/pygsti/modelmembers/povms/effect.py index 72920f973..e1e0d46ee 100644 --- a/pygsti/modelmembers/povms/effect.py +++ b/pygsti/modelmembers/povms/effect.py @@ -16,6 +16,8 @@ from pygsti.tools import optools as _ot from pygsti.baseobjs.opcalc import bulk_eval_compact_polynomials_complex as _bulk_eval_compact_polynomials_complex +from typing import Any + class POVMEffect(_modelmember.ModelMember): """ @@ -117,13 +119,11 @@ def set_time(self, t): ## PUT term calc methods here if appropriate... - def frobeniusdist_squared(self, other_spam_vec, transform=None, - inv_transform=None): + def frobeniusdist_squared(self, other_spam_vec, transform=None, inv_transform=None) -> _np.floating[Any]: """ - Return the squared frobenius difference between this operation and `other_spam_vec`. + Return the squared frobenius difference between this effect and `other_spam_vec`. - Optionally transforms this vector first using `transform` and - `inv_transform`. + Optionally transforms this vector first using `transform`. Parameters ---------- @@ -134,18 +134,17 @@ def frobeniusdist_squared(self, other_spam_vec, transform=None, Transformation matrix. inv_transform : numpy.ndarray, optional - Inverse of `tranform`. + Ignored. (We keep this as a positional argument for consistency with + the frobeniusdist_squared method of pyGSTi's LinearOperator objects.) Returns ------- float """ vec = self.to_dense() - if transform is None: - return _ot.frobeniusdist_squared(vec, other_spam_vec.to_dense()) - else: - return _ot.frobeniusdist_squared(_np.dot(_np.transpose(transform), - vec), other_spam_vec.to_dense()) + if transform is not None: + vec = transform.T @ vec + return _ot.frobeniusdist_squared(vec, other_spam_vec.to_dense()) def residuals(self, other_spam_vec, transform=None, inv_transform=None): """ diff --git a/pygsti/modelmembers/povms/marginalizedpovm.py b/pygsti/modelmembers/povms/marginalizedpovm.py index e52c1e159..b46775d75 100644 --- a/pygsti/modelmembers/povms/marginalizedpovm.py +++ b/pygsti/modelmembers/povms/marginalizedpovm.py @@ -184,7 +184,7 @@ def __getitem__(self, key): if key in self: # calls __contains__ to efficiently check for membership #create effect vector now that it's been requested (lazy creation) - #FUTURE: maybe have a "SumPOVMEffect" that can add spamvecs to preserve paramterization and avoid dense reps + #FUTURE: maybe have a "SumPOVMEffect" that can add spamvecs to preserve parameterization and avoid dense reps effect_vec = None # Note: currently all marginalized POVMs are *static*, since # we don't have a good general way to add parameterized effect vectors. diff --git a/pygsti/modelmembers/states/densestate.py b/pygsti/modelmembers/states/densestate.py index 714b39913..53e12bcbf 100644 --- a/pygsti/modelmembers/states/densestate.py +++ b/pygsti/modelmembers/states/densestate.py @@ -243,9 +243,33 @@ def to_memoized_dict(self, mmg_memo): @classmethod def _from_memoized_dict(cls, mm_dict, serial_memo): vec = cls._decodemx(mm_dict['dense_superket_vector']) - state_space = _statespace.StateSpace.from_nice_serialization(mm_dict['state_space']) - basis = _Basis.from_nice_serialization(mm_dict['basis']) if (mm_dict['basis'] is not None) else None - return cls(vec, basis, mm_dict['evotype'], state_space) + try: + state_space = _statespace.StateSpace.from_nice_serialization(mm_dict['state_space']) + basis = _Basis.from_nice_serialization(mm_dict['basis']) if (mm_dict['basis'] is not None) else None + return cls(vec, basis, mm_dict['evotype'], state_space) + except AssertionError as e: + """ + This codepath can get hit when deserializing TPPOVM or UnconstrainedPOVM objects. + + More specifically, it can get hit when objects other than POVMEffect vectors were + passed to the constructors of these classes. When that happens, their base class + constructor (for BasePOVM) constructs the effects as FullPOVMEffect objects (which in + turn rely on FullState and then DenseState) with None passed for the basis argument. + Somewhere downstream that None gets cast to a Basis object where `.dim == 0`, which + is inconsistent with the array representation of the effect having length > 0. + + In an ideal world we'd have POVMs enforce not only non-None bases but also *consistent* + bases for all constituent effects. For now, our fix is to just try and recover the basis + from serial_memo. (Perhaps not-coincidentally, this function wasn't using the serial_memo + argument before this code path was added.) + """ + se = str(e) + if 'Basis object has unexpected dimension' in se and len(serial_memo) > 0: + member = list(serial_memo.values())[0] + basis = member.parent.basis + state_space = basis.state_space + return cls(vec, basis, mm_dict['evotype'], state_space) + raise e def _is_similar(self, other, rtol, atol): """ diff --git a/pygsti/modelmembers/states/state.py b/pygsti/modelmembers/states/state.py index da693a957..6cea48cae 100644 --- a/pygsti/modelmembers/states/state.py +++ b/pygsti/modelmembers/states/state.py @@ -19,6 +19,8 @@ from pygsti.tools import optools as _ot from pygsti import SpaceT +from typing import Any + class State(_modelmember.ModelMember): """ @@ -298,13 +300,11 @@ def taylor_order_terms_above_mag(self, order, max_polynomial_vars, min_term_mag) terms_at_order = [t.copy_with_magnitude(abs(coeff)) for coeff, t in zip(coeffs, terms_at_order)] return [t for t in terms_at_order if t.magnitude >= min_term_mag] - def frobeniusdist_squared(self, other_spam_vec, transform=None, - inv_transform=None): + def frobeniusdist_squared(self, other_spam_vec, transform=None, inv_transform=None) -> _np.floating[Any]: """ Return the squared frobenius difference between this operation and `other_spam_vec`. - Optionally transforms this vector first using `transform` and - `inv_transform`. + Optionally transforms this vector first using `inv_transform`. Parameters ---------- @@ -312,21 +312,21 @@ def frobeniusdist_squared(self, other_spam_vec, transform=None, The other spam vector transform : numpy.ndarray, optional - Transformation matrix. + Ignored. (We keep this as a positional argument for consistency with + the frobeniusdist_squared method of pyGSTi's LinearOperator objects.) inv_transform : numpy.ndarray, optional - Inverse of `tranform`. + Named for consistency with the frobeniusdist_squared method of pyGSTi's + LinearOperator objects. Returns ------- float """ vec = self.to_dense("minimal") - if inv_transform is None: - return _ot.frobeniusdist_squared(vec, other_spam_vec.to_dense("minimal")) - else: - return _ot.frobeniusdist_squared(_np.dot(inv_transform, vec), - other_spam_vec.to_dense("minimal")) + if inv_transform is not None: + vec = inv_transform @ vec + return _ot.frobeniusdist_squared(vec, other_spam_vec.to_dense("minimal")) def residuals(self, other_spam_vec, transform=None, inv_transform=None): """ diff --git a/pygsti/models/explicitcalc.py b/pygsti/models/explicitcalc.py index 14803141b..6717e1f66 100644 --- a/pygsti/models/explicitcalc.py +++ b/pygsti/models/explicitcalc.py @@ -16,10 +16,15 @@ import warnings as _warnings import numpy as _np +import scipy.linalg as _la from pygsti.baseobjs import basisconstructors as _bc from pygsti.tools import matrixtools as _mt +from typing import Union + +TransformMxPair = tuple[_np.ndarray, Union[_np.ndarray, _mt.IdentityOperator]] + # Tolerace for matrix_rank when finding rank of a *normalized* projection # matrix. This is a legitimate tolerace since a normalized projection matrix # should have just 0 and 1 eigenvalues, and thus a tolerace << 1.0 will work @@ -32,7 +37,7 @@ class ExplicitOpModelCalc(object): Performs calculations with explicitly-represented objects. This class performs calculations with *simplified* objects (so don't - need to worry abount POVMs or Instruments, just preps, ops, & effects), + need to worry about POVMs or Instruments, just preps, ops, & effects), but, unlike forward simulators, these calculations require knowledge of *all* of the possible operations in each category (not just the ones in a given circuti). As such, instances of `ExplicitOpModelCalc` are almost always @@ -41,7 +46,7 @@ class ExplicitOpModelCalc(object): Parameters ---------- dim : int - The dimenstion of the Hilbert-Schmidt space upon which the + The dimension of the Hilbert-Schmidt space upon which the various operators act. simplified_preps : dict @@ -65,7 +70,7 @@ def __init__(self, dim, simplified_preps, simplified_ops, simplified_effects, np Parameters ---------- dim : int - The dimenstion of the Hilbert-Schmidt space upon which the + The dimension of the Hilbert-Schmidt space upon which the various operators act. simplified_preps, simplified_ops, simplified_effects : dict @@ -105,7 +110,7 @@ def copy(self): """ return ExplicitOpModelCalc(self.dim, self.preps, self.operations, self.effects, self.Np, self.interposer) - def frobeniusdist(self, other_calc, transform_mx=None, + def frobeniusdist(self, other_calc, transform_mx: Union[None, _np.ndarray, TransformMxPair]=None, item_weights=None, normalize=True): """ Compute the weighted frobenius norm of the difference between this calc object and `other_calc`. @@ -120,12 +125,17 @@ def frobeniusdist(self, other_calc, transform_mx=None, other_calc : ForwardSimulator the other gate calculator to difference against. - transform_mx : numpy array, optional - if not None, transform this model by - G => inv(transform_mx) * G * transform_mx, for each operation matrix G - (and similar for rho and E vectors) before taking the difference. - This transformation is applied only for the difference and does - not alter the values stored in this model. + transform_mx : numpy array or tuple, optional + If transform_mx is a numpy array, then for each operation matrix + G we implicitly consider the transformed quantity + G => inv(transform_mx) * G * transform_mx + Similar transformations are applied for effect vectors. + This transformation is applied only for the difference and does not + alter the values stored in this model. + + If transform_mx is a tuple, then its first entry should be a numpy ndarray + that will be interpreted as transform_mx in the usual sense, and its second + entry will be syntactically substituted for inv(transform_mx). item_weights : dict, optional Dictionary of weighting factors for individual gates and spam @@ -146,51 +156,37 @@ def frobeniusdist(self, other_calc, transform_mx=None, ------- float """ - d = 0; T = transform_mx + if transform_mx is None: + P, invP = None, None + # ^ It would be equivalent to use P = invP = _mt.IdentityOperator() + elif isinstance(transform_mx, _np.ndarray): + P = transform_mx + invP = _np.linalg.pinv(P) + else: + P, invP = transform_mx + assert _mt.is_operatorlike(P) + assert _mt.is_operatorlike(invP) + + d = 0.0 nSummands = 0.0 if item_weights is None: item_weights = {} opWeight = item_weights.get('gates', 1.0) spamWeight = item_weights.get('spam', 1.0) - if T is not None: - Ti = _np.linalg.inv(T) # TODO: generalize inverse op (call T.inverse() if T were a "transform" object?) - for opLabel, gate in self.operations.items(): - wt = item_weights.get(opLabel, opWeight) - d += wt * gate.frobeniusdist_squared( - other_calc.operations[opLabel], T, Ti) - nSummands += wt * (gate.dim)**2 - - for lbl, rhoV in self.preps.items(): - wt = item_weights.get(lbl, spamWeight) - d += wt * rhoV.frobeniusdist_squared(other_calc.preps[lbl], T, Ti) - nSummands += wt * rhoV.dim - - for lbl, Evec in self.effects.items(): - wt = item_weights.get(lbl, spamWeight) - d += wt * Evec.frobeniusdist_squared(other_calc.effects[lbl], T, Ti) - nSummands += wt * Evec.dim + for opLabel, gate in self.operations.items(): + wt = item_weights.get(opLabel, opWeight) + d += wt * gate.frobeniusdist_squared(other_calc.operations[opLabel], P, invP) + nSummands += wt * (gate.dim)**2 - else: - for opLabel, gate in self.operations.items(): - wt = item_weights.get(opLabel, opWeight) - d += wt * gate.frobeniusdist_squared(other_calc.operations[opLabel]) - nSummands += wt * (gate.dim)**2 - - for lbl, rhoV in self.preps.items(): - wt = item_weights.get(lbl, spamWeight) - d += wt * rhoV.frobeniusdist_squared(other_calc.preps[lbl]) - nSummands += wt * rhoV.dim - - for lbl, Evec in self.effects.items(): - wt = item_weights.get(lbl, spamWeight) - d += wt * Evec.frobeniusdist_squared(other_calc.effects[lbl]) - nSummands += wt * Evec.dim - - #Temporary: check that this function can be computed by - # calling residuals - replace with this later. - resids, chk_nSummands = self.residuals(other_calc, transform_mx, item_weights) - assert(_np.isclose(_np.sum(resids**2), d)) - assert(_np.isclose(chk_nSummands, nSummands)) + for lbl, rhoV in self.preps.items(): + wt = item_weights.get(lbl, spamWeight) + d += wt * rhoV.frobeniusdist_squared(other_calc.preps[lbl], P, invP) + nSummands += wt * rhoV.dim + + for lbl, Evec in self.effects.items(): + wt = item_weights.get(lbl, spamWeight) + d += wt * Evec.frobeniusdist_squared(other_calc.effects[lbl], P, invP) + nSummands += wt * Evec.dim if normalize and nSummands > 0: return _np.sqrt(d / nSummands) @@ -241,47 +237,30 @@ def residuals(self, other_calc, transform_mx=None, item_weights=None): sqrt_itemWeights = {k: _np.sqrt(v) for k, v in item_weights.items()} opWeight = sqrt_itemWeights.get('gates', 1.0) spamWeight = sqrt_itemWeights.get('spam', 1.0) - - if T is not None: - Ti = _np.linalg.inv(T) # TODO: generalize inverse op (call T.inverse() if T were a "transform" object?) - for opLabel, gate in self.operations.items(): - wt = sqrt_itemWeights.get(opLabel, opWeight) - resids.append( - wt * gate.residuals( - other_calc.operations[opLabel], T, Ti)) - nSummands += wt**2 * (gate.dim)**2 - - for lbl, rhoV in self.preps.items(): - wt = sqrt_itemWeights.get(lbl, spamWeight) - resids.append( - wt * rhoV.residuals(other_calc.preps[lbl], T, Ti)) - nSummands += wt**2 * rhoV.dim - - for lbl, Evec in self.effects.items(): - wt = sqrt_itemWeights.get(lbl, spamWeight) - resids.append( - wt * Evec.residuals(other_calc.effects[lbl], T, Ti)) - - nSummands += wt**2 * Evec.dim - - else: - for opLabel, gate in self.operations.items(): - wt = sqrt_itemWeights.get(opLabel, opWeight) - resids.append( - wt * gate.residuals(other_calc.operations[opLabel])) - nSummands += wt**2 * (gate.dim)**2 - - for lbl, rhoV in self.preps.items(): - wt = sqrt_itemWeights.get(lbl, spamWeight) - resids.append( - wt * rhoV.residuals(other_calc.preps[lbl])) - nSummands += wt**2 * rhoV.dim - - for lbl, Evec in self.effects.items(): - wt = sqrt_itemWeights.get(lbl, spamWeight) - resids.append( - wt * Evec.residuals(other_calc.effects[lbl])) - nSummands += wt**2 * Evec.dim + Ti = None if T is None else _np.linalg.pinv(T) + # ^ TODO: generalize inverse op (call T.inverse() if T were a "transform" object?) + + for opLabel, gate in self.operations.items(): + wt = sqrt_itemWeights.get(opLabel, opWeight) + other_gate = other_calc.operations[opLabel] + resid = wt * gate.residuals(other_gate, T, Ti) + resids.append(resid) + nSummands += wt**2 * (gate.dim)**2 + + for lbl, rhoV in self.preps.items(): + wt = sqrt_itemWeights.get(lbl, spamWeight) + other_prep = other_calc.preps[lbl] + resid = wt * rhoV.residuals(other_prep, T, Ti) + resids.append(resid) + nSummands += wt**2 * rhoV.dim + + for lbl, Evec in self.effects.items(): + wt = sqrt_itemWeights.get(lbl, spamWeight) + other_effect = other_calc.effects[lbl] + resid = wt * Evec.residuals(other_effect, T, Ti) + resids.append(resid) + + nSummands += wt**2 * Evec.dim resids = [r.ravel() for r in resids] resids = _np.concatenate(resids) @@ -532,7 +511,7 @@ def _buildup_dpg(self): # #use acton... maybe throw error if dim is too large (maybe above?) # deriv = _np.zeros((dim,dim),'d') # uv = _np.zeros((dim,1),'d') # unit vec - # for k in range(dim): #FUTURE: could optimize this by bookeeping and pulling this loop outward + # for k in range(dim): #FUTURE: could optimize this by bookkeeping and pulling this loop outward # uv[k] = 1.0; Guv = gate.acton(uv); uv[k] = 0.0 #get k-th col of operation matrix # # termA_mn = sum( U_mk*Gkn ) so U locks m=i,k=j => termA_in = 1.0*Gjn # # termB_mn = sum( Gmk*U_kn ) so U locks k=i,n=j => termB_mj = 1.0*Gmi @@ -590,7 +569,7 @@ def nongauge_and_gauge_spaces(self, item_weights=None, non_gauge_mix_mx=None): orthog_to = gauge_space + _np.dot(nonGaugeDirections, non_gauge_mix_mx) # add non-gauge components in # dims: (nParams,n_gauge_params) + (nParams,n_non_gauge_params) * (n_non_gauge_params,n_gauge_params) # non_gauge_mix_mx is a (n_non_gauge_params,n_gauge_params) matrix whose i-th column specifies the - # coefficents to multipy each of the non-gauge directions by before adding them to the i-th + # coefficients to multiply each of the non-gauge directions by before adding them to the i-th # direction to project out (i.e. what were the pure gauge directions). elif item_weights is not None: @@ -609,19 +588,23 @@ def nongauge_and_gauge_spaces(self, item_weights=None, non_gauge_mix_mx=None): else: orthog_to = gauge_space - #OLD: nongauge_space = _mt.nullspace(orthog_to.T) #cols are non-gauge directions - nongauge_space = _mt.nullspace_qr(orthog_to.T) # cols are non-gauge directions - # print("DB: nullspace of gen_dG (shape = %s, rank=%d) = %s" \ - # % (str(gen_dG.shape),_np.linalg.matrix_rank(gen_dG),str(gen_ndG.shape))) + u,s,_ = _la.svd(orthog_to, full_matrices=True, compute_uv=True) + TOL = 1e-7 + numerical_rank = _np.count_nonzero(s >= TOL*s[0]) + gauge_space = u[:, :numerical_rank] + nongauge_space = u[:, numerical_rank:] + + """NOTE: our use of SVD above. - #REMOVE - ## reduce gen_dG if it doesn't have full rank - #u, s, vh = _np.linalg.svd(gen_dG, full_matrices=False) - #rank = _np.count_nonzero(s > P_RANK_TOL) - #if rank < gen_dG.shape[1]: - # gen_dG = u[:, 0:rank] + We used to use nullspace_qr instead of SVD, and we used to have a dimension check + at the end of this function to ensure that the span of gauge and non-gauge space + gave us the full space in which these subspaces lived. + + We switched to SVD because, if item_weights are not None, it's possible that the + `metric` matrix is singular, which causes orthog_to to be rank-deficient (i.e., + neither its rows nor its columns are linearly independent). + """ - assert(nongauge_space.shape[0] == gauge_space.shape[0] == nongauge_space.shape[1] + gauge_space.shape[1]) return nongauge_space, gauge_space #UNUSED - just used for checking understanding of where the nonzero logL Hessian on gauge space comes from. @@ -715,7 +698,7 @@ def nongauge_projector(self, item_weights=None, non_gauge_mix_mx=None): item_weights : dict, optional Dictionary of weighting factors for individual gates and spam operators. Keys can be gate, state preparation, POVM effect, spam labels, or the - special strings "gates" or "spam" whic represent the entire set of gate + special strings "gates" or "spam" which represent the entire set of gate or SPAM operators, respectively. Values are floating point numbers. These weights define the metric used to compute the non-gauge space, *orthogonal* the gauge space, that is projected onto. @@ -752,7 +735,7 @@ def nongauge_projector(self, item_weights=None, non_gauge_mix_mx=None): # An element of the gauge group can be written gg = exp(-K), where K is a n x n matrix. If K is # hermitian then gg is unitary, but this need not be the case. A gauge transform acts on a # gatset via Model => gg^-1 G gg, i.e. G => exp(K) G exp(-K). We care about the action of - # infinitesimal gauge tranformations (b/c the *derivative* vectors span the tangent space), + # infinitesimal gauge transformations (b/c the *derivative* vectors span the tangent space), # which act as: # G => (I+K) G (I-K) = G + [K,G] + ignored(K^2), where [K,G] := KG-GK # @@ -796,18 +779,18 @@ def nongauge_projector(self, item_weights=None, non_gauge_mix_mx=None): # x[0:nOpParams] is the basis vector for the intersection space within the gate parameter space, # that is, the analogue of dParams_ij in the single-dG_ij introduction above. # - # Still, we just substitue these dParams_ij vectors (as many as the nullspace dimension) for dG_ij + # Still, we just substitute these dParams_ij vectors (as many as the nullspace dimension) for dG_ij # above to get the general case projector. gen_ndG, _ = self.nongauge_and_gauge_spaces(item_weights, non_gauge_mix_mx) - # ORIG WAY: use psuedo-inverse to normalize projector. Ran into problems where + # ORIG WAY: use pseudo-inverse to normalize projector. Ran into problems where # default rcond == 1e-15 didn't work for 2-qubit case, but still more stable than inv method below P = _np.dot(gen_ndG, _np.transpose(gen_ndG)) # almost a projector, but cols of dG are not orthonormal Pp = _np.dot(_np.linalg.pinv(P, rcond=1e-7), P) # make P into a true projector (onto gauge space) # ALT WAY: use inverse of dG^T*dG to normalize projector (see wikipedia on projectors, dG => A) # This *should* give the same thing as above, but numerical differences indicate the pinv method - # is prefereable (so long as rcond=1e-7 is ok in general...) + # is preferable (so long as rcond=1e-7 is ok in general...) # Check: P'*P' = (dG (dGT dG)^1 dGT)(dG (dGT dG)^-1 dGT) = (dG (dGT dG)^1 dGT) = P' #invGG = _np.linalg.inv(_np.dot(_np.transpose(gen_ndG), gen_ndG)) #Pp_alt = _np.dot(gen_ndG, _np.dot(invGG, _np.transpose(gen_ndG))) # a true projector (onto gauge space) diff --git a/pygsti/models/explicitmodel.py b/pygsti/models/explicitmodel.py index 52756cc5a..149c70f28 100644 --- a/pygsti/models/explicitmodel.py +++ b/pygsti/models/explicitmodel.py @@ -24,7 +24,10 @@ from pygsti.models.layerrules import LayerRules as _LayerRules from pygsti.models.modelparaminterposer import ModelParamsInterposer as _ModelParamsInterposer from pygsti.models.fogistore import FirstOrderGaugeInvariantStore as _FirstOrderGaugeInvariantStore -from pygsti.models.gaugegroup import GaugeGroup as _GaugeGroup +from pygsti.models.gaugegroup import ( + GaugeGroup as _GaugeGroup, + GaugeGroupElement as _GaugeGroupElement +) from pygsti.forwardsims.forwardsim import ForwardSimulator as _FSim from pygsti.forwardsims import matrixforwardsim as _matrixfwdsim from pygsti.modelmembers import instruments as _instrument @@ -591,7 +594,7 @@ def compute_nongauge_projector(self, item_weights=None, non_gauge_mix_mx=None): item_weights : dict, optional Dictionary of weighting factors for individual gates and spam operators. Keys can be gate, state preparation, POVM effect, spam labels, or the - special strings "gates" or "spam" whic represent the entire set of gate + special strings "gates" or "spam" which represent the entire set of gate or SPAM operators, respectively. Values are floating point numbers. These weights define the metric used to compute the non-gauge space, *orthogonal* the gauge space, that is projected onto. @@ -1690,6 +1693,72 @@ def _op_decomposition(self, op_label): return self.operations[op_label], self.operations[op_label] +def transform_composed_model(mdl: ExplicitOpModel, s : _GaugeGroupElement) -> ExplicitOpModel: + """ + Return a copy of `mdl` whose members have been gauge-transformed by `s`, + while retaining the parameterization of `mdl`. + + This function's implementation requires that `mdl` use ComposedState for + stateprep and ComposedPOVM for measurements. It does NOT require that + operations be represented with ComposedOp. It ignores any factories that + might be present in mdl. + """ + from pygsti.models import ExplicitOpModel + assert isinstance(mdl, ExplicitOpModel) + if len(mdl.factories) > 0: + _warnings.warn('The returned model will not retain the factories in mdl.') + if len(mdl.instruments) > 0: + raise NotImplementedError('Models with instruments are not supported.') + + oldmdl = mdl + + def mycopy(_m): + # This function is a hack. It makes us robust to errors + # arising from copy.deepcopy in (re)linking model members + # to parent model objects. + s = _m.to_nice_serialization() + t = ExplicitOpModel.from_nice_serialization(s) + return t + + mdl = mycopy(oldmdl) + + from pygsti.modelmembers.operations import ComposedOp, StaticArbitraryOp + from pygsti.modelmembers.povms import ComposedPOVM + from pygsti.modelmembers.states import ComposedState + + U = StaticArbitraryOp(s.transform_matrix, basis=oldmdl.basis) + invU = StaticArbitraryOp(s.transform_matrix_inverse, basis=oldmdl.basis) + + # NOTE: the operations passed to ComposedOp are interpreted in + # reverse order. For example, ComposedOp([X,Y,Z]) is applied to + # a vector v as Z @ Y @ X @ v. + + for key, rho in oldmdl.preps.items(): + # replace each ComposedState superket `rhoVec` with `invU @ rhoVec`; + # do this by packing invU into a new ComposedState's error map. + assert isinstance(rho, ComposedState) + static_rho = rho.state_vec + errmap = ComposedOp([rho.error_map, invU]) + mdl.preps[key] = ComposedState(static_rho, errmap) + + for key, povm in oldmdl.povms.items(): + # replace each ComposedPOVM `p` with another ComposedPOVM `q`, where + # effects `Evec` belonging to `p` are mapped to effects `EVec @ U` + # belonging to `q`. Do this by packing U into the error map of `q`. + assert isinstance(povm, ComposedPOVM) + static_povm = povm.base_povm + errmap = ComposedOp([U, povm.error_map]) + mdl.povms[key] = ComposedPOVM(errmap, static_povm, mx_basis=oldmdl.basis) + + for key, op in oldmdl.operations.items(): + # replace each operation `G` with `invU @ G @ U`. + op_s = ComposedOp([U, op, invU]) + mdl.operations[key] = op_s + + mdl._clean_paramvec() # transform may leave dirty members + return mdl + + class ExplicitLayerRules(_LayerRules): """ Rule: layer must appear explicitly as a "primitive op" """ def prep_layer_operator(self, model, layerlbl, caches): diff --git a/pygsti/optimize/optimize.py b/pygsti/optimize/optimize.py index a0bed3ef5..738fea37a 100644 --- a/pygsti/optimize/optimize.py +++ b/pygsti/optimize/optimize.py @@ -149,15 +149,16 @@ def _basin_callback(x, f, accept): elif method == 'evolve': solution = _fmin_evolutionary(fn, x0, num_generations=maxiter, num_individuals=500, printer=printer) -# elif method == 'homebrew': -# solution = fmin_homebrew(fn, x0, maxiter) - else: - #Set options for different algorithms + # Set options for different algorithms opts = {'maxiter': maxiter, 'disp': False} - if method == "BFGS": opts['gtol'] = tol # gradient norm tolerance - elif method == "L-BFGS-B": opts['gtol'] = opts['ftol'] = tol # gradient norm and fractional y-tolerance - elif method == "Nelder-Mead": opts['maxfev'] = maxfev # max fn evals (note: ftol and xtol can also be set) + if method == "BFGS": + opts['gtol'] = tol + elif method == "L-BFGS-B": + opts['gtol'] = opts['ftol'] = tol # gradient norm and fractional y-tolerance + opts.pop('disp') # SciPy deprecated this for L-BFGS-B. + elif method == "Nelder-Mead": + opts['maxfev'] = maxfev # max fn evals (note: ftol and xtol can also be set) solution = _spo.minimize(fn, x0, options=opts, method=method, tol=tol, callback=callback, jac=jac) diff --git a/pygsti/protocols/confidenceregionfactory.py b/pygsti/protocols/confidenceregionfactory.py index e3845f025..85b3bad6b 100644 --- a/pygsti/protocols/confidenceregionfactory.py +++ b/pygsti/protocols/confidenceregionfactory.py @@ -16,6 +16,7 @@ import warnings as _warnings import numpy as _np +import scipy.linalg as _la import scipy.stats as _stats from pygsti import optimize as _opt @@ -654,16 +655,21 @@ def _project_hessian(self, hessian, nongauge_space, gauge_space, gradient=None): # to transform H -> H' in another coordinate system v -> w = B @ v: # v.T @ H @ v = some 2nd deriv = v.T @ B.T @ H' @ B @ v in another basis # so H' = invB.T @ H @ invB - assert(_np.allclose(hessian, hessian.T)) + TOL = 1e-7 + assert(_la.norm(hessian.imag) == 0) + sym_err_abs = _la.norm(hessian - hessian.T) + sym_err_rel = sym_err_abs / _la.norm(hessian) + assert(sym_err_rel < TOL) + hessian += hessian.T + hessian /= 2 invB = _np.concatenate([nongauge_space, gauge_space], axis=1) # takes (nongauge,guage) -> orig coords B = _np.linalg.inv(invB) # takes orig -> (nongauge,gauge) coords Hprime = invB.T @ hessian @ invB - #assert(_np.allclose(Hprime, Hprime.T)) # doesn't handle large magnituge Hessians well - assert(_np.linalg.norm(Hprime - Hprime.T) / _np.linalg.norm(Hprime) < 1e-7) + assert(_la.norm(Hprime.imag) == 0) if gradient is not None: # Check that Hprime is block-diagonal -- off-diag should be ~O(gradient) coupling = Hprime[0:nongauge_space.shape[1], nongauge_space.shape[1]:] - if _np.linalg.norm(coupling) / (1e-6 + _np.linalg.norm(gradient)) > 5: + if _np.linalg.norm(coupling) / (10*TOL + _np.linalg.norm(gradient)) > 5: _warnings.warn("Gauge-nongauge mixed partials have unusually high magnitude: \n" + "|off-diag blk| = %.2g should be ~ |gradient| = %.2g" % (_np.linalg.norm(coupling), _np.linalg.norm(gradient))) @@ -1013,47 +1019,7 @@ def retrieve_profile_likelihood_confidence_intervals(self, label=None, component raise ValueError(("Invalid item label (%s) for computing" % label) + "profile likelihood confidence intervals") - def compute_confidence_interval(self, fn_obj, eps=1e-7, - return_fn_val=False, verbosity=0): - """ - Compute the confidence interval for an arbitrary function. - - This "function", however, must be encapsulated as a - `ModelFunction` object, which allows it to neatly specify - what its dependencies are and allows it to compaute finite- - different derivatives more efficiently. - - Parameters - ---------- - fn_obj : ModelFunction - An object representing the function to evaluate. The - returned confidence interval is based on linearizing this function - and propagating the model-space confidence region. - - eps : float, optional - Step size used when taking finite-difference derivatives of fnOfOp. - - return_fn_val : bool, optional - If True, return the value of fnOfOp along with it's confidence - region half-widths. - - verbosity : int, optional - Specifies level of detail in standard output. - - Returns - ------- - df : float or numpy array - Half-widths of confidence intervals for each of the elements - in the float or array returned by fnOfOp. Thus, shape of - df matches that returned by fnOfOp. - f0 : float or numpy array - Only returned when return_fn_val == True. Value of fnOfOp - at the gate specified by op_label. - """ - - nParams = self.model.num_params - f0 = fn_obj.evaluate(self.model) # function value at "base point" - + def compute_grad_f(self, fn_obj, f0, nParams, eps=1e-7): #Get finite difference derivative gradF that is shape (nParams, ) gradF = _create_empty_grad_f(f0, nParams) @@ -1105,6 +1071,51 @@ def compute_confidence_interval(self, fn_obj, eps=1e-7, assert(_np.linalg.norm(_np.imag(f - f0)) < 1e-12 or _np.iscomplexobj(gradF) ), "gradF seems to be the wrong type!" gradF[igp] = _np.real_if_close(f - f0) / eps + return gradF + + def compute_confidence_interval(self, fn_obj, eps=1e-7, + return_fn_val=False, verbosity=0): + """ + Compute the confidence interval for an arbitrary function. + + This "function", however, must be encapsulated as a + `ModelFunction` object, which allows it to neatly specify + what its dependencies are and allows it to compaute finite- + different derivatives more efficiently. + + Parameters + ---------- + fn_obj : ModelFunction + An object representing the function to evaluate. The + returned confidence interval is based on linearizing this function + and propagating the model-space confidence region. + + eps : float, optional + Step size used when taking finite-difference derivatives of fnOfOp. + + return_fn_val : bool, optional + If True, return the value of fnOfOp along with it's confidence + region half-widths. + + verbosity : int, optional + Specifies level of detail in standard output. + + Returns + ------- + df : float or numpy array + Half-widths of confidence intervals for each of the elements + in the float or array returned by fnOfOp. Thus, shape of + df matches that returned by fnOfOp. + f0 : float or numpy array + Only returned when return_fn_val == True. Value of fnOfOp + at the gate specified by op_label. + """ + + nParams = self.model.num_params + f0 = fn_obj.evaluate(self.model) # function value at "base point" + + #Get finite difference derivative gradF that is shape (nParams, ) + gradF = self.compute_grad_f(fn_obj, f0, nParams, eps) return self._compute_return_from_grad_f(gradF, f0, return_fn_val, verbosity) diff --git a/pygsti/protocols/gst.py b/pygsti/protocols/gst.py index 2353bd31f..a770172d3 100644 --- a/pygsti/protocols/gst.py +++ b/pygsti/protocols/gst.py @@ -39,7 +39,11 @@ from pygsti.modelmembers import operations as _op from pygsti.models import Model as _Model from pygsti.models.explicitmodel import ExplicitOpModel as _ExplicitOpModel -from pygsti.models.gaugegroup import GaugeGroup as _GaugeGroup, GaugeGroupElement as _GaugeGroupElement +from pygsti.models.gaugegroup import ( + GaugeGroup as _GaugeGroup, + GaugeGroupElement as _GaugeGroupElement, + FullGaugeGroupElement as _FullGaugeGroupElement +) from pygsti.objectivefns import objectivefns as _objfns, wildcardbudget as _wild from pygsti.circuits.circuitlist import CircuitList as _CircuitList from pygsti.baseobjs.resourceallocation import ResourceAllocation as _ResourceAllocation @@ -3262,7 +3266,44 @@ def __str__(self): s += " " + "\n ".join(list(self.estimates.keys())) + "\n" s += "\n" return s + + +def _add_param_preserving_gauge_opt(results: ModelEstimateResults, est_key: str, gop_params: GSTGaugeOptSuite, verbosity: int = 0): + """ + This function adds one or more models to the `results[est_key].estimates` dict, in a way + that's similar to _add_gauge_opt. It requires that the model + + results[est_key].models['final iteration estimate'] + + use ComposedState for stateprep and ComposedPOVM for measurement. + It starts by calling _add_gauge_opt with (results, est_key, gop_params, verbosity), + then it extracts the gauge group elements from each step of the gauge-optimization + suite, computes their composition, and finally applies a parameterization-preserving + gauge transformation with that composition. + """ + from pygsti.models.explicitmodel import transform_composed_model + est = results.estimates[est_key] + seed_mdl = est.models['final iteration estimate'] + seed_mdl = _ExplicitOpModel.from_nice_serialization(seed_mdl._to_nice_serialization()) + _add_gauge_opt(results, est_key, gop_params, seed_mdl, unreliable_ops=(), verbosity=verbosity) + # ^ That can convert to whatever parameterization it wants + # It'll write to est._gaugeopt_suite. + for gop_name, gop_dictorlist in est._gaugeopt_suite.gaugeopt_argument_dicts.items(): + if isinstance(gop_dictorlist, list): + ggel_mx = _np.eye(seed_mdl.basis.dim) + for sub_gopdict in gop_dictorlist: + assert isinstance(sub_gopdict, dict) + assert '_gaugeGroupEl' in sub_gopdict + ggel_mx = ggel_mx @ sub_gopdict['_gaugeGroupEl'].transform_matrix + ggel = _FullGaugeGroupElement(ggel_mx) + else: + ggel = gop_dictorlist['_gaugeGroupEl'] + model_implicit_gauge = transform_composed_model(est.models['final iteration estimate'], ggel) + est.models[gop_name] = model_implicit_gauge + return + + class GateSetTomographyCheckpoint(_proto.ProtocolCheckpoint): """ A class for storing intermediate results associated with running @@ -3328,7 +3369,6 @@ def _from_nice_serialization(cls, state): # memo holds already de-serialized ob final_objfn, name) - class StandardGSTCheckpoint(_proto.ProtocolCheckpoint): """ A class for storing intermediate results associated with running diff --git a/pygsti/report/factory.py b/pygsti/report/factory.py index af52d5e60..c8781f53f 100644 --- a/pygsti/report/factory.py +++ b/pygsti/report/factory.py @@ -184,8 +184,7 @@ def _create_master_switchboard(ws, results_dict, confidence_level, Ls = None for results in results_dict.values(): - est_labels = _add_new_estimate_labels(est_labels, results.estimates, - combine_robust) + est_labels = _add_new_estimate_labels(est_labels, results.estimates, combine_robust) loc_Ls = results.circuit_lists['final'].xs \ if isinstance(results.circuit_lists['final'], _PlaquetteGridCircuitStructure) else [0] Ls = _add_new_labels(Ls, loc_Ls) @@ -311,10 +310,8 @@ def _create_master_switchboard(ws, results_dict, confidence_level, else: est_modvi = est - switchBd.objfn_builder[d, i] = est.parameters.get( - 'final_objfn_builder', _objfns.ObjectiveFunctionBuilder.create_from('logl')) - switchBd.objfn_builder_modvi[d, i] = est_modvi.parameters.get( - 'final_objfn_builder', _objfns.ObjectiveFunctionBuilder.create_from('logl')) + switchBd.objfn_builder[d, i] = est.parameters.get('final_objfn_builder', _objfns.ObjectiveFunctionBuilder.create_from('logl')) + switchBd.objfn_builder_modvi[d, i] = _objfns.ObjectiveFunctionBuilder.create_from('logl') switchBd.params[d, i] = est.parameters #add the final mdc store @@ -1233,6 +1230,9 @@ def construct_standard_report(results, title="auto", ws = ws or _ws.Workspace() advanced_options = advanced_options or {} + n_leak = advanced_options.get('n_leak', 0) + # ^ It would be preferable to store n_leak in a Basis object, or something similar. + # We're using this for now since it's simple and gets the job done. linlogPercentile = advanced_options.get('linlog percentile', 5) nmthreshold = advanced_options.get('nmthreshold', DEFAULT_NONMARK_ERRBAR_THRESHOLD) embed_figures = advanced_options.get('embed_figures', True) @@ -1324,7 +1324,11 @@ def construct_standard_report(results, title="auto", try: idt_results = _construct_idtresults(idtIdleOp, idtPauliDicts, results, printer) except Exception as e: - _warnings.warn("Idle tomography failed:\n" + str(e)) + if isinstance(e, ValueError) and 'Expected matrix of shape' in str(e): + msg = "Idle tomography skipped. Currently, this is only supported for 2-level systems." + printer.log(msg) + else: + _warnings.warn("Idle tomography unexpectedly failed:\n" + str(e)) idt_results = {} if len(idt_results) > 0: sections.append(_section.IdleTomographySection()) @@ -1375,7 +1379,8 @@ def construct_standard_report(results, title="auto", 'gauge_opt_labels': tuple(gauge_opt_labels), 'max_lengths': tuple(Ls), 'switchbd_maxlengths': tuple(swLs), - 'show_unmodeled_error': bool('ShowUnmodeledError' in flags) + 'show_unmodeled_error': bool('ShowUnmodeledError' in flags), + 'n_leak' : n_leak } templates = dict( diff --git a/pygsti/report/reportables.py b/pygsti/report/reportables.py index cd1b403ee..909af876c 100644 --- a/pygsti/report/reportables.py +++ b/pygsti/report/reportables.py @@ -25,7 +25,11 @@ from pygsti.report import modelfunction as _modf from pygsti import algorithms as _alg from pygsti import tools as _tools -from pygsti.baseobjs.basis import Basis as _Basis, DirectSumBasis as _DirectSumBasis +from pygsti.baseobjs.basis import ( + Basis as _Basis, + DirectSumBasis as _DirectSumBasis, + BuiltinBasis as _BuiltinBasis +) from pygsti.baseobjs.label import Label as _Lbl from pygsti.baseobjs.errorgenlabel import LocalElementaryErrorgenLabel as _LEEL from pygsti.modelmembers.operations.lindbladcoefficients import LindbladCoefficientBlock as _LindbladCoefficientBlock @@ -1012,6 +1016,88 @@ def maximum_trace_dist(gate, mx_basis): Maximum_trace_dist = _modf.opfn_factory(maximum_trace_dist) # init args == (model, op_label) +def leaky_maximum_trace_dist(gate, mx_basis): + closestUOpMx = _alg.find_closest_unitary_opmx(gate) + _tools.jamiolkowski_iso(closestUOpMx, mx_basis, mx_basis) + n_leak = 1 + return _tools.subspace_jtracedist(gate, closestUOpMx, mx_basis, n_leak) + +Leaky_maximum_trace_dist = _modf.opfn_factory(leaky_maximum_trace_dist) + +def diamonddist_to_leakfree_cptp(op, ignore, mx_basis): + import pygsti.tools.sdptools as _sdps + prob, _, solvers = _sdps.diamond_distance_projection_model( + op, mx_basis, leakfree=True, seepfree=False, n_leak=1, cptp=True, subspace_diamond=False + ) + for s in solvers: + try: + prob.solve(solver=s) + return prob.value + except _sdps.cp.SolverError: + continue + return -1 + +Diamonddist_to_leakfree_cptp = _modf.opsfn_factory(diamonddist_to_leakfree_cptp) + +def subspace_diamonddist_to_leakfree_cptp(op, ignore, mx_basis): + import pygsti.tools.sdptools as _sdps + prob, _, solvers = _sdps.diamond_distance_projection_model( + op, mx_basis, leakfree=True, seepfree=False, n_leak=1, cptp=True, subspace_diamond=True + ) + for s in solvers: + try: + prob.solve(solver=s) + return prob.value + except _sdps.cp.SolverError: + continue + return -1 + +SubspaceDiamonddist_to_leakfree_cptp = _modf.opsfn_factory(subspace_diamonddist_to_leakfree_cptp) + +def subspace_diamonddist(op_a, op_b, basis): + dim_mixed = op_a.shape[0] + dim_pure = int(dim_mixed**0.5) + dim_pure_compsub = dim_pure - 1 + from pygsti.tools.leakage import leading_dxd_submatrix_basis_vectors + U = leading_dxd_submatrix_basis_vectors(dim_pure_compsub, dim_pure, basis) + P = U @ U.T.conj() + assert _np.linalg.norm(P - P.real) < 1e-10 + P = P.real + from pygsti.tools.optools import diamonddist + return diamonddist(op_a @ P, op_b @ P, basis) / 2 + +SubspaceDiamonddist = _modf.opsfn_factory(subspace_diamonddist) + +def pergate_leakrate_reduction(op, ignore, mx_basis, reduction): + assert op.shape == (9, 9) + lfb = _BuiltinBasis('l2p1', 9) + op_lfb = _tools.change_basis(op, mx_basis, lfb) + elinds = lfb.elindlookup + compinds = [elinds[sslbl] for sslbl in ['I','X','Y','Z'] ] + leakage_effect_superket = op_lfb[elinds['L'], compinds] + leakage_effect = _tools.vec_to_stdmx(leakage_effect_superket, 'pp') + leakage_rates = _np.linalg.eigvalsh(leakage_effect) + return reduction(leakage_rates) + +def pergate_leakrate_max(op, ignore, mx_basis): + return pergate_leakrate_reduction(op, ignore, mx_basis, max) + +def pergate_leakrate_min(op, ignore, mx_basis): + return pergate_leakrate_reduction(op, ignore, mx_basis, min) + +PerGateLeakRateMax = _modf.opsfn_factory(pergate_leakrate_max) +PerGateLeakRateMin = _modf.opsfn_factory(pergate_leakrate_min) + +def pergate_seeprate(op, ignore, mx_basis): + assert op.shape == (9, 9) + lfb = _BuiltinBasis('l2p1', 9) + op_lfb = _tools.change_basis(op, mx_basis, lfb) + elinds = lfb.elindlookup + seeprate = op_lfb[elinds['I'], elinds['L']] + return seeprate + +PerGateSeepRate = _modf.opsfn_factory(pergate_seeprate) + def angles_btwn_rotn_axes(model): """ @@ -1083,6 +1169,12 @@ def entanglement_fidelity(a, b, mx_basis): Entanglement_fidelity = _modf.opsfn_factory(entanglement_fidelity) # init args == (model1, model2, op_label) +def subspace_entanglement_fidelity(a, b, mx_basis): + n_leak = 1 + return _tools.subspace_entanglement_fidelity(a, b, mx_basis, n_leak) + +Subspace_entanglement_fidelity = _modf.opsfn_factory(subspace_entanglement_fidelity) + def entanglement_infidelity(a, b, mx_basis): """ @@ -1109,6 +1201,11 @@ def entanglement_infidelity(a, b, mx_basis): Entanglement_infidelity = _modf.opsfn_factory(entanglement_infidelity) # init args == (model1, model2, op_label) +def leaky_entanglement_infidelity(a, b, mx_basis): + return 1 - subspace_entanglement_fidelity(a, b, mx_basis) + +Leaky_entanglement_infidelity = _modf.opsfn_factory(leaky_entanglement_infidelity) + def closest_unitary_fidelity(a, b, mx_basis): # assume vary model1, model2 fixed """ @@ -1175,6 +1272,14 @@ def frobenius_diff(a, b, mx_basis): # assume vary model1, model2 fixed # init args == (model1, model2, op_label) +def leaky_gate_frob_dist(a, b, mx_basis): + n_leak = 1 + return _tools.subspace_superop_fro_dist(a, b, mx_basis, n_leak) + + +Leaky_gate_frob_dist = _modf.opsfn_factory(leaky_gate_frob_dist) + + def jtrace_diff(a, b, mx_basis): # assume vary model1, model2 fixed """ Jamiolkowski trace distance between a and b @@ -1200,6 +1305,12 @@ def jtrace_diff(a, b, mx_basis): # assume vary model1, model2 fixed Jt_diff = _modf.opsfn_factory(jtrace_diff) # init args == (model1, model2, op_label) +def leaky_jtrace_diff(a, b, mx_basis): + n_leak = 1 + return _tools.subspace_jtracedist(a, b, mx_basis, n_leak) + +Leaky_Jt_diff = _modf.opsfn_factory(leaky_jtrace_diff) + if _CVXPY_AVAILABLE: @@ -2435,6 +2546,17 @@ def info_of_opfn_by_name(name): - "frob" : frobenius distance - "unmodeled" : unmodeled "wildcard" budget + - "sub-inf" : subspace entanglement infidelity + - "sub-trace" : subspace trace distance + - "sub-diamond" : subspace diamond distance + - "sub-frob" : subspace Frobenius distance + - "plf-diamond" : diamond distance to the set of leakage-free CPTP maps + - "plf-sub-diamond" : subspace diamond distance to the set of leakage-free CPTP maps + - "leak-rate-max" : maximum transport of population out of the computational subspace + - "leak-rate-min" : minimum transport of population out of the computational subspace, assuming all population starts in that subspace. + - "seep-rate" : maximum transport of population into the computational subspace + + Returns ------- nicename : str @@ -2476,7 +2598,20 @@ def info_of_opfn_by_name(name): "frob": ("Frobenius|Distance", "sqrt( sum( (A_ij - B_ij)^2 ) )"), "unmodeled": ("Un-modeled|Error", - "The per-operation budget used to account for un-modeled errors (model violation)") + "The per-operation budget used to account for un-modeled errors (model violation)"), + "sub-inf": ("Entanglement|Infidelity (subspace)", + "1.0 - , where |psi> is the standard test state supported U x U, and U denotes the computational subspace."), + "sub-trace" : ("1/2 Trace|Distance (subspace)", + "0.5 | Chi(A|U) - Chi(B|U) |_tr, where U is the computational subspace and Chi(G|U) is the image of (1 x G) on the standard test state supported on U x U."), + "sub-diamond": ("1/2 Diamond|Distance (subspace)", + "0.5 sup | (1 x (A-B))(rho) |_tr, rho in subspace"), + "sub-frob": ("Frobenius|Distance (subspace)", + "| (A - B) P |_Fro, where P is the projector onto the computational subspace."), + "plf-diamond" : ("1/2 Diamond|Distance to|No-leak Channel", "0.5 min |A - X|_◆, X is CPTP and leakage-free."), + "plf-sub-diamond" : ("Leakage:|1/2 Min|Diamond|Dist.|to a No-leak|Channel|", "0.5 min |A - X|_{subspace ◆}, X is CPTP and leakage-free."), + "seep-rate": ("Seepage:|Max TOP", "This gate's maximum transport of population INTO the computational subspace. This is only an error metric insofar as seepage implies leakage under unitary dynamics."), + "leak-rate-max": ("Leakage:|Max TOP", "This gate's maximum transport of population OUT of the computational subspace."), + "leak-rate-min": ("Leakage:|Min TOP|(subspace)", "This gate's minimum transport of population OUT of the computational subspace, assuming the entire population starts in the computational subspace.") } if name in info: return info[name] @@ -2520,6 +2655,9 @@ def evaluate_opfn_by_name(name, model, target_model, op_label_or_string, if name == "inf": fn = Entanglement_infidelity if b else \ Circuit_entanglement_infidelity + elif name == "sub-inf": + assert b + fn = Leaky_entanglement_infidelity elif name == "agi": fn = Avg_gate_infidelity if b else \ Circuit_avg_gate_infidelity @@ -2529,9 +2667,30 @@ def evaluate_opfn_by_name(name, model, target_model, op_label_or_string, elif name == "trace": fn = Jt_diff if b else \ Circuit_jt_diff + elif name == "sub-trace": + assert b + fn = Leaky_Jt_diff elif name == "diamond": fn = HalfDiamondNorm if b else \ CircuitHalfDiamondNorm + elif name == 'sub-diamond': + assert b + fn = SubspaceDiamonddist + elif name == 'plf-sub-diamond': + assert b + fn = SubspaceDiamonddist_to_leakfree_cptp + elif name == 'plf-diamond': + assert b + fn = Diamonddist_to_leakfree_cptp + elif name == 'leak-rate-max': + assert b + fn = PerGateLeakRateMax + elif name == 'leak-rate-min': + assert b + fn = PerGateLeakRateMin + elif name == 'seep-rate': + assert b + fn = PerGateSeepRate elif name == "nuinf": fn = Nonunitary_entanglement_infidelity if b else \ Circuit_nonunitary_entanglement_infidelity @@ -2556,6 +2715,9 @@ def evaluate_opfn_by_name(name, model, target_model, op_label_or_string, elif name == "evnudiamond": fn = Eigenvalue_nonunitary_diamondnorm if b else \ Circuit_eigenvalue_nonunitary_diamondnorm + elif name == "sub-frob": + assert b + fn = Leaky_gate_frob_dist elif name == "frob": fn = Fro_diff if b else \ Circuit_fro_diff diff --git a/pygsti/report/section/gauge.py b/pygsti/report/section/gauge.py index 5608d2435..219bba41f 100644 --- a/pygsti/report/section/gauge.py +++ b/pygsti/report/section/gauge.py @@ -167,9 +167,13 @@ def final_model_spam_vs_target_table(workspace, switchboard=None, confidence_lev @_Section.figure_factory(4) def final_gates_vs_target_table_gauge_var(workspace, switchboard=None, confidence_level=None, ci_brevity=1, **kwargs): + if kwargs.get('n_leak', 0) == 0: + display = ('inf', 'agi', 'geni', 'trace', 'diamond', 'nuinf', 'nuagi') + else: + display = ('sub-inf', 'sub-trace', 'sub-diamond', 'plf-sub-diamond', 'leak-rate-max', 'leak-rate-min', 'seep-rate' ) return workspace.GatesVsTargetTable( switchboard.mdl_final, switchboard.mdl_target, _cri(1, switchboard, confidence_level, ci_brevity), - display=('inf', 'agi', 'geni', 'trace', 'diamond', 'nuinf', 'nuagi') + display=display ) @_Section.figure_factory(3) diff --git a/pygsti/report/section/summary.py b/pygsti/report/section/summary.py index 21a64a1cd..4d6764763 100644 --- a/pygsti/report/section/summary.py +++ b/pygsti/report/section/summary.py @@ -38,7 +38,10 @@ def final_model_fit_histogram(workspace, switchboard=None, linlog_percentile=5, @_Section.figure_factory() def final_gates_vs_target_table_insummary(workspace, switchboard=None, confidence_level=None, ci_brevity=1, show_unmodeled_error=False, **kwargs): - summary_display = ('inf', 'trace', 'diamond', 'geni', 'evinf', 'evdiamond') + if kwargs.get('n_leak', 0) == 0: + summary_display = ('inf', 'trace', 'diamond', 'evinf', 'evdiamond') + else: + summary_display = ('sub-inf', 'sub-trace', 'sub-diamond', 'plf-sub-diamond', 'leak-rate-max') wildcardBudget = None if show_unmodeled_error: summary_display += ('unmodeled',) diff --git a/pygsti/report/templates/standard_html_report/tabs/GaugeVariants.html b/pygsti/report/templates/standard_html_report/tabs/GaugeVariants.html index a83d32fc2..2398b2907 100644 --- a/pygsti/report/templates/standard_html_report/tabs/GaugeVariants.html +++ b/pygsti/report/templates/standard_html_report/tabs/GaugeVariants.html @@ -7,8 +7,8 @@

Gauge Variant Error Metrics

-
Individual gate error metrics This table presents various (gauge-variant) metrics that quantify errors in each individual estimated logic gate, with respect to the ideal target gates. Note that "Entanglement infidelity" and "Average gate infidelity" are two common definitions of process fidelity, and related by a constant dimensional factor. A description of each metric can be found by hovering the pointer over the column header.
- {{ final_gates_vs_target_table_gauge_var|render }} +
Individual gate error metrics This table presents various (gauge-variant) metrics that quantify errors in each individual estimated logic gate, with respect to the ideal target gates. Note that "Entanglement infidelity" and "Average gate infidelity" are two common definitions of process fidelity, and related by a constant dimensional factor. A description of each metric can be found by hovering the pointer over the column header.
+ {{ final_gates_vs_target_table_gauge_var|render }}
diff --git a/pygsti/report/templates/standard_pdf_report.tex b/pygsti/report/templates/standard_pdf_report.tex index 9f2b455a9..e75188da0 100644 --- a/pygsti/report/templates/standard_pdf_report.tex +++ b/pygsti/report/templates/standard_pdf_report.tex @@ -409,7 +409,7 @@ \section{Gauge-dependent Error Metrics} \begin{adjustbox}{max width=\textwidth} \putfield{final_gates_vs_target_table_gauge_var}{} \end{adjustbox} - \caption{\textbf{Individual gate error metrics} This table presents various (gauge-variant) metrics that quantify errors in each individual estimated logic gate, with respect to the ideal target gates. Note that "Entanglement infidelity" and "Average gate infidelity" are two common definitions of process fidelity, and related by a constant dimensional factor. A description of each metric can be found by hovering the pointer over the column header.} + \caption{\textbf{Individual gate error metrics} This table presents various (gauge-variant) metrics that quantify errors in each individual estimated logic gate, with respect to the ideal target gates. Note that "Entanglement infidelity" and "Average gate infidelity" are two common definitions of process fidelity, and related by a constant dimensional factor. A description of each metric can be found by hovering the pointer over the column header.} \end{center} \end{table} diff --git a/pygsti/tools/__init__.py b/pygsti/tools/__init__.py index 52829ca9a..b631e8ea3 100644 --- a/pygsti/tools/__init__.py +++ b/pygsti/tools/__init__.py @@ -19,6 +19,7 @@ from .jamiolkowski import * from .legacytools import * from .likelihoodfns import * +from .leakage import * from .lindbladtools import * from .listtools import * from .matrixmod2 import * diff --git a/pygsti/tools/jamiolkowski.py b/pygsti/tools/jamiolkowski.py index c53ec9d7a..ea964c396 100644 --- a/pygsti/tools/jamiolkowski.py +++ b/pygsti/tools/jamiolkowski.py @@ -10,12 +10,20 @@ # http://www.apache.org/licenses/LICENSE-2.0 or in the LICENSE file in the root pyGSTi directory. #*************************************************************************************************** +from __future__ import annotations + import numpy as _np from pygsti.tools import basistools as _bt -from pygsti.tools import matrixtools as _mt from pygsti.baseobjs.basis import Basis as _Basis -from pygsti import SpaceT +from typing import Union, TYPE_CHECKING + +if TYPE_CHECKING: + from cvxpy import Expression + +BasisLike = Union[_Basis, str] + + # Gate Mx G: rho --> G rho where G and rho are in the Pauli basis (by definition/convention) # noqa # vec(rhoS) --> GStd vec(rhoS) where GS and rhoS are in the std basis, GS = PtoS * G * StoP # noqa @@ -64,13 +72,14 @@ # J(Phi) = sum_(0 Union[_np.ndarray, Expression]: """ - Given a operation matrix, return the corresponding Choi matrix that is normalized to have trace == 1. + Return a Choi matrix (in the choi_mx_basis) for operation_mx, when operation_mx + is interpreted in the op_mx_basis. Parameters ---------- - operation_mx : numpy array + operation_mx : numpy array or cvxpy Expression the operation matrix to compute Choi matrix of. op_mx_basis : Basis object @@ -83,12 +92,23 @@ def jamiolkowski_iso(operation_mx, op_mx_basis='pp', choi_mx_basis='pp'): values are Matrix-unit (std), Gell-Mann (gm), Pauli-product (pp), and Qutrit (qt) (or a custom basis object). + normalized : bool + If normalized=True, then this function maps trace-preserving + operation matrices to trace-1 Choi matrices. + Returns ------- - numpy array - the Choi matrix, normalized to have trace == 1, in the desired basis. + numpy array or cvxpy Expression + the Choi matrix, in the desired basis. """ - operation_mx = _np.asarray(operation_mx) + try: + import cvxpy as cp + is_cvxpy_expression = isinstance(operation_mx, cp.Expression) + except ImportError: + is_cvxpy_expression = False + + if not is_cvxpy_expression: + operation_mx = _np.asarray(operation_mx) op_mx_basis = _bt.create_basis_for_matrix(operation_mx, op_mx_basis) opMxInStdBasis = _bt.change_basis(operation_mx, op_mx_basis, op_mx_basis.create_equivalent('std')) @@ -113,25 +133,33 @@ def jamiolkowski_iso(operation_mx, op_mx_basis='pp', choi_mx_basis='pp'): M = len(BVec) # can be < N if basis has multiple block dims assert(M == N), 'Expected {}, got {}'.format(M, N) - choiMx = _np.empty((N, N), 'complex') + opMxInStdBasis_vec = opMxInStdBasis.flatten() + # ^ use flatten, not ravel, in case we're using a CVXPY Expression. + choiMx_rows = [] for i in range(M): + rows = [] for j in range(M): BiBj = _np.kron(BVec[i], _np.conjugate(BVec[j])) - num = _np.vdot(BiBj, opMxInStdBasis) - den = _np.linalg.norm(BiBj) ** 2 - choiMx[i, j] = num / den - + BiBj /= _np.linalg.norm(BiBj) ** 2 + rows.append(BiBj.conj().ravel()) + choiMx_rows.append( _np.array(rows) @ opMxInStdBasis_vec ) + if is_cvxpy_expression: + choiMx = cp.vstack(choiMx_rows) + else: + choiMx = _np.vstack(choiMx_rows) # This construction results in a Jmx with trace == dim(H) = sqrt(operation_mx.shape[0]) # (dimension of density matrix) but we'd like a Jmx with trace == 1, so normalize: - choiMx_normalized = choiMx / dmDim - return choiMx_normalized + if normalized: + choiMx /= dmDim + return choiMx # GStd = sum_ij Jij (BSi x BSj^*) -def jamiolkowski_iso_inv(choi_mx, choi_mx_basis='pp', op_mx_basis='pp'): +def jamiolkowski_iso_inv(choi_mx: Union[_np.ndarray, Expression], choi_mx_basis: BasisLike ='pp', op_mx_basis: BasisLike ='pp', normalized: bool = True) -> Union[_np.ndarray, Expression]: """ - Given a choi matrix, return the corresponding operation matrix. + Given a choi matrix (interpreted in choi_mx_basis), return the corresponding + operation matrix (in op_mx_basis). This function performs the inverse of :func:`jamiolkowski_iso`. @@ -150,6 +178,10 @@ def jamiolkowski_iso_inv(choi_mx, choi_mx_basis='pp', op_mx_basis='pp'): values are Matrix-unit (std), Gell-Mann (gm), Pauli-product (pp), and Qutrit (qt) (or a custom basis object). + normalized : bool + If normalized=True, then we assume choi_mx was computed with the + convention that trace-preserving maps have trace-1 choi matrices. + Returns ------- numpy array @@ -167,7 +199,10 @@ def jamiolkowski_iso_inv(choi_mx, choi_mx_basis='pp', op_mx_basis='pp'): assert(len(BVec) == N) # make sure the number of basis matrices matches the dim of the choi matrix given # Invert normalization - choiMx_unnorm = choi_mx * dmDim + if normalized: + choiMx_unnorm = choi_mx * dmDim + else: + choiMx_unnorm = choi_mx opMxInStdBasis = _np.zeros((N, N), 'complex') # in matrix unit basis of entire density matrix for i in range(N): @@ -187,13 +222,14 @@ def jamiolkowski_iso_inv(choi_mx, choi_mx_basis='pp', op_mx_basis='pp'): return _bt.change_basis(opMxInStdBasis, op_mx_basis.create_equivalent('std'), op_mx_basis) -def fast_jamiolkowski_iso_std(operation_mx, op_mx_basis): +def fast_jamiolkowski_iso_std(operation_mx: _np.ndarray, op_mx_basis: BasisLike, normalized: bool = True) -> _np.ndarray: """ - The corresponding Choi matrix in the standard basis that is normalized to have trace == 1. + Returns the standard-basis representation of the Choi matrix for operation_mx, + where operation_mx is interpreted in op_mx_basis. This routine *only* computes the case of the Choi matrix being in the standard (matrix unit) basis, but does so more quickly than - :func:`jamiolkowski_iso` and so is particuarly useful when only the + :func:`jamiolkowski_iso` and so is particularly useful when only the eigenvalues of the Choi matrix are needed. Parameters @@ -206,6 +242,10 @@ def fast_jamiolkowski_iso_std(operation_mx, op_mx_basis): values are Matrix-unit (std), Gell-Mann (gm), Pauli-product (pp), and Qutrit (qt) (or a custom basis object). + normalized : bool + If normalized=True, then this function maps trace-preserving + operation matrices to trace-1 Choi matrices. + Returns ------- numpy array @@ -215,11 +255,13 @@ def fast_jamiolkowski_iso_std(operation_mx, op_mx_basis): #first, get operation matrix into std basis operation_mx = _np.asarray(operation_mx) op_mx_basis = _bt.create_basis_for_matrix(operation_mx, op_mx_basis) - opMxInStdBasis = _bt.change_basis(operation_mx, op_mx_basis, op_mx_basis.create_equivalent('std')) + temp = op_mx_basis.create_equivalent('std') + opMxInStdBasis = _bt.change_basis(operation_mx, op_mx_basis, temp) #expand operation matrix so it acts on entire space of dmDim x dmDim density matrices - opMxInStdBasis = _bt.resize_std_mx(opMxInStdBasis, 'expand', op_mx_basis.create_equivalent('std'), - op_mx_basis.create_simple_equivalent('std')) + temp1 = op_mx_basis.create_equivalent('std') + temp2 = op_mx_basis.create_simple_equivalent('std') + opMxInStdBasis = _bt.resize_std_mx(opMxInStdBasis, 'expand',temp1, temp2) #Shuffle indices to go from process matrix to Jamiolkowski matrix (they vectorize differently) N2 = opMxInStdBasis.shape[0]; N = int(_np.sqrt(N2)) @@ -230,13 +272,15 @@ def fast_jamiolkowski_iso_std(operation_mx, op_mx_basis): # This construction results in a Jmx with trace == dim(H) = sqrt(gateMxInPauliBasis.shape[0]) # but we'd like a Jmx with trace == 1, so normalize: - Jmx_norm = Jmx / N - return Jmx_norm + if normalized: + Jmx /= N + return Jmx -def fast_jamiolkowski_iso_std_inv(choi_mx, op_mx_basis): +def fast_jamiolkowski_iso_std_inv(choi_mx: _np.ndarray, op_mx_basis: BasisLike, normalized: bool = True) -> _np.ndarray: """ - Given a choi matrix in the standard basis, return the corresponding operation matrix. + Given a choi matrix in the standard basis, return the corresponding + operation matrix (in op_mx_basis). This function performs the inverse of :func:`fast_jamiolkowski_iso_std`. @@ -251,6 +295,10 @@ def fast_jamiolkowski_iso_std_inv(choi_mx, op_mx_basis): values are Matrix-unit (std), Gell-Mann (gm), Pauli-product (pp), and Qutrit (qt) (or a custom basis object). + normalized : bool + If normalized=True, then we assume choi_mx was computed with the + convention that trace-preserving maps have trace-1 choi matrices. + Returns ------- numpy array @@ -260,7 +308,10 @@ def fast_jamiolkowski_iso_std_inv(choi_mx, op_mx_basis): #Shuffle indices to go from process matrix to Jamiolkowski matrix (they vectorize differently) N2 = choi_mx.shape[0]; N = int(_np.sqrt(N2)) assert(N * N == N2) # make sure N2 is a perfect square - opMxInStdBasis = choi_mx.reshape((N, N, N, N)) * N + opMxInStdBasis = choi_mx.reshape((N, N, N, N)) + if normalized: + opMxInStdBasis = N * opMxInStdBasis + # ^ Don't do this in-place. opMxInStdBasis = _np.swapaxes(opMxInStdBasis, 1, 2).ravel() opMxInStdBasis = opMxInStdBasis.reshape((N2, N2)) op_mx_basis = _bt.create_basis_for_matrix(opMxInStdBasis, op_mx_basis) @@ -299,7 +350,7 @@ def sum_of_negative_choi_eigenvalues(model, weights=None): """ Compute the amount of non-CP-ness of a model. - This is defined (somewhat arbitarily) by summing the negative + This is defined (somewhat arbitrarily) by summing the negative eigenvalues of the Choi matrix for each gate in `model`. Parameters @@ -330,7 +381,7 @@ def sums_of_negative_choi_eigenvalues(model): """ Compute the amount of non-CP-ness of a model. - This is defined (somewhat arbitarily) by summing the negative + This is defined (somewhat arbitrarily) by summing the negative eigenvalues of the Choi matrix for each gate in model separately. This function is different from :func:`sum_of_negative_choi_eigenvalues` in that it returns sums separately for each operation of `model`. @@ -359,7 +410,7 @@ def sums_of_negative_choi_eigenvalues(model): def magnitudes_of_negative_choi_eigenvalues(model): """ - Compute the magnitudes of the negative eigenvalues of the Choi matricies for each gate in `model`. + Compute the magnitudes of the negative eigenvalues of the Choi matrices for each gate in `model`. Parameters ---------- diff --git a/pygsti/tools/leakage.py b/pygsti/tools/leakage.py new file mode 100644 index 000000000..d5dbe8f31 --- /dev/null +++ b/pygsti/tools/leakage.py @@ -0,0 +1,481 @@ +#*************************************************************************************************** +# Copyright 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. +#*************************************************************************************************** + +from __future__ import annotations + +from functools import lru_cache + +import copy +from pygsti.tools import optools as pgot +from pygsti.tools import basistools as pgbt +from pygsti.tools.basistools import stdmx_to_vec +from pygsti.processors import QubitProcessorSpec +from pygsti.baseobjs.basis import TensorProdBasis, Basis, BuiltinBasis +import numpy as np +import warnings + +from typing import Union, Dict, Optional, List, Any, TYPE_CHECKING +if TYPE_CHECKING: + from pygsti.protocols.gst import ModelEstimateResults, GSTGaugeOptSuite + from pygsti.models import ExplicitOpModel + + +# Question: include the parenthetical in the heading below? +NOTATION = \ +""" +Default notation (deferential to text above) +-------------------------------------------- + * H is a complex Hilbert space equipped with the standard basis. + + * C, the computational subspace, is the complex-linear span of the first dim(C) + standard basis vectors of H. + + * Given a complex Hilbert space, U, we write M[U] to denote the space of linear + operators from U to U. Elements of M[U] have natural matrix representations. + + * Given a space of linear operators, L, we write S[L] for the set of linear + transformations ("superoperators") from L to L. + + Matrix representations for elements of S[L] are only meaningful in the + presence of a designated basis for L. + + If elements of L are naturally expressed as matrices, then a basis for L + lets us identify elements of L with vectors of length dim(L). + + * If U denotes a complex Hilbert space (e.g., U=H or U=C), then we abbreviate + S[M[U]] by S[U]. +\n +""" + +def set_docstring(docstr): + def assign(fn): + fn.__doc__ = docstr + return fn + return assign + + +# MARK: metrics + +@lru_cache +@set_docstring( +""" +Here, H has dimension dim, and C ⊂ H has co-dimension n_leak. + +This function returns a rank-1 density matrix, rho_t = |psi> np.ndarray: + temp = np.eye(dim, dtype=np.complex128) + if n_leak > 0: + temp[-n_leak:,-n_leak:] = 0.0 + temp /= np.sqrt(dim - n_leak) + psi = pgbt.stdmx_to_stdvec(temp).ravel() + # + # |psi> = (|00> + |11> + ... + |dim - n_leak - 1>) / sqrt(dim - n_leak). + # + rho_test = np.outer(psi, psi) + return rho_test + + +@set_docstring( +""" +The pair (op_x, op_y) represent some superoperators (X, Y) in S[H], using op_basis. + +Let rho_t = tensorized_teststate_density(dim(H), n_leak), and set I to the identity +operator in S[H]. + +This function returns a triplet, consisting of + * tb : a tensor product basis for S[H]⨂S[H], + * a vector representation of X⨂I(rho_t) in the basis tb, and + * a vector representation of Y⨂I(rho_t) in the basis tb. + +*Warning!* At present, this function can only be used for gates over a single system +(e.g., a single qubit), not for tensor products of such systems. +""" + NOTATION) +def apply_tensorized_to_teststate(op_x: np.ndarray, op_y, op_basis: np.ndarray, n_leak: int=0) -> tuple[TensorProdBasis, np.ndarray, np.ndarray]: + dim = int(np.sqrt(op_x.shape[0])) + assert op_x.shape == (dim**2, dim**2) + assert op_y.shape == (dim**2, dim**2) + # op_x and op_y act on M[H]. + # + # We care about op_x and op_y only up to their action on the subspace + # U = {rho in M[H] : = 0 for all i >= dim - n_leak }. + # + # It's easier to talk about this subspace (and related subspaces) if op_x and op_y are in + # the standard basis. So the first thing we do is convert to that basis. + std_basis = BuiltinBasis('std', dim**2) + op_x_std = pgbt.change_basis(op_x, op_basis, std_basis) + op_y_std = pgbt.change_basis(op_y, op_basis, std_basis) + + # Our next step is to construct lifted operators "lift_op_x" and "lift_op_y" that act on the + # tensor product space M[H]⨂M[H] according to the identities + # + # lift_op_x( sigma \otimes rho ) = op_x(sigma) \otimes rho + # lift_op_y( sigma \otimes rho ) = op_y(sigma) \otimes rho + # + # for all sigma, rho in M[H]. The way we do this implicitly fixes a basis for S[H]⨂S[H] as + # the tensor product basis. We'll make that explicit later on. + idle_gate = np.eye(dim**2, dtype=np.complex128) + lift_op_x_std = np.kron(op_x_std, idle_gate) + lift_op_y_std = np.kron(op_y_std, idle_gate) + + # Now we'll compare these lifted operators by how they act on specific state in M[H]⨂M[H]. + rho_test = tensorized_teststate_density(dim, n_leak) + + # lift_op_x and lift_op_y only act on states in their superket representations, so we convert + # rho_test to a superket representation in the induced tensor product basis for S[H]⨂S[H]. + ten_std_basis = TensorProdBasis((std_basis, std_basis)) + rho_test_superket = pgbt.stdmx_to_vec(rho_test, ten_std_basis).ravel() + + temp1 = lift_op_x_std @ rho_test_superket + temp2 = lift_op_y_std @ rho_test_superket + + return ten_std_basis, temp1, temp2 + + +@lru_cache +def leading_dxd_submatrix_basis_vectors(d: int, n: int, current_basis: Basis): + """ + Let "H" denote n^2 dimensional Hilbert-Schdmit space, and let "U" denote the d^2 + dimensional subspace of H spanned by vectors whose Hermitian matrix representations + are zero outside the leading d-by-d submatrix. + + This function returns a column-unitary matrix "B" where P = B B^{\\dagger} is the + orthogonal projector from H to U with respect to current_basis. We return B rather + than P only because it's simpler to get P from B than it is to get B from P. + + See below for this function's original use-case. + + Raison d'etre + ------------- + Suppose we canonically measure the distance between two process matrices (M1, M2) by + + D(M1, M2; H) = max || (M1 - M2) v || + v is in H, (Eq. 1) + tr(v) = 1, + v is positive + + for some norm || * ||. Suppose also that we want an analog of this distance when + (M1, M2) are restricted to the linear subspace U consisting of all vectors in H + whose matrix representations are zero outside of their leading d-by-d submatrix. + + One natural way to do this is via the function D(M1, M2; U) -- i.e., just replace + H in (Eq. 1) with the subspace U. Using P to denote the orthogonal projector onto U, + we claim that we can evaluate this function via the identity + + D(M1, M2; U) = D(M1 P, M2 P; H). (Eq. 2) + + To see why this is the case, consider a positive vector v and its projection u = P v. + Since a vector is called positive whenever its Hermitian matrix representation is positive + semidefinite (PSD), we need to show that u is positive. This can be seen by considering + block 2-by-2 partitions of the matrix representations of (u,v), where the leading block + is d-by-d: + + mat(v) = [x11, x12] and mat(u) = [x11, 0] + [x21, x22] [ 0, 0]. + + In particular, u is positive if and only if x11 is PSD, and x11 must be PSD for v + to be positive. Furthermore, positivity of v requires that x22 is PSD, which implies + + 0 <= tr(u) = tr(x11) <= tr(v). + + Given this, it is easy to establish (Eq 2.) by considering how the following pair + of problems have the same optimal objective function value + + max || (M1 - M2) P v || and max || (M1 - M2) P v || + mat(v) = [x11, x12] mat(v) = [x11, x12] + [x21, x22] [x21, x22] + mat(v) is PSD x11 is PSD + tr(x11) + tr(x22) = 1 tr(x11) <= 1. + + In fact, this can be taken a little further! The whole argument goes through unchanged + if, instead of starting with the objective function || (M1 - M2) v ||, we started with + f((M1 - M2) v) and f satisfied the property that f(c v) >= f(v) whenever c is a scalar + greater than or equal to one. + """ + assert d <= n + current_basis = Basis.cast(current_basis, dim=n**2) + X = current_basis.create_transform_matrix('std') + X = X.T.conj() + if d == n: + return X + # we have to select a proper subset of columns in current_basis + std_basis = BuiltinBasis(name='std', dim_or_statespace=n**2) + label2ind = std_basis.elindlookup + basis_ind = [] + for i in range(d): + for j in range(d): + ell = f"({i},{j})" + basis_ind.append(label2ind[ell]) + basis_ind = np.array(basis_ind) + submatrix_basis_vectors = X[:, basis_ind] + return submatrix_basis_vectors + + +CHOI_INDUCED_METRIC = \ +""" +The pair (op_x, op_y) represent some superoperators (X, Y) in S[H], using op_basis. + +Let rho_t = tensorized_teststate_density(dim(H), n_leak), and set I to the identity +operator in S[H]. + +This function returns the %s between X⨂I(rho_t) and Y⨂I(rho_t). + +*Warning!* At present, this function can only be used for gates over a single system +(e.g., a single qubit), not for tensor products of such systems. +""" + NOTATION + + +@set_docstring(CHOI_INDUCED_METRIC % 'entanglement fidelity') +def subspace_entanglement_fidelity(op_x, op_y, op_basis, n_leak=0): + ten_std_basis, temp1, temp2 = apply_tensorized_to_teststate(op_x, op_y, op_basis, n_leak) + temp1_mx = pgbt.vec_to_stdmx(temp1, ten_std_basis, keep_complex=True) + temp2_mx = pgbt.vec_to_stdmx(temp2, ten_std_basis, keep_complex=True) + ent_fid = pgot.fidelity(temp1_mx, temp2_mx) + return ent_fid + + +@set_docstring(CHOI_INDUCED_METRIC % 'jamiolkowski trace distance') +def subspace_jtracedist(op_x, op_y, op_basis, n_leak=0): + ten_std_basis, temp1, temp2 = apply_tensorized_to_teststate(op_x, op_y, op_basis, n_leak) + temp1_mx = pgbt.vec_to_stdmx(temp1, ten_std_basis, keep_complex=True) + temp2_mx = pgbt.vec_to_stdmx(temp2, ten_std_basis, keep_complex=True) + j_dist = pgot.tracedist(temp1_mx, temp2_mx) + return j_dist + + +@set_docstring( +""" +The pair (op_x, op_y) represent some superoperators (X, Y) in S[H], using op_basis. + +We return the Frobenius distance between op_x @ P and op_y @ P, where P is the +projector onto the computational subspace (i.e., C), of co-dimension n_leak. + +*Warning!* At present, this function can only be used for gates over a single system +(e.g., a single qubit), not for tensor products of such systems. +""" + NOTATION) +def subspace_superop_fro_dist(op_x, op_y, op_basis, n_leak=0): + diff = op_x - op_y + if n_leak == 0: + return np.linalg.norm(diff, 'fro') + + d = int(np.sqrt(op_x.shape[0])) + assert op_x.shape == op_y.shape == (d**2, d**2) + B = leading_dxd_submatrix_basis_vectors(d-n_leak, d, op_basis) + P = B @ B.T.conj() + return np.linalg.norm(diff @ P) + + +# MARK: model construction + +def leaky_qubit_model_from_pspec(ps_2level: QubitProcessorSpec, mx_basis: Union[str, Basis]='l2p1', levels_readout_zero=(0,)) -> ExplicitOpModel: + """ + Return an ExplicitOpModel `m` whose (ideal) gates act on three-dimensional Hilbert space and whose members + are represented in `mx_basis`, constructed as follows: + + The Hermitian matrix representation of m['rho0'] is the 3-by-3 matrix with a 1 in the upper-left + corner and all other entries equal to zero. + + Operations in `m` are defined by taking each 2-by-2 unitary `u2` from ps_2level, and promoting it + to a 3-by-3 unitary according to + + u3 = [u2[0, 0], u2[0, 1], 0] + [u2[1, 0], u2[1, 1], 0] + [ 0, 0, 1] + + m['Mdefault'] has two effects, labeled "0" and "1". If E0 is the Hermitian matrix representation of + effect "0", then E0[i,i]=1 for all i in levels_readout_zero, and E0 is zero in all other components. + + This function might be called in a workflow like the following: + + from pygsti.models import create_explicit_model + from pygsti.algorithms import find_fiducials, find_germs + from pygsti.protocols import StandardGST, StandardGSTDesign, ProtocolData + + # Step 1: Make the experiment design for the 1-qubit system. + tm_2level = create_explicit_model( ps_2level, ideal_spam_type='CPTPLND', ideal_gate_type='CPTPLND' ) + fids = find_fiducials( tm_2level ) + germs = find_germs( tm_2level ) + lengths = [1, 2, 4, 8, 16, 32] + design = StandardGSTDesign( tm_2level, fids[0], fids[1], germs, lengths ) + + # Step 2: ... run the experiment specified by "design"; store results in a directory "dir" ... + + # Step 3: read in the experimental data and run GST. + pd = ProtocolData.from_dir(dir) + tm_3level = leaky_qubit_model_from_pspec( ps_2level, basis='l2p1' ) + gst = StandardGST( modes=('CPTPLND',), target_model=tm_3level, verbosity=4 ) + res = gst.run(pd) + """ + from pygsti.models.explicitmodel import ExplicitOpModel + from pygsti.baseobjs.statespace import ExplicitStateSpace + from pygsti.modelmembers.povms import UnconstrainedPOVM + from pygsti.modelmembers.states import FullState + assert ps_2level.num_qubits == 1 + + if isinstance(mx_basis, str): + mx_basis = BuiltinBasis(mx_basis, 9) + assert isinstance(mx_basis, Basis) + + Us = ps_2level.gate_unitaries + rho0 = np.array([[1, 0, 0], [0, 0, 0], [0, 0, 0]], complex) + E0 = np.zeros((3, 3)) + E0[levels_readout_zero, levels_readout_zero] = 1 + E1 = np.eye(3, dtype=complex) - E0 + + ss = ExplicitStateSpace([0],[3]) + tm_3level = ExplicitOpModel(ss, mx_basis) # type: ignore + tm_3level.preps['rho0'] = FullState(stdmx_to_vec(rho0, mx_basis)) + tm_3level.povms['Mdefault'] = UnconstrainedPOVM( + [("0", stdmx_to_vec(E0, mx_basis)), ("1", stdmx_to_vec(E1, mx_basis))], evotype="default", + ) + + def u2x2_to_9x9_superoperator(u2x2): + u3x3 = np.eye(3, dtype=np.complex128) + u3x3[:2,:2] = u2x2 + superop_std = pgot.unitary_to_std_process_mx(u3x3) + superop = pgbt.change_basis(superop_std, 'std', mx_basis) + return superop + + for gatename, unitary in Us.items(): + gatekey = (gatename, 0) if gatename != '{idle}' else ('Gi',0) + tm_3level.operations[gatekey] = u2x2_to_9x9_superoperator(unitary) + + tm_3level.sim = 'map' # can use 'matrix', if that's preferred for whatever reason. + return tm_3level + + +# MARK: gauge optimization + +def lagoified_gopparams_dicts(gopparams_dicts: List[Dict]) -> List[Dict]: + """ + Return a modified deep-copy of gopparams_dicts, where the modifications + strip out gauge optimization over the TPSpam gauge group and apply a + few other changes appropriate for leakage modeling. + """ + gopparams_dicts = [gp for gp in gopparams_dicts if 'TPSpam' not in str(type(gp['_gaugeGroupEl']))] + gopparams_dicts = copy.deepcopy(gopparams_dicts) + for inner_dict in gopparams_dicts: + inner_dict['n_leak'] = 1 + # ^ This function could accept n_leak as an argument instead. However, + # downstream functions for gauge optimization only support n_leak=0 or 1. + # + # When n_leak=1 we use subspace-restricted loss functions that only care + # about mismatches between an estimate and a target when restricted to the + # computational subspace. We have code for evaluating the loss functions + # themselves, but not their gradients. + inner_dict['method'] = 'L-BFGS-B' + # ^ We need this optimizer because it doesn't require a gradient oracle. + inner_dict['gates_metric'] = 'frobenius squared' + inner_dict['spam_metric'] = 'frobenius squared' + # ^ We have other metrics for leakage-aware gauge optimization, but they + # aren't really useful. + inner_dict['convert_model_to']['to_type'] = 'full' + # ^ The natural basis for Hilbert-Schmidt space in leakage modeling doesn't + # have the identity matrix as its first element. This means we can't use + # the full TP parameterization. There's no real harm in taking "full" as + # our default because add_lago_models uses parameterization-preserving + # gauge optimization. + return gopparams_dicts + + +def std_lago_gopsuite(model: ExplicitOpModel) -> dict[str, list[dict]]: + """ + Return a dictionary of the form {'LAGO': v}, where v is a + list-of-dicts representation of a gauge optimization suite. + + We construct v by getting the list-of-dicts representation of the + "stdgaugeopt" suite for `model`, and then changing some of its + options to be suitable for leakage-aware gauge optimization. These + changes are made in the `lagoified_gopparams_dicts` function. + """ + from pygsti.protocols.gst import GSTGaugeOptSuite + std_gop_suite = GSTGaugeOptSuite(gaugeopt_suite_names=('stdgaugeopt',)) + std_gos_lods = std_gop_suite.to_dictionary(model)['stdgaugeopt'] # list of dictionaries + lago_gos_lods = lagoified_gopparams_dicts(std_gos_lods) + gop_params = {'LAGO': lago_gos_lods} + return gop_params + + +def add_lago_models(results: ModelEstimateResults, est_key: Optional[str] = None, gos: Optional[GSTGaugeOptSuite] = None, verbosity: int = 0): + """ + Update each estimate in results.estimates (or just results.estimates[est_key], + if est_key is not None) with a model obtained by parameterization-preserving + gauge optimization. + + If no gauge optimization suite is provided, then we construct one by starting + with "stdgaugeopt" suite and then changing the settings to be suitable for + leakage-aware gauge optimization. + """ + from pygsti.protocols.gst import GSTGaugeOptSuite, _add_param_preserving_gauge_opt + if gos is None: + existing_est = results.estimates[est_key] + std_gos_lods = existing_est.goparameters['stdgaugeopt'] + lago_gos_lods = lagoified_gopparams_dicts(std_gos_lods) + gop_params = {'LAGO': lago_gos_lods} + gos = GSTGaugeOptSuite(gaugeopt_argument_dicts=gop_params) + if isinstance(est_key, str): + _add_param_preserving_gauge_opt(results, est_key, gos, verbosity) + elif est_key is None: + for est_key in results.estimates.keys(): + add_lago_models(results, est_key, gos, verbosity) + else: + raise ValueError() + return + + +# MARK: reports + +def construct_leakage_report( + results : ModelEstimateResults, + title : str = 'auto', + extra_report_kwargs : Optional[dict[str,Any]] = None, + gaugeopt_verbosity : int = 0, + ): + """ + This is a small wrapper around construct_standard_report. It generates a Report object + with leakage analysis, and returns that object along with a copy of ``results`` which + contains gauge-optimized models created during leakage analysis. + + Notes + ----- + The special gauge optimization performed by this function uses the unitary gauge group, + and uses a modified version of the Frobenius distance loss function. The modification + reflects how the target gates in a leakage model are only _really_ defined on the + computational subspace. + """ + + if extra_report_kwargs is None: + extra_report_kwargs = {'title': title} + + if extra_report_kwargs.get('title', title) != title: + # Yes, we let you pass keyword arguments through to construct_standard_report, + # but that doesn't mean you should use that ability to pass in contradictory + # arguments to this function. + ktitle = extra_report_kwargs['title'] + msg = f"Replacing report title in extra_report_kwargs ({ktitle}) " + msg += f"with this function's title argument ({title})." + warnings.warn(msg) + extra_report_kwargs['title'] = title + + results = copy.deepcopy(results) + est_key = results.estimates.keys() + for ek in est_key: + assert isinstance(ek, str) + add_lago_models(results, ek, verbosity=gaugeopt_verbosity) + from pygsti.report import construct_standard_report + report = construct_standard_report( + results, advanced_options={'n_leak': 1}, **extra_report_kwargs + ) + return report, results + +print() diff --git a/pygsti/tools/matrixtools.py b/pygsti/tools/matrixtools.py index c6fbc29c9..210ef23cd 100644 --- a/pygsti/tools/matrixtools.py +++ b/pygsti/tools/matrixtools.py @@ -13,6 +13,7 @@ import functools as _functools import itertools as _itertools import warnings as _warnings +import typing as _typing import numpy as _np import scipy.linalg as _spl @@ -2394,3 +2395,42 @@ def sign_fix_qr(q, r, tol=1e-6): qq[:, i] = -q[:, i] rr[i, :] = -r[i, :] return qq, rr + + +def is_operatorlike(obj: _typing.Any) -> bool: + return hasattr(obj, '__matmul__') and hasattr(obj, '__rmatmul__') and hasattr(obj, 'T') and hasattr(obj, 'conj') + + +class IdentityOperator: + """ + A representation of the identity operator on any and all vector spaces. + """ + + __array_priority__ = 101 + # ^ The __array_priority__ static class variable determines how infix + # operators (most notably, @) are dispatched. + # + # If Python initially dispatches the ndarray implementation of an infix + # operator, then the ndarray implementation will first check if the other + # operand has __array_priority__ greater than zero. If it does, then a + # function call like ndarray.__matmul__(array, other) will return the + # result of other.__rmatmul__(array). + # + # For our purposes, it should suffice to set __array_priority__ to 1. + # We set it to 101 in case someone does something weird and passes in + # a numpy matrix (which behaves like an ndarray in many respects, and + # has __array_priority__ of 10). + # + + def __matmul__(self, other: _typing.Any) -> _typing.Any: + return other + + def __rmatmul__(self, other: _typing.Any) -> _typing.Any: + return other + + @property + def T(self): + return self + + def conj(self): + return self diff --git a/pygsti/tools/optools.py b/pygsti/tools/optools.py index a4f2feec1..5dab09e8a 100644 --- a/pygsti/tools/optools.py +++ b/pygsti/tools/optools.py @@ -10,7 +10,8 @@ # http://www.apache.org/licenses/LICENSE-2.0 or in the LICENSE file in the root pyGSTi directory. #*************************************************************************************************** -import collections as _collections +from __future__ import annotations + import warnings as _warnings import numpy as _np @@ -23,12 +24,22 @@ from pygsti.tools import jamiolkowski as _jam from pygsti.tools import lindbladtools as _lt from pygsti.tools import matrixtools as _mt -from pygsti.baseobjs.basis import Basis as _Basis, ExplicitBasis as _ExplicitBasis, DirectSumBasis as _DirectSumBasis +from pygsti.tools import sdptools as _sdps +from pygsti.baseobjs import basis as _pgb +from pygsti.baseobjs.basis import ( + Basis as _Basis, + BuiltinBasis as _BuiltinBasis, + DirectSumBasis as _DirectSumBasis, + TensorProdBasis as _TensorProdBasis +) from pygsti.baseobjs.label import Label as _Label from pygsti.baseobjs.errorgenlabel import LocalElementaryErrorgenLabel as _LocalElementaryErrorgenLabel from pygsti.tools.legacytools import deprecate as _deprecated_fn from pygsti import SpaceT +from typing import Union, Optional, Any +BasisLike = Union[str, _Basis] + IMAG_TOL = 1e-7 # tolerance for imaginary part being considered zero @@ -173,7 +184,7 @@ def psd_square_root(mat): return f -def frobeniusdist(a, b): +def frobeniusdist(a, b) -> _np.floating[Any]: """ Returns the frobenius distance between arrays: ||a - b||_Fro. @@ -296,25 +307,17 @@ def diamonddist(a, b, mx_basis='pp', return_x=False): Only returned if `return_x = True`. Encodes the state rho, such that `dm = trace( |(J(a)-J(b)).T * W| )`. """ - mx_basis = _bt.create_basis_for_matrix(a, mx_basis) - - # currently cvxpy is only needed for this function, so don't import until here + assert _sdps.CVXPY_ENABLED import cvxpy as _cp + assert a.shape == b.shape + mx_basis = _bt.create_basis_for_matrix(a, mx_basis) + c = b - a # _jam code below assumes *un-normalized* Jamiol-isomorphism. - # It will convert a & b to a "single-block" basis representation - # when mx_basis has multiple blocks. So after we call it, we need - # to multiply by mx dimension (`smallDim`). - JAstd = _jam.fast_jamiolkowski_iso_std(a, mx_basis) - JBstd = _jam.fast_jamiolkowski_iso_std(b, mx_basis) - dim = JAstd.shape[0] - smallDim = int(_np.sqrt(dim)) - JAstd *= smallDim - JBstd *= smallDim - assert(dim == JAstd.shape[1] == JBstd.shape[0] == JBstd.shape[1]) - - J = JBstd - JAstd - prob, vars = _diamond_norm_model(dim, smallDim, J) + # 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) objective_val = -2 varvals = [_np.zeros_like(J), None, None] @@ -341,62 +344,6 @@ def diamonddist(a, b, mx_basis='pp', return_x=False): return objective_val -def _diamond_norm_model(dim, smallDim, J): - # return a model for computing the diamond norm. - # - # Uses the primal SDP from arXiv:1207.5726v2, Sec 3.2 - # - # Maximize 1/2 ( < J(phi), X > + < J(phi).dag, X.dag > ) - # Subject to [[ I otimes rho0, X ], - # [ X.dag , I otimes rho1]] >> 0 - # rho0, rho1 are density matrices - # X is linear operator - - import cvxpy as _cp - - rho0 = _cp.Variable((smallDim, smallDim), name='rho0', hermitian=True) - rho1 = _cp.Variable((smallDim, smallDim), name='rho1', hermitian=True) - X = _cp.Variable((dim, dim), name='X', complex=True) - Y = _cp.real(X) - Z = _cp.imag(X) - - K = J.real - L = J.imag - if hasattr(_cp, 'scalar_product'): - objective_expr = _cp.scalar_product(K, Y) + _cp.scalar_product(L, Z) - else: - Kf = K.flatten(order='F') - Yf = Y.flatten(order='F') - Lf = L.flatten(order='F') - Zf = Z.flatten(order='F') - objective_expr = Kf @ Yf + Lf @ Zf - - objective = _cp.Maximize(objective_expr) - - ident = _np.identity(smallDim, 'd') - kr_tau0 = _cp.kron(ident, _cp.imag(rho0)) - kr_tau1 = _cp.kron(ident, _cp.imag(rho1)) - kr_sig0 = _cp.kron(ident, _cp.real(rho0)) - kr_sig1 = _cp.kron(ident, _cp.real(rho1)) - - block_11 = _cp.bmat([[kr_sig0 , Y ], - [ Y.T , kr_sig1]]) - block_21 = _cp.bmat([[kr_tau0 , Z ], - [ -Z.T , kr_tau1]]) - block_12 = block_21.T - mat_joint = _cp.bmat([[block_11, block_12], - [block_21, block_11]]) - constraints = [ - mat_joint >> 0, - rho0 >> 0, - rho1 >> 0, - _cp.trace(rho0) == 1., - _cp.trace(rho1) == 1. - ] - prob = _cp.Problem(objective, constraints) - return prob, [X, rho0, rho1] - - def jtracedist(a, b, mx_basis='pp'): # Jamiolkowski trace distance: Tr(|J(a)-J(b)|) """ Compute the Jamiolkowski trace distance between operation matrices. @@ -430,7 +377,7 @@ def jtracedist(a, b, mx_basis='pp'): # Jamiolkowski trace distance: Tr(|J(a)-J return tracedist(JA, JB) -def entanglement_fidelity(a, b, mx_basis='pp', is_tp=None, is_unitary=None): +def entanglement_fidelity(a, b, mx_basis: BasisLike='pp', is_tp=None, is_unitary=None): """ Returns the "entanglement" process fidelity between gate matrices. @@ -519,6 +466,33 @@ def is_tp_fn(x): return fidelity(JA, JB) +def tensorized_with_eye(op: _np.ndarray, op_basis: _Basis, ten_basis: Optional[_Basis]=None, std_basis: Optional[_Basis]=None, ten_std_basis: Optional[_Basis]=None): + """ + `op` is a linear operator on a Hilbert-Schmidt space, S, represented as a matrix in the `op_basis` basis. + + Returns a matrix representing a linear operator *** AND *** a basis in which to interpret that matrix. + The linear operator itself is the tensor product operator `op`⨂I_S, where I_S is the identity on S. + The basis is either ten_basis (if provided) or _TensorProdBasis((op_basis, op_basis)). + + Notes + ----- + We accept std_basis and ten_std_basis for efficiency reasons. + * If std_basis is provided, it must be equivalent to `_BuiltinBasis('std', op_basis.size)`. + * If ten_std_basis is provided, it must be equivalent to `_TensorProdBasis((std_basis, std_basis))`. + """ + if ten_basis is None: + ten_basis = _TensorProdBasis((op_basis, op_basis)) + if std_basis is None: + std_basis = _BuiltinBasis('std', op_basis.size) + if ten_std_basis is None: + ten_std_basis = _TensorProdBasis((std_basis, std_basis)) + op_std = _bt.change_basis(op, op_basis, std_basis) + eye = _np.eye(op_basis.size) + ten_op_std = _np.kron(op_std, eye) + ten_op = _bt.change_basis(ten_op_std, ten_std_basis, ten_basis) + return ten_op, ten_basis + + def average_gate_fidelity(a, b, mx_basis='pp', is_tp=None, is_unitary=None): """ Computes the average gate fidelity (AGF) between two gates. @@ -571,7 +545,7 @@ def average_gate_fidelity(a, b, mx_basis='pp', is_tp=None, is_unitary=None): d = int(round(_np.sqrt(a.shape[0]))) PF = entanglement_fidelity(a, b, mx_basis, is_tp, is_unitary) AGF = (d * PF + 1) / (1 + d) - return float(AGF) + return AGF def average_gate_infidelity(a, b, mx_basis='pp', is_tp=None, is_unitary=None): @@ -619,7 +593,7 @@ def average_gate_infidelity(a, b, mx_basis='pp', is_tp=None, is_unitary=None): return 1 - average_gate_fidelity(a, b, mx_basis, is_tp, is_unitary) -def entanglement_infidelity(a, b, mx_basis='pp', is_tp=None, is_unitary=None): +def entanglement_infidelity(a, b, mx_basis: BasisLike = 'pp', is_tp=None, is_unitary=None): """ Returns the entanglement infidelity (EI) between gate matrices. @@ -661,7 +635,7 @@ def entanglement_infidelity(a, b, mx_basis='pp', is_tp=None, is_unitary=None): EI : float The EI of a to b. """ - return 1 - float(entanglement_fidelity(a, b, mx_basis, is_tp, is_unitary)) + return 1 - entanglement_fidelity(a, b, mx_basis, is_tp, is_unitary) def generator_infidelity(a, b, mx_basis = 'pp'): """ @@ -1096,9 +1070,9 @@ def instrument_diamonddist(a, b, mx_basis): #decompose operation matrix into axis of rotation, etc def decompose_gate_matrix(operation_mx): """ - Decompse a gate matrix into fixed points, axes of rotation, angles of rotation, and decay rates. + Decompose a gate matrix into fixed points, axes of rotation, angles of rotation, and decay rates. - This funtion computes how the action of a operation matrix can be is + This function computes how the action of a operation matrix can be is decomposed into fixed points, axes of rotation, angles of rotation, and decays. Also determines whether a gate appears to be valid and/or unitary. @@ -1757,7 +1731,7 @@ def elementary_errorgens_dual(dim, typ, basis): def extract_elementary_errorgen_coefficients(errorgen, elementary_errorgen_labels, elementary_errorgen_basis='PP', errorgen_basis='pp', return_projected_errorgen=False): """ - Extract a dictionary of elemenary error generator coefficients and rates fromt he specified dense error generator + Extract a dictionary of elemenary error generator coefficients and rates from the specified dense error generator matrix. Parameters @@ -1875,7 +1849,7 @@ def project_errorgen(errorgen, elementary_errorgen_type, elementary_errorgen_bas projection values themseves. return_scale_fctr : bool, optional - If True, also return the scaling factor that was used to multply the + If True, also return the scaling factor that was used to multiply the projections onto *normalized* error generators to get the returned values. @@ -1889,7 +1863,7 @@ def project_errorgen(errorgen, elementary_errorgen_type, elementary_errorgen_bas Only returned when `return_generators == True`. An array of shape (#basis-els,op_dim,op_dim) such that `generators[i]` is the generator corresponding to the i-th basis element. Note - that these matricies are in the *std* (matrix unit) basis. + that these matrices are in the *std* (matrix unit) basis. """ if isinstance(errorgen_basis, _Basis): @@ -2206,7 +2180,7 @@ def rotation_gate_mx(r, mx_basis="gm"): Parameters ---------- r : tuple - A tuple of coeffiecients, one per non-identity + A tuple of coefficients, one per non-identity Pauli-product basis element mx_basis : {'std', 'gm', 'pp', 'qt'} or Basis object @@ -2360,7 +2334,7 @@ def project_model(model, target_model, lnd_error_gen_cp, targetOp, basis, gen_type) NpDict['LND'] += HBlk.block_data.size + otherBlk.block_data.size - #Collect and return requrested results: + #Collect and return requested results: ret_gs = [gsDict[p] for p in projectiontypes] ret_Nps = [NpDict[p] for p in projectiontypes] return ret_gs, ret_Nps @@ -2704,7 +2678,7 @@ def is_valid_lindblad_paramtype(typ): Parameters ---------- typ : str - A paramterization type. + A parameterization type. Returns ------- diff --git a/pygsti/tools/sdptools.py b/pygsti/tools/sdptools.py new file mode 100644 index 000000000..071195a20 --- /dev/null +++ b/pygsti/tools/sdptools.py @@ -0,0 +1,233 @@ +""" +Functions for constructing semidefinite programming models +""" +#*************************************************************************************************** +# Copyright 2024 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +# Under the terms of Contract DE-NA0003525 with NTESS, the U.S. Government retains certain rights +# in this software. +# Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except +# in compliance with the License. You may obtain a copy of the License at +# http://www.apache.org/licenses/LICENSE-2.0 or in the LICENSE file in the root pyGSTi directory. +#*************************************************************************************************** + +from __future__ import annotations + +import numpy as np + +from typing import Union, List, Tuple, TYPE_CHECKING +if TYPE_CHECKING: + import cvxpy as cp + ExpressionLike = Union[cp.Expression, np.ndarray] + +from pygsti.baseobjs import Basis +from pygsti.tools.basistools import stdmx_to_vec +from pygsti.tools.jamiolkowski import jamiolkowski_iso + +BasisLike = Union[str, Basis] + +try: + import cvxpy as cp + CVXPY_ENABLED = True + mosek_warning_pattern = ".*Incorrect array format causing data to be copied*" + import warnings + warnings.filterwarnings('ignore', mosek_warning_pattern) +except: + CVXPY_ENABLED = False + + +def diamond_norm_model_jamiolkowski(J: ExpressionLike) -> tuple[cp.Problem, List[cp.Variable]]: + # return a model for computing the diamond norm. + # + # Uses the primal SDP from arXiv:1207.5726v2, Sec 3.2 + # + # Throughout comments in this function, "A.dag" is the + # Hermitian adjoint (complex conjugate transpose). + # + # Maximize 1/2 ( < J, X > + < J.dag, X.dag > ) + # Subject to [[ I otimes rho0, X ], + # [ X.dag , I otimes rho1]] >> 0 + # rho0, rho1 are density matrices + # X is linear operator + # + # ".dag" returns the adjoint. + dim = J.shape[0] + smallDim = int(np.sqrt(dim)) + assert dim == smallDim**2 + + rho0 = cp.Variable((smallDim, smallDim), name='rho0', hermitian=True) + rho1 = cp.Variable((smallDim, smallDim), name='rho1', hermitian=True) + X = cp.Variable((dim, dim), name='X', complex=True) + Y = cp.real(X) + Z = cp.imag(X) + # = J.dag.ravel() @ X.ravel() + # = J.ravel() @ X.dag.ravel() = conj() + # + # ---> real() = 1/2 ( + ) + # ---> can skip the factor 1/2 if we just form real() directly. + # + + K = J.real + L = J.imag + if hasattr(cp, 'scalar_product'): + objective_expr = cp.scalar_product(K, Y) + cp.scalar_product(L, Z) + else: + Kf = K.flatten(order='F') + Yf = Y.flatten(order='F') + Lf = L.flatten(order='F') + Zf = Z.flatten(order='F') + objective_expr = Kf @ Yf + Lf @ Zf + + objective = cp.Maximize(objective_expr) + + ident = np.identity(smallDim, 'd') + kr_tau0 = cp.kron(ident, cp.imag(rho0)) + kr_tau1 = cp.kron(ident, cp.imag(rho1)) + kr_sig0 = cp.kron(ident, cp.real(rho0)) + kr_sig1 = cp.kron(ident, cp.real(rho1)) + + block_11 = cp.bmat([[kr_sig0 , Y ], + [ Y.T , kr_sig1]]) + block_21 = cp.bmat([[kr_tau0 , Z ], + [ -Z.T , kr_tau1]]) + block_12 = block_21.T + mat_joint = cp.bmat([[block_11, block_12], + [block_21, block_11]]) + constraints = [ + mat_joint >> 0, + rho0 >> 0, + rho1 >> 0, + cp.trace(rho0) == 1., + cp.trace(rho1) == 1. + ] + prob = cp.Problem(objective, constraints) + return prob, [X, rho0, rho1] + + +def diamond_norm_canon(arg : cp.Expression, basis) -> Tuple[cp.Expression, List[cp.Constraint]]: + """ + This more or less implements canonicalization of the nonlinear expression + \\|arg\\|_{\\diamond} into CVXPY Constraints and a representation of its epigraph. + The canonicalization isn't quite "complete" in CVXPY's usual sense, which would + require that the epigraph is affine and that no structured variables (like + Hermitian matrices) are used. + """ + constraints = [] + d = arg.shape[0] + small_d = int(np.sqrt(d)) + assert d == small_d**2 + assert arg.shape == (d, d) + Jarg = jamiolkowski_iso(arg, basis, basis, normalized=False) + Y0 = cp.Variable(shape=(d, d), hermitian=True) + Y1 = cp.Variable(shape=(d, d), hermitian=True) + bmat = cp.bmat([ + [ Y0 , -Jarg], + [-Jarg.T.conj(), Y1 ] + ]) + constraints.append(bmat >> 0) + TrX_Y0 = cp.partial_trace(Y0, [small_d, small_d], 0) + TrX_Y1 = cp.partial_trace(Y1, [small_d, small_d], 0) + expr0 = cp.lambda_max(TrX_Y0) + expr1 = cp.lambda_max(TrX_Y1) + epi = (expr0 + expr1)/2 + return epi, constraints + + +def cptp_superop_variable(purestate_dim: int, basis: BasisLike) -> Tuple[cp.Expression, List[cp.Constraint]]: + d = purestate_dim ** 2 + basis = Basis.cast(basis, d) + constraints = [] + if basis.first_element_is_identity: + toprow = np.zeros((1,d)) + toprow[0,0] = 1 + X_free = cp.Variable((d-1, d)) + X = cp.vstack((toprow, X_free)) + else: + X = cp.Variable((d, d)) + matI = np.eye(purestate_dim) + vecI = stdmx_to_vec(matI, basis) + constraints.append(X.T @ vecI == vecI) + """ + Let X be the process matrix for a gate "G". We have + tr(G(rho)) = < I, G(rho) > + = < vec(I), vec(G(rho)) > + = < vec(I), X @ vec(rho) > + = < X.T @ vec(I), vec(rho) >. + Therefore tr(G(rho)) = tr(rho) for all rho iff X.T @ vec(I) == vec(I). + """ + J = jamiolkowski_iso(X, basis, basis, normalized=True) + constraints.append(J >> 0) + return X, constraints + + +def diamond_distance_projection_model(superop: np.ndarray, basis: Basis, leakfree: bool=False, seepfree: bool=False, n_leak: int=0, cptp: bool=True, subspace_diamond: bool=False): + assert CVXPY_ENABLED + dim_mixed = superop.shape[0] + dim_pure = int(np.sqrt(dim_mixed)) + assert dim_pure**2 == dim_mixed + constraints = [] + if cptp: + proj_superop, cons = cptp_superop_variable(dim_pure, basis) + constraints.extend(cons) + else: + proj_superop = cp.Variable((dim_mixed, dim_mixed)) + diamondnorm_arg = superop - proj_superop + if (leakfree or seepfree or subspace_diamond): + assert n_leak == 1 + from pygsti.tools.leakage import leading_dxd_submatrix_basis_vectors + dim_pure_compsub = dim_pure - n_leak + U = leading_dxd_submatrix_basis_vectors(dim_pure_compsub, dim_pure, basis) + P = U @ U.T.conj() + I = np.eye(dim_mixed) + if leakfree: + constraints.append( (I - P) @ proj_superop @ U == 0 ) + if seepfree: + constraints.append( U.T @ proj_superop @ (I - P) == 0 ) + if subspace_diamond: + diamondnorm_arg = diamondnorm_arg @ P + expr, cons = diamond_norm_canon(diamondnorm_arg, basis) + objective = cp.Minimize(expr / 2) + # ^ We define the diamond distance between two channels as + # 1/2 the diamond norm of their difference. + constraints.extend(cons) + problem = cp.Problem(objective, constraints) + viable_solvers = [solver for solver in ['MOSEK', 'CLARABEL', 'CVXOPT'] if solver in cp.installed_solvers()] + return problem, proj_superop, viable_solvers + + +def root_fidelity_canon(sigma: cp.Expression, rho: cp.Expression) -> Tuple[cp.Expression, List[cp.Constraint]]: + """ + pyGSTi defines fidelity as + + F(sigma, rho) = tr([sigma^{1/2} rho sigma^{1/2}]^{1/2})^2. + + Others (including Neilson and Chuang, Sect. 9.2.2) define it without the + square on the trace. We'll call the unsquared version the *root fidelity,* + and denote it by + + \\sqrt{F}(sigma, rho) = (F(sigma, rho))^{1/2}. + + The root fidelity is jointly concave (Neilson and Chuang, Exercise 9.19). + In fact, it admits the following semidefinite programming characterization + + \\sqrt{F}(sigma, rho) = Maximize real(tr(X)) + s.t. [[sigma, X],[X.T.conj(), rho]] >> 0 + + -- see Section 7.1.3 of Killoran's PhD thesis, "Entanglement quantification + and quantum benchmarking of optical communication devices." + + This function returns a pair (expr, constraints) where expr is the hypograph + variable for \\sqrt{F}(sigma, rho) and constraints is a list of CVXPY Constraint + objects used in the semidefinite representation of the hypograph. + """ + t = cp.Variable() + d = sigma.shape[0] + X = cp.Variable(shape=(d, d), complex=True) + bmat = cp.hermitian_wrap(cp.bmat([ + [ sigma, X ], + [ X.T.conj(), rho ] + ])) + constraints = [ + bmat >> 0, + cp.trace(cp.real(X)) >= t + ] + return t, constraints diff --git a/test/unit/algorithms/test_gaugeopt_correctness.py b/test/unit/algorithms/test_gaugeopt_correctness.py new file mode 100644 index 000000000..af4b39c70 --- /dev/null +++ b/test/unit/algorithms/test_gaugeopt_correctness.py @@ -0,0 +1,223 @@ + +from pygsti.models.gaugegroup import TPGaugeGroup, TPGaugeGroupElement, \ + UnitaryGaugeGroup, UnitaryGaugeGroupElement,\ + FullGaugeGroup, FullGaugeGroupElement +from ..util import BaseCase + +import time +from collections import OrderedDict +import numpy as np +from pygsti.modelpacks import smq1Q_XYI +import scipy.linalg as la +import pygsti.tools.optools as pgo +import pygsti.baseobjs.label as pgl +import pygsti.baseobjs.basisconstructors as pgbc +import pygsti.algorithms.gaugeopt as gop + + + +def gate_metrics_dict(model, target): + metrics = {'infids': OrderedDict(), 'frodists': OrderedDict(), 'tracedists': OrderedDict()} + for lbl in model.operations.keys(): + model_gate = model[lbl] + target_gate = target[lbl] + metrics['infids'][lbl] = pgo.entanglement_infidelity(model_gate, target_gate) + metrics['frodists'][lbl] = pgo.frobeniusdist(model_gate, target_gate) + metrics['tracedists'][lbl] = pgo.tracedist(model_gate, target_gate) + return metrics + + +def check_gate_metrics_are_nontrivial(metrics, tol): + """ + "metrics" is a dict of the kind produced by gate_metrics_dict(model, target). + + This function checks if all values in that dict are strictly greater than tol, + with the exception of those corresponding to the idle gate. + + This makes sure that gauge optimization for "model" matching "target" is a + nontrivial problem. + """ + idle_label = pgl.Label(tuple()) + for lbl, val in metrics['infids'].items(): + if lbl == idle_label: + continue + assert val > tol, f"{val} is at most {tol}, failure for gate {lbl} w.r.t. infidelity." + for lbl, val in metrics['frodists'].items(): + if lbl == idle_label: + continue + assert val > tol, f"{val} is at most {tol}, failure for gate {lbl} w.r.t. frobenius distance." + for lbl, val in metrics['tracedists'].items(): + if lbl == idle_label: + continue + assert val > tol, f"{val} is at most {tol}, failure for gate {lbl} w.r.r. trace distance." + return + + +def check_gate_metrics_near_zero(metrics, tol): + """ + "metrics" is a dict of the kind produced by gate_metrics_dict(model, target). + + This function should be called to check if gauge optimization for "model" + matching "target" succeeded. + """ + for lbl, val in metrics['infids'].items(): + assert val <= tol, f"{val} exceeds {tol}, failure for gate {lbl} w.r.t infidelity." + for lbl, val in metrics['frodists'].items(): + assert val <= tol, f"{val} exceeds {tol}, failure for gate {lbl} w.r.t. frobenius distance." + for lbl, val in metrics['tracedists'].items(): + assert val <= tol, f"{val} exceeds {tol}, failure for gate {lbl} w.r.t. trace distance." + return + + + +class sm1QXYI_FindPerfectGauge: + + @property + def default_gauge_group(self): + """ + Instance of an object that subclasses GaugeGroup. + The state space for this gauge group is set based + on that of "self.target". + """ + return self._default_gauge_group + + @property + def gauge_grp_el_class(self): + """ + Handle for some subclass of GaugeGroupElement. + Must be compatible with self.default_gauge_group. + """ + return self._gauge_grp_el_class + + @property + def gauge_space_name(self) -> str: + """ + Only used for printing purposes in "self._main_tester". + """ + return self._gauge_space_name + + def setUp(self): + np.random.seed(0) + target = smq1Q_XYI.target_model() + self.target = target.depolarize(op_noise=0.01, spam_noise=0.001) + self.model = None + self._default_gauge_group = None + self._gauge_grp_el_class = None + self._gauge_space_name = '' + self.N_REPS = 1 + raise NotImplementedError() + + def _gauge_transform_model(self, seed): + self.model = self.target.copy() + self.model.default_gauge_group = self.default_gauge_group + np.random.seed(seed) + self.U = la.expm(np.random.randn()/2 * -1j * (pgbc.sigmax + pgbc.sigmaz)/np.sqrt(2)) + self.gauge_grp_el = self.gauge_grp_el_class(pgo.unitary_to_pauligate(self.U)) + self.model.transform_inplace(self.gauge_grp_el) + self.metrics_before = gate_metrics_dict(self.model, self.target) + check_gate_metrics_are_nontrivial(self.metrics_before, tol=1e-2) + return + + def _main_tester(self, method, seed, test_tol, alg_tol, gop_objective): + assert gop_objective in {'frobenius', 'frobenius_squared', 'tracedist', 'fidelity'} + self._gauge_transform_model(seed) + tic = time.time() + newmodel = gop.gaugeopt_to_target(self.model, self.target, method=method, tol=alg_tol, spam_metric=gop_objective, gates_metric=gop_objective) + toc = time.time() + dt = toc - tic + metrics_after = gate_metrics_dict(newmodel, self.target) + check_gate_metrics_near_zero(metrics_after, test_tol) + return dt + + def test_lbfgs_frodist(self): + self.setUp() + times = [] + for seed in range(self.N_REPS): + dt = self._main_tester('L-BFGS-B', seed, test_tol=1e-6, alg_tol=1e-8, gop_objective='frobenius') + times.append(dt) + print(f'GaugeOpt over {self.gauge_space_name} w.r.t. Frobenius dist, using L-BFGS: {times}.') + return + + def test_lbfgs_tracedist(self): + self.setUp() + times = [] + for seed in range(self.N_REPS): + dt = self._main_tester('L-BFGS-B', seed, test_tol=1e-6, alg_tol=1e-8, gop_objective='tracedist') + times.append(dt) + print(f'GaugeOpt over {self.gauge_space_name} w.r.t. trace dist, using L-BFGS: {times}.') + return + + def test_ls(self): + self.setUp() + times = [] + for seed in range(self.N_REPS): + dt = self._main_tester('ls', seed, test_tol=1e-6, alg_tol=1e-15, gop_objective='frobenius') + times.append(dt) + print(f'GaugeOpt over {self.gauge_space_name} w.r.t. Frobenius dist, using LS : {times}.') + return + + pass + + +""" +This class' test for minimizing trace distance with L-BFGS takes excessively long (about 5 seconds). +""" +class FindPerfectGauge_TPGroup_Tester(BaseCase, sm1QXYI_FindPerfectGauge): + + def setUp(self): + np.random.seed(0) + target = smq1Q_XYI.target_model() + self.target = target.depolarize(op_noise=0.01, spam_noise=0.001) + self.model = None + self._default_gauge_group = TPGaugeGroup(self.target.state_space, self.target.basis) + self._gauge_grp_el_class = TPGaugeGroupElement + self._gauge_space_name = 'TPGaugeGroup' + self.N_REPS = 1 + return + + def gauge_transform_model(self, seed): + sm1QXYI_FindPerfectGauge.gauge_transform_model(self, seed) + # ^ That sets self.model + return + + +class FindPerfectGauge_UnitaryGroup_Tester(BaseCase, sm1QXYI_FindPerfectGauge): + + def setUp(self): + np.random.seed(0) + target = smq1Q_XYI.target_model() + self.target = target.depolarize(op_noise=0.01, spam_noise=0.001) + self.model = None + self._default_gauge_group = UnitaryGaugeGroup(self.target.state_space, self.target.basis) + self._gauge_grp_el_class = UnitaryGaugeGroupElement + self._gauge_space_name = 'UnitaryGaugeGroup' + self.N_REPS = 3 + # ^ This is fast, so we can afford more tests. + return + + def gauge_transform_model(self, seed): + sm1QXYI_FindPerfectGauge.gauge_transform_model(self, seed) + # ^ That sets self.model + return + + +""" +This class' test for minimizing trace distance with L-BFGS takes excessively long (about 5 seconds). +""" +class FindPerfectGauge_FullGroup_Tester(BaseCase, sm1QXYI_FindPerfectGauge): + + def setUp(self): + np.random.seed(0) + target = smq1Q_XYI.target_model() + self.target = target.depolarize(op_noise=0.01, spam_noise=0.001) + self.model = None + self._default_gauge_group = FullGaugeGroup(self.target.state_space, self.target.basis) + self._gauge_grp_el_class = FullGaugeGroupElement + self._gauge_space_name = 'FullGaugeGroup' + self.N_REPS = 1 + return + + def gauge_transform_model(self, seed): + sm1QXYI_FindPerfectGauge.gauge_transform_model(self, seed) + # ^ That sets self.model + return diff --git a/test/unit/algorithms/test_gaugeopt.py b/test/unit/algorithms/test_gaugeopt_noerrors.py similarity index 100% rename from test/unit/algorithms/test_gaugeopt.py rename to test/unit/algorithms/test_gaugeopt_noerrors.py diff --git a/test/unit/tools/test_leakage_gst_pipeline.py b/test/unit/tools/test_leakage_gst_pipeline.py new file mode 100644 index 000000000..967d0e394 --- /dev/null +++ b/test/unit/tools/test_leakage_gst_pipeline.py @@ -0,0 +1,28 @@ + + +from pygsti.modelpacks import smq1Q_XY +from pygsti.tools.leakage import leaky_qubit_model_from_pspec, construct_leakage_report +from pygsti.data import simulate_data +from pygsti.protocols import StandardGST, ProtocolData +import unittest + + +class TestLeakageGSTPipeline(unittest.TestCase): + """ + This is a system-integration smoke test for our leakage modeling abilities. + """ + + def test_smoke(self): + # This is adapted from the Leakage-automagic ipython notebook. + mp = smq1Q_XY + ed = mp.create_gst_experiment_design(max_max_length=32) + tm3 = leaky_qubit_model_from_pspec(mp.processor_spec(), mx_basis='l2p1') + ds = simulate_data(tm3, ed.all_circuits_needing_data, num_samples=1000, seed=1997) + gst = StandardGST( modes=('CPTPLND',), target_model=tm3, verbosity=2) + pd = ProtocolData(ed, ds) + res = gst.run(pd) + _, updated_res = construct_leakage_report(res, title='easy leakage analysis!') + est = updated_res.estimates['CPTPLND'] + assert 'LAGO' in est.models + # ^ That's the leakage-aware version of 'stdgaugeopt' + return diff --git a/test/unit/tools/test_optools.py b/test/unit/tools/test_optools.py index af0a17d49..63743c2b1 100644 --- a/test/unit/tools/test_optools.py +++ b/test/unit/tools/test_optools.py @@ -4,12 +4,15 @@ import sys import numpy as np import scipy +import scipy.linalg as la from pygsti.baseobjs.basis import Basis from pygsti.baseobjs.errorgenlabel import LocalElementaryErrorgenLabel as LEEL import pygsti.tools.basistools as bt import pygsti.tools.lindbladtools as lt import pygsti.tools.optools as ot +import pygsti.tools.sdptools as sdps +import pygsti.tools.leakage as pgleak from pygsti.modelmembers.operations.lindbladcoefficients import LindbladCoefficientBlock from pygsti.modelpacks.legacy import std2Q_XXYYII from ..util import BaseCase, needs_cvxpy @@ -377,18 +380,62 @@ def setUp(self): [-0.35432747-0.27939404j, -0.02266757+0.71502652j, -0.27452307+0.07511567j, 0.35432747+0.27939404j], [ 0.71538573+0.j, 0.2680266 +0.36300238j, 0.2680266 -0.36300238j, 0.28461427+0.j]]) + def test_frobenius_distance(self): + self.assertAlmostEqual(ot.frobeniusdist(self.A, self.A), 0.0) + self.assertAlmostEqual(ot.frobeniusdist(self.A, self.B), 0.6204836823) + + self.assertAlmostEqual(ot.frobeniusdist_squared(self.A, self.A), 0.0) + self.assertAlmostEqual(ot.frobeniusdist_squared(self.A, self.B), 0.385) + def test_jtrace_distance(self): val = ot.jtracedist(self.A_TP, self.A_TP, mx_basis="pp") self.assertAlmostEqual(val, 0.0) val = ot.jtracedist(self.A_TP, self.B_unitary, mx_basis="pp") self.assertGreaterEqual(val, 0.5) + def validate_primal_diamondnorm(self, objective_val, operator, rho0, rho1, places=4): + d = rho0.shape[0] + assert rho0.shape == (d, d) + assert rho1.shape == (d, d) + assert operator.shape == (d**2, d**2) + left = np.kron(np.eye(d), la.sqrtm(rho0)) + right = np.kron(np.eye(d), la.sqrtm(rho1)) + Jop = ot._jam.jamiolkowski_iso(operator, 'pp', 'std', normalized=False) + arg = left @ Jop @ right + s = la.svdvals(arg) + expected_val = np.sum(s) + self.assertAlmostEqual(objective_val, expected_val, places=places) + return + @needs_cvxpy def test_diamond_distance(self): + mosek_warning_pattern = ".*Incorrect array format causing data to be copied*" + import warnings + warnings.filterwarnings('ignore', mosek_warning_pattern) val = ot.diamonddist(self.A_TP, self.A_TP, mx_basis="pp") self.assertAlmostEqual(val, 0.0) - val = ot.diamonddist(self.A_TP, self.B_unitary, mx_basis="pp") - self.assertGreaterEqual(val, 0.7) + objective_val, modelvars = ot.diamonddist(self.A_TP, self.B_unitary, mx_basis="pp", return_x=True) + rho0, rho1 = modelvars[1:] + self.validate_primal_diamondnorm(objective_val, self.A_TP - self.B_unitary, rho0, rho1) + self.assertGreaterEqual(objective_val, 0.7) + return + + @needs_cvxpy + def test_diamond_norm_epigraph(self): + mosek_warning_pattern = ".*Incorrect array format causing data to be copied*" + import warnings + warnings.filterwarnings('ignore', mosek_warning_pattern) + val0 = ot.diamonddist(self.A_TP, self.B_unitary, mx_basis="pp") + delta0 = self.A_TP - self.B_unitary + delta1 = bt.change_basis(delta0, 'pp', 'std') + import cvxpy as cp + obj_expr, constraints = sdps.diamond_norm_canon(delta1, 'std') + objective = cp.Minimize(obj_expr) + problem = cp.Problem(objective, constraints) + val1 = problem.solve(verbose=True, solver='CLARABEL') + self.assertGreaterEqual(val0, 0.7) + self.assertAlmostEqual(val0, val1, places=4) + return def test_entanglement_fidelity(self): fidelity_TP_unitary= ot.entanglement_fidelity(self.A_TP, self.B_unitary, is_tp=True, is_unitary=True) @@ -402,6 +449,19 @@ def test_entanglement_fidelity(self): self.assertAlmostEqual(fidelity_TP_unitary_jam, expect) self.assertAlmostEqual(fidelity_TP_unitary_std, expect) + def test_subspace_entanglement_fidelity(self): + fidelity_TP_unitary= pgleak.subspace_entanglement_fidelity(self.A_TP, self.B_unitary, 'pp') + fidelity_TP_unitary_no_flag= pgleak.subspace_entanglement_fidelity(self.A_TP, self.B_unitary, 'pp') + fidelity_TP_unitary_jam= pgleak.subspace_entanglement_fidelity(self.A_TP, self.B_unitary, 'pp') + fidelity_TP_unitary_std= pgleak.subspace_entanglement_fidelity(self.A_TP_std, self.B_unitary_std, op_basis='std') + + expect = 0.4804724656092404 + self.assertAlmostEqual(fidelity_TP_unitary, expect) + self.assertAlmostEqual(fidelity_TP_unitary_no_flag, expect) + self.assertAlmostEqual(fidelity_TP_unitary_jam, expect) + self.assertAlmostEqual(fidelity_TP_unitary_std, expect) + pass + def test_fidelity_upper_bound(self): np.random.seed(0) Q = np.linalg.qr(np.random.randn(4,4) + 1j*np.random.randn(4,4))[0]