From 5fda87bc4fd940585ac2b67610e6127eca420282 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Fri, 17 Oct 2025 15:25:56 -0700 Subject: [PATCH 1/9] couple of bugfixes in leakage.py (make leaky_qubit_model_from_pspec robust to different labelings of the idenity gate, and make the returned model`s default_gauge_group have a direct sum structure; make the default leakage gaugeopt suite 2-stage, where the second stage is over the default gauge group in the target model) --- pygsti/modelmembers/operations/linearop.py | 16 +- pygsti/modelmembers/povms/effect.py | 6 +- pygsti/protocols/gst.py | 280 +++++++++------------ pygsti/tools/leakage.py | 116 +++++++-- 4 files changed, 225 insertions(+), 193 deletions(-) diff --git a/pygsti/modelmembers/operations/linearop.py b/pygsti/modelmembers/operations/linearop.py index 7a323145c..b24cc07cc 100644 --- a/pygsti/modelmembers/operations/linearop.py +++ b/pygsti/modelmembers/operations/linearop.py @@ -15,6 +15,7 @@ from pygsti.baseobjs.opcalc import bulk_eval_compact_polynomials_complex as _bulk_eval_compact_polynomials_complex from pygsti.modelmembers import modelmember as _modelmember from pygsti.tools import optools as _ot +from pygsti.tools import matrixtools as _mt from pygsti import SpaceT from typing import Any @@ -416,11 +417,14 @@ def frobeniusdist_squared(self, other_op, transform=None, inv_transform=None) -> float """ self_mx = self.to_dense("minimal") + other_mx = other_op.to_dense("minimal") if transform is not None: self_mx = self_mx @ transform + if isinstance(inv_transform, _mt.IdentityOperator): + other_mx = other_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")) + return _ot.frobeniusdist_squared(self_mx, other_mx) def frobeniusdist(self, other_op, transform=None, inv_transform=None): @@ -506,9 +510,8 @@ def jtracedist(self, other_op, transform=None, inv_transform=None): if transform is None and inv_transform is None: return _ot.jtracedist(self.to_dense("minimal"), other_op.to_dense("minimal")) else: - return _ot.jtracedist(_np.dot( - inv_transform, _np.dot(self.to_dense("minimal"), transform)), - other_op.to_dense("minimal")) + arg = inv_transform @ self.to_dense("minimal") @ transform + return _ot.jtracedist(arg, other_op.to_dense("minimal")) def diamonddist(self, other_op, transform=None, inv_transform=None): """ @@ -535,9 +538,8 @@ def diamonddist(self, other_op, transform=None, inv_transform=None): if transform is None and inv_transform is None: return _ot.diamonddist(self.to_dense("minimal"), other_op.to_dense("minimal")) else: - return _ot.diamonddist(_np.dot( - inv_transform, _np.dot(self.to_dense("minimal"), transform)), - other_op.to_dense("minimal")) + arg = arg = inv_transform @ self.to_dense("minimal") @ transform + return _ot.diamonddist(arg, other_op.to_dense("minimal")) def transform_inplace(self, s): """ diff --git a/pygsti/modelmembers/povms/effect.py b/pygsti/modelmembers/povms/effect.py index e1e0d46ee..b1aa31546 100644 --- a/pygsti/modelmembers/povms/effect.py +++ b/pygsti/modelmembers/povms/effect.py @@ -14,6 +14,7 @@ from pygsti.modelmembers import modelmember as _modelmember from pygsti.tools import optools as _ot +from pygsti.tools import matrixtools as _mt from pygsti.baseobjs.opcalc import bulk_eval_compact_polynomials_complex as _bulk_eval_compact_polynomials_complex from typing import Any @@ -142,9 +143,12 @@ def frobeniusdist_squared(self, other_spam_vec, transform=None, inv_transform=No float """ vec = self.to_dense() + other_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()) + if isinstance(inv_transform, _mt.IdentityOperator): + other_vec = transform.T @ other_vec + return _ot.frobeniusdist_squared(vec, other_vec) def residuals(self, other_spam_vec, transform=None, inv_transform=None): """ diff --git a/pygsti/protocols/gst.py b/pygsti/protocols/gst.py index 83680912f..9e28a0093 100644 --- a/pygsti/protocols/gst.py +++ b/pygsti/protocols/gst.py @@ -1,3 +1,4 @@ + """ GST Protocol objects """ @@ -1511,6 +1512,7 @@ def run(self, data, memlimit=None, comm=None, checkpoint=None, checkpoint_path=N return ret + class LinearGateSetTomography(_proto.Protocol): """ The linear gate set tomography protocol. @@ -1971,13 +1973,6 @@ def _load_pspec(processorspec_filename_or_obj): return processorspec_filename_or_obj -def _load_model(model_filename_or_obj): - if isinstance(model_filename_or_obj, str): - return _io.load_model(model_filename_or_obj) - else: - return model_filename_or_obj # assume a Model object - - def _load_pspec_or_model(processorspec_or_model_filename_or_obj): if isinstance(processorspec_or_model_filename_or_obj, str): # if a filename is given, just try to load a processor spec (can't load a model file yet) @@ -2010,25 +2005,6 @@ def _load_fiducials_and_germs(prep_fiducial_list_or_filename, return prep_fiducials, meas_fiducials, germs -def _load_dataset(data_filename_or_set, comm, verbosity): - """Loads a DataSet from the data_filename_or_set argument of functions in this module.""" - printer = _baseobjs.VerbosityPrinter.create_printer(verbosity, comm) - if isinstance(data_filename_or_set, str): - if comm is None or comm.Get_rank() == 0: - if _os.path.splitext(data_filename_or_set)[1] == ".pkl": - with open(data_filename_or_set, 'rb') as pklfile: - ds = _pickle.load(pklfile) - else: - ds = _io.read_dataset(data_filename_or_set, True, "aggregate", printer) - if comm is not None: comm.bcast(ds, root=0) - else: - ds = comm.bcast(None, root=0) - else: - ds = data_filename_or_set # assume a Dataset object - - return ds - - def _add_gaugeopt_and_badfit(results, estlbl, target_model, gaugeopt_suite, unreliable_ops, badfit_options, optimizer, resource_alloc, printer): tref = _time.time() @@ -2243,19 +2219,16 @@ def _add_badfit_estimates(results, base_estimate_label, badfit_options, #If this estimate is the target model then skip adding the diamond distance wildcard. if base_estimate_label != 'Target': try: - budget = _compute_wildcard_budget_1d_model(base_estimate, objfn_cache, mdc_objfn, parameters, - badfit_options, printer - 1, gaugeopt_suite) - + budget = _compute_wildcard_budget_1d_model(base_estimate, mdc_objfn, printer - 1, gaugeopt_suite) #if gaugeopt_labels is None then budget with be a PrimitiveOpsSingleScaleWildcardBudget #if it was non-empty though then we'll instead have budge as a dictionary with keys #corresponding to the elements of gaugeopt_labels. In this case let's make #base_estimate.extra_parameters['wildcard1d' + "_unmodeled_error"] a dictionary of #the serialized PrimitiveOpsSingleScaleWildcardBudget elements - if gaugeopt_suite is not None and gaugeopt_suite.gaugeopt_suite_names is not None: - gaugeopt_labels = gaugeopt_suite.gaugeopt_suite_names + if gaugeopt_suite is not None: + gaugeopt_labels = budget.keys() base_estimate.extra_parameters['wildcard1d' + "_unmodeled_error"] = {lbl: budget[lbl].to_nice_serialization() for lbl in gaugeopt_labels} - base_estimate.extra_parameters['wildcard1d' + "_unmodeled_active_constraints"] \ - = None + base_estimate.extra_parameters['wildcard1d' + "_unmodeled_active_constraints"] = None base_estimate.extra_parameters["unmodeled_error"] = {lbl: budget[lbl].to_nice_serialization() for lbl in gaugeopt_labels} base_estimate.extra_parameters["unmodeled_active_constraints"] = None @@ -2316,7 +2289,7 @@ def _add_badfit_estimates(results, base_estimate_label, badfit_options, gauge_opt_params.copy(), go_gs_final, gokey, comm, printer - 1) -def _compute_wildcard_budget_1d_model(estimate, objfn_cache, mdc_objfn, parameters, badfit_options, verbosity, gaugeopt_suite=None): +def _compute_wildcard_budget_1d_model(estimate, mdc_objfn, verbosity, gaugeopt_suite): """ Create a wildcard budget for a model estimate. This version of the function produces a wildcard estimate using the model introduced in https://doi.org/10.1038/s41534-023-00764-y. @@ -2326,20 +2299,12 @@ def _compute_wildcard_budget_1d_model(estimate, objfn_cache, mdc_objfn, paramete estimate : Estimate The estimate object containing the model and data to be used. - objfn_cache : ObjectiveFunctionCache - A cache for objective function evaluations. - mdc_objfn : ModelDatasetCircuitsStore An object that stores the model, dataset, and circuits to be used in the computation. parameters : dict Various parameters of the estimate at hand. - badfit_options : GSTBadFitOptions, optional - Options specifying what post-processing actions should be performed when - a fit is unsatisfactory. Contains detailed parameters for wildcard budget - creation. - verbosity : int, optional Level of detail printed to stdout. @@ -2351,12 +2316,10 @@ def _compute_wildcard_budget_1d_model(estimate, objfn_cache, mdc_objfn, paramete Returns ------- PrimitiveOpsWildcardBudget or dict - The computed wildcard budget. If gauge optimization is applied, a dictionary of - budgets keyed by gauge optimization labels is returned. + A dictionary of budgets keyed by gauge optimization labels. """ printer = _baseobjs.VerbosityPrinter.create_printer(verbosity, mdc_objfn.resource_alloc) - badfit_options = GSTBadFitOptions.cast(badfit_options) model = mdc_objfn.model ds = mdc_objfn.dataset global_circuits_to_use = mdc_objfn.global_circuits @@ -2378,46 +2341,54 @@ def _compute_wildcard_budget_1d_model(estimate, objfn_cache, mdc_objfn, paramete two_dlogl_threshold = _chi2.ppf(1 - percentile, max(ds_dof - nparams, 1)) redbox_threshold = _chi2.ppf(1 - percentile / nboxes, 1) - ref, reference_name = _compute_1d_reference_values_and_name(estimate, badfit_options, gaugeopt_suite) - - if gaugeopt_suite is None or gaugeopt_suite.gaugeopt_suite_names is None: - gaugeopt_labels = None - primitive_ops = list(ref.keys()) - if sum([v**2 for v in ref.values()]) < 1e-4: - _warnings.warn("Reference values for 1D wildcard budget are all near-zero!" - " This usually indicates an incorrect target model and will likely cause problems computing alpha.") - - else: - gaugeopt_labels = gaugeopt_suite.gaugeopt_suite_names - primitive_ops = list(ref[list(gaugeopt_labels)[0]].keys()) - if sum([v**2 for v in ref[list(gaugeopt_labels)[0]].values()]) < 1e-4: - _warnings.warn("Reference values for 1D wildcard budget are all near-zero!" - " This usually indicates an incorrect target model and will likely cause problems computing alpha.") - - if gaugeopt_labels is None: - wcm = _wild.PrimitiveOpsSingleScaleWildcardBudget(primitive_ops, [ref[k] for k in primitive_ops], - reference_name=reference_name) - _opt.wildcardopt.optimize_wildcard_bisect_alpha(wcm, mdc_objfn, two_dlogl_threshold, redbox_threshold, printer, + target = gaugeopt_suite.gaugeopt_target + if target is None: + target = estimate.models['target'] + + # Getting the dict of gauge-optimized models is very annoying. + if gaugeopt_suite.gaugeopt_argument_dicts is None: + gaugeopt_suite.gaugeopt_argument_dicts = dict() + # ^ We make that assignment because _compute_1d_reference_values_and_name + # shouldn't have to include reassign. + gop_dicts = gaugeopt_suite.gaugeopt_argument_dicts + gop_names = gaugeopt_suite.gaugeopt_suite_names + if gop_names is None: + gop_names = [gn for gn in gop_dicts if gn in estimate.models] + + goped_models = {lbl: estimate.models[lbl] for lbl in gop_names} + if len(goped_models) == 0: + final_iter_model = estimate.models['final iteration estimate'] + if isinstance(final_iter_model, _ExplicitOpModel): + goped_model = _alg.gaugeopt_to_target(final_iter_model, target) + goped_models['simple gaugeopt'] = goped_model + else: + goped_models['no gaugeopt'] = final_iter_model + + from pygsti.baseobjs.label import Label + ref : dict[str,dict[Label, float]] = _compute_1d_reference_values_and_name(target, goped_models, gaugeopt_suite) + reference_name = 'diamond distance' + + gaugeopt_labels = list(ref.keys()) + primitive_ops = list(ref[gaugeopt_labels[0]].keys()) + if sum([v**2 for v in ref[gaugeopt_labels[0]].values()]) < 1e-4: + _warnings.warn("Reference values for 1D wildcard budget are all near-zero!" + " This usually indicates an incorrect target model and will likely cause problems computing alpha.") + + wcm = {lbl:_wild.PrimitiveOpsSingleScaleWildcardBudget(primitive_ops, [ref[lbl][k] for k in primitive_ops], + reference_name=reference_name) for lbl in gaugeopt_labels} + for budget_to_optimize in wcm.values(): + _opt.wildcardopt.optimize_wildcard_bisect_alpha(budget_to_optimize, mdc_objfn, two_dlogl_threshold, redbox_threshold, printer, guess=0.1, tol=1e-3) # results in optimized wcm - else: - wcm = {lbl:_wild.PrimitiveOpsSingleScaleWildcardBudget(primitive_ops, [ref[lbl][k] for k in primitive_ops], - reference_name=reference_name) for lbl in gaugeopt_labels} - for budget_to_optimize in wcm.values(): - _opt.wildcardopt.optimize_wildcard_bisect_alpha(budget_to_optimize, mdc_objfn, two_dlogl_threshold, redbox_threshold, printer, - guess=0.1, tol=1e-3) # results in optimized wcm return wcm -def _compute_1d_reference_values_and_name(estimate, badfit_options, gaugeopt_suite=None): +def _compute_1d_reference_values_and_name(target_model, gopped_models, gaugeopt_suite): """ - Compute the reference values and name for the 1D wildcard budget. + Compute the reference values the 1D wildcard budget. Parameters ---------- - estimate : Estimate - The estimate object containing the models to be used. - badfit_options : GSTBadFitOptions Options specifying what post-processing actions should be performed when a fit is unsatisfactory. Contains detailed parameters for wildcard budget @@ -2429,101 +2400,75 @@ def _compute_1d_reference_values_and_name(estimate, badfit_options, gaugeopt_sui Returns ------- - dict or dict of dicts - The computed reference values for the 1D wildcard budget. If gauge optimization is applied, - a dictionary of reference values keyed by gauge optimization labels is returned. - - str - The name of the reference metric used. + dict of dicts + A dictionary of 1D wildcard budget reference values, keyed by gauge optimization labels. """ - if badfit_options.wildcard1d_reference == 'diamond distance': - if gaugeopt_suite is None or gaugeopt_suite.gaugeopt_suite_names is None: - final_model = estimate.models['final iteration estimate'] - target_model = estimate.models['target'] - - if isinstance(final_model, _ExplicitOpModel): - gaugeopt_model = _alg.gaugeopt_to_target(final_model, target_model) - operations_dict = gaugeopt_model.operations - targetops_dict = target_model.operations - preps_dict = gaugeopt_model.preps - targetpreps_dict = target_model.preps - povmops_dict = gaugeopt_model.povms - insts_dict = gaugeopt_model.instruments - targetinsts_dict = target_model.instruments - + assert isinstance(gopped_models, dict) + assert all([isinstance(k, str) for k in gopped_models]) + + def _memberdicts(_m) -> tuple[dict,dict,dict,dict]: + if isinstance(_m, _ExplicitOpModel): + ops = _m.operations + preps = _m.preps + povms = _m.povms + insts = _m.instruments + else: + ops = _m.operation_blks['gates'] + preps = _m.prep_blks['layers'] + povms = {} # disabled until povm_diamonddist can handle non-ExplicitOpModels. + insts = _m.instrument_blks['layers'] + return ops, preps, povms, insts + + target_ops, target_preps, target_povms, target_insts = _memberdicts(target_model) + + dd = {lbl: dict() for lbl in gopped_models} + for lbl, gaugeopt_model in gopped_models.items(): + + argdicts = gaugeopt_suite.gaugeopt_argument_dicts.get(lbl, dict()) + n_leak: int = 0 + if isinstance(argdicts, list) and len(argdicts) > 0: + n_leak = argdicts[0].get('n_leak', n_leak) + + basis = gaugeopt_model.basis + udim = int(_np.round(_np.sqrt(basis.dim))) + I = _tools.matrixtools.IdentityOperator() + if n_leak == 0: + P = I + elif n_leak > 0: + U = _tools.leading_dxd_submatrix_basis_vectors(udim - n_leak, udim, basis) + P = U @ U.T.conj() + P = P.real + + ops, preps, _, insts = _memberdicts(gaugeopt_model) + mx_basis = gaugeopt_model.basis + + for key, op in ops.items(): + dd[lbl][key] : float = 0.5 * _tools.diamonddist(op.to_dense() @ P, target_ops[key].to_dense() @ P, mx_basis=mx_basis) # type: ignore + if dd[lbl][key] < 0: # indicates that diamonddist failed (cvxpy failure) + _warnings.warn(("Diamond distance failed to compute %s reference value for 1D wildcard budget!" + " Falling back to trace distance.") % str(key)) + dd[lbl][key] = _tools.jtracedist(op.to_dense(), target_ops[key].to_dense()) + + for key, op in insts.items(): + inst_dd : float = 0.5* _tools.instrument_diamonddist(op, target_insts[key], mx_basis=mx_basis) # type: ignore + if inst_dd < 0: # indicates that instrument_diamonddist failed + _warnings.warn(("Diamond distance failed to compute %s reference value for 1D wildcard budget!" + "No fallback presently available for instruments, so skipping.") % str(key)) else: - # Local/cloud noise models don't have default_gauge_group attribute and can't be gauge - # optimized - at least not easily. - gaugeopt_model = final_model - operations_dict = gaugeopt_model.operation_blks['gates'] - targetops_dict = target_model.operation_blks['gates'] - preps_dict = gaugeopt_model.prep_blks['layers'] - targetpreps_dict = target_model.prep_blks['layers'] - povmops_dict = {} # HACK - need to rewrite povm_diamonddist below to work - insts_dict = gaugeopt_model.instrument_blks['layers'] - targetinsts_dict = target_model.instrument_blks['layers'] - - dd = {} - for key, op in operations_dict.items(): - dd[key] = 0.5 * _tools.diamonddist(op.to_dense(), targetops_dict[key].to_dense()) - if dd[key] < 0: # indicates that diamonddist failed (cvxpy failure) - _warnings.warn(("Diamond distance failed to compute %s reference value for 1D wildcard budget!" - " Falling back to trace distance.") % str(key)) - dd[key] = _tools.jtracedist(op.to_dense(), targetops_dict[key].to_dense()) - - for key, op in insts_dict.items(): - inst_dd = 0.5* _tools.instrument_diamonddist(op, targetinsts_dict[key]) - if inst_dd < 0: # indicates that instrument_diamonddist failed - _warnings.warn(("Diamond distance failed to compute %s reference value for 1D wildcard budget!" - "No fallback presently available for instruments, so skipping.") % str(key)) - else: - dd[key] = inst_dd - - spamdd = {} - for key, op in preps_dict.items(): - spamdd[key] = _tools.tracedist(_tools.vec_to_stdmx(op.to_dense(), 'pp'), - _tools.vec_to_stdmx(targetpreps_dict[key].to_dense(), 'pp')) + dd[lbl][key] = inst_dd - for key in povmops_dict.keys(): - spamdd[key] = 0.5 * _tools.optools.povm_diamonddist(gaugeopt_model, target_model, key) + spamdd = {} + for key, op in preps.items(): + rho1 = _tools.vec_to_stdmx(op.to_dense(), basis=mx_basis) + rho2 = _tools.vec_to_stdmx(target_preps[key].to_dense(), basis=mx_basis) + spamdd[key] = _tools.tracedist(rho1, rho2) - dd['SPAM'] = sum(spamdd.values()) + for key in target_povms: + spamdd[key] = 0.5 * _tools.optools.povm_diamonddist(gaugeopt_model, target_model, key) # type: ignore - else: - #pull the target model from the gaugeopt suite. The GSTGaugeOptSuite doesn't currently - #support multiple targets, so this is safe for now. - target_model = gaugeopt_suite.gaugeopt_target - gaugeopt_models = [estimate.models[lbl] for lbl in gaugeopt_suite.gaugeopt_suite_names] - dd = {lbl: {} for lbl in gaugeopt_suite.gaugeopt_suite_names} - for gaugeopt_model, lbl in zip(gaugeopt_models, gaugeopt_suite.gaugeopt_suite_names): - for key, op in gaugeopt_model.operations.items(): - dd[lbl][key] = 0.5 * _tools.diamonddist(op.to_dense(), target_model.operations[key].to_dense()) - if dd[lbl][key] < 0: # indicates that diamonddist failed (cvxpy failure) - _warnings.warn(("Diamond distance failed to compute %s reference value for 1D wildcard budget!" - " Falling back to trace distance.") % str(key)) - dd[lbl][key] = _tools.jtracedist(op.to_dense(), target_model.operations[key].to_dense()) - - for key, op in gaugeopt_model.instruments.items(): - inst_dd = 0.5* _tools.instrument_diamonddist(op, target_model.instruments[key], gaugeopt_model.basis) - if inst_dd < 0: # indicates that instrument_diamonddist failed - _warnings.warn(("Diamond distance failed to compute %s reference value for 1D wildcard budget!" - "No fallback presently available for instruments, so skipping.") % str(key)) - else: - dd[lbl][key] = inst_dd - spamdd = {} - for key, op in gaugeopt_model.preps.items(): - spamdd[key] = _tools.tracedist(_tools.vec_to_stdmx(op.to_dense(), 'pp'), - _tools.vec_to_stdmx(target_model.preps[key].to_dense(), 'pp')) - - for key in gaugeopt_model.povms.keys(): - spamdd[key] = 0.5 * _tools.optools.povm_diamonddist(gaugeopt_model, target_model, key) + dd[lbl]['SPAM'] = sum(spamdd.values()) - dd[lbl]['SPAM'] = sum(spamdd.values()) - - return dd, 'diamond distance' - else: - raise ValueError("Invalid wildcard1d_reference value (%s) in bad-fit options!" - % str(badfit_options.wildcard1d_reference)) + return dd def _compute_robust_scaling(scale_typ, objfn_cache, mdc_objfn): @@ -3301,6 +3246,11 @@ def _add_param_preserving_gauge_opt(results: ModelEstimateResults, est_key: str, ggel = gop_dictorlist['_gaugeGroupEl'] model_implicit_gauge = transform_composed_model(est.models['final iteration estimate'], ggel) est.models[gop_name] = model_implicit_gauge + protocol = est.parameters['protocol'] + badfit_options = protocol.badfit_options + optimizer = protocol.optimizer + resource_alloc = _baseobjs.resourceallocation.ResourceAllocation() + _add_badfit_estimates(results, est_key, badfit_options, optimizer, resource_alloc, verbosity, gaugeopt_suite=gop_params) return diff --git a/pygsti/tools/leakage.py b/pygsti/tools/leakage.py index f748f6a2b..3494bfa17 100644 --- a/pygsti/tools/leakage.py +++ b/pygsti/tools/leakage.py @@ -14,9 +14,14 @@ import copy from pygsti.tools import optools as pgot from pygsti.tools import basistools as pgbt +from pygsti.tools import matrixtools as pgmt from pygsti.tools.basistools import stdmx_to_vec +from pygsti.baseobjs import Label from pygsti.baseobjs.basis import TensorProdBasis, Basis, BuiltinBasis import numpy as np +import scipy.linalg as la +import warnings + import warnings from typing import Union, Dict, Optional, List, Any, TYPE_CHECKING @@ -278,7 +283,10 @@ def subspace_superop_fro_dist(op_x, op_y, op_basis, n_leak=0): # MARK: model construction -def leaky_qubit_model_from_pspec(ps_2level: QubitProcessorSpec, mx_basis: Union[str, Basis]='l2p1', levels_readout_zero=(0,)) -> ExplicitOpModel: +def leaky_qubit_model_from_pspec( + ps_2level: QubitProcessorSpec, mx_basis: Union[str, Basis]='l2p1', + levels_readout_zero=(0,), default_idle_gatename: Label = Label(()) + ) -> 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: @@ -322,10 +330,14 @@ def leaky_qubit_model_from_pspec(ps_2level: QubitProcessorSpec, mx_basis: Union[ from pygsti.modelmembers.povms import UnconstrainedPOVM from pygsti.modelmembers.states import FullState assert ps_2level.num_qubits == 1 + if ps_2level.idle_gate_names == ['{idle}']: + ps_2level.rename_gate_inplace('{idle}', default_idle_gatename) if isinstance(mx_basis, str): mx_basis = BuiltinBasis(mx_basis, 9) assert isinstance(mx_basis, Basis) + + ql = ps_2level.qubit_labels[0] Us = ps_2level.gate_unitaries rho0 = np.array([[1, 0, 0], [0, 0, 0], [0, 0, 0]], complex) @@ -333,7 +345,7 @@ def leaky_qubit_model_from_pspec(ps_2level: QubitProcessorSpec, mx_basis: Union[ E0[levels_readout_zero, levels_readout_zero] = 1 E1 = np.eye(3, dtype=complex) - E0 - ss = ExplicitStateSpace([0],[3]) + ss = ExplicitStateSpace([ql],[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( @@ -348,9 +360,16 @@ def u2x2_to_9x9_superoperator(u2x2): return superop for gatename, unitary in Us.items(): - gatekey = (gatename, 0) if gatename != '{idle}' else ('Gi',0) + gatekey = gatename if isinstance(gatename, Label) else Label((gatename, ql)) tm_3level.operations[gatekey] = u2x2_to_9x9_superoperator(unitary) + from pygsti.models.gaugegroup import UnitaryGaugeGroup, DirectSumUnitaryGroup, TrivialGaugeGroup, _ExplicitStateSpace + ss_comp = _ExplicitStateSpace(ps_2level.qubit_labels, [2]) + ss_leak = _ExplicitStateSpace(['L'], [1]) + g_comp = UnitaryGaugeGroup(ss_comp, 'pp') + g_leak = TrivialGaugeGroup(ss_leak) + g_full = DirectSumUnitaryGroup((g_comp, g_leak), mx_basis) + tm_3level.default_gauge_group = g_full tm_3level.sim = 'map' # can use 'matrix', if that's preferred for whatever reason. return tm_3level @@ -359,10 +378,46 @@ def u2x2_to_9x9_superoperator(u2x2): 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. + goppparams_dicts is a list-of-dicts (LoDs) representation of a gauge optimization suite + suitable for models without leakage (e.g., a model of a 2-level system). + + This function returns a new gauge optimization suite (also in the LoDs representation) + by applying leakage-specific modifications to a deep-copy of gopparams_dicts. + + Example + ------- + Suppose we have a ModelEstimateResults object called `results` that includes a + CPTPLND estimate, and we want to update the models of that estimate to include + two types of leakage-aware gauge optimization. + + # + # Step 1: get the input to this function + # + estimate = results.estimates['CPTPLND'] + model = estimate.models['target'] + stdsuite = GSTGaugeOptSuite(gaugeopt_suite_names=('stdgaugeopt',)) + gopparams_dicts = stdsuite.to_dictionary(model)['stdgaugeopt'] + + # + # Step 2: use this function to build our GSTGaugeOptSuite. + # + s = lagoified_gopparams_dicts(gopparams_dicts) + t = lagoified_gopparams_dicts(gopparams_dicts) + t[-1]['gates_metric'] = 'fidelity' + t[-1]['spam_metric'] = 'fidelity' + specification = {'LAGO-std': s,'LAGO-custom': t} + gos = GSTGaugeOptSuite(gaugeopt_argument_dicts=specification) + + # + # Step 3: updating `estimate` requires that we modify `results`. + # + add_lago_models(results, 'CPTPLND', gos) + + After those lines execute, the `estimate.models` dict will have two new + key-value pairs, where the keys are 'LAGO-std' and 'LAGO-custom'. """ + from pygsti.models.gaugegroup import UnitaryGaugeGroup + tm = gopparams_dicts[0]['target_model'] 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: @@ -374,18 +429,36 @@ def lagoified_gopparams_dicts(gopparams_dicts: List[Dict]) -> List[Dict]: # 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['item_weights'] = {'gates': 0.0, 'spam': 1.0} + gg = UnitaryGaugeGroup(tm.basis.state_space, tm.basis) + inner_dict['gauge_group'] = gg + inner_dict['_gaugeGroupEl'] = gg.compute_element(gg.initial_params) + # ^ We start with gauge optimization over the full unitary group, minimizing + # SPAM differences between the estimate and the target on the computational + # subspace. Our last step of gauge optimization (which is after this loop) + # includes gates. + inner_dict['method'] = 'L-BFGS-B' + # ^ We need this optimizer because it doesn't require a gradient oracle. 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. + inner_dict = inner_dict.copy() + gg = tm.default_gauge_group + # ^ The most likely scenario is that gg is a DirectSumGaugeGroup, consisting + # of a UnitaryGaugeGroup and a TrivialGaugeGroup. Rather than hard-code + # that choice, we go with the default gauge group of the target model. + if gg is not None: + inner_dict['gauge_group'] = gg + inner_dict['_gaugeGroupEl'] = gg.compute_element(gg.initial_params) + inner_dict['gates_metric'] = 'frobenius squared' + inner_dict['spam_metric'] = 'frobenius squared' + inner_dict['item_weights'] = {'gates': 1.0, 'spam': 1.0} + gopparams_dicts.append(inner_dict) return gopparams_dicts @@ -411,20 +484,23 @@ def add_lago_models(results: ModelEstimateResults, est_key: Optional[str] = None """ 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. + + If no gauge optimization suite is provided, then we construct one by making + appropriate modifications to either the estimate's existing 'stdgaugeopt' suite + (if that exists) or to the 'stdgaugeopt' suite induced by the target model. """ 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): + if gos is None: + existing_est = results.estimates[est_key] + if 'stdgaugeopt' in existing_est.goparameters: + std_gos_lods = existing_est.goparameters['stdgaugeopt'] + lago_gos_lods = lagoified_gopparams_dicts(std_gos_lods) + gop_params = {'LAGO': lago_gos_lods} + else: + gop_params = std_lago_gopsuite(results.estimates[est_key].models['target']) + gos = GSTGaugeOptSuite(gaugeopt_argument_dicts=gop_params) _add_param_preserving_gauge_opt(results, est_key, gos, verbosity) elif est_key is None: for est_key in results.estimates.keys(): From 70314b0f59e5d792ef7a62a4dcc6c63621c0c6dc Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Wed, 22 Oct 2025 11:07:40 -0700 Subject: [PATCH 2/9] remove spurious second `arg = ` --- pygsti/modelmembers/operations/linearop.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pygsti/modelmembers/operations/linearop.py b/pygsti/modelmembers/operations/linearop.py index b24cc07cc..102b284bd 100644 --- a/pygsti/modelmembers/operations/linearop.py +++ b/pygsti/modelmembers/operations/linearop.py @@ -538,7 +538,7 @@ def diamonddist(self, other_op, transform=None, inv_transform=None): if transform is None and inv_transform is None: return _ot.diamonddist(self.to_dense("minimal"), other_op.to_dense("minimal")) else: - arg = arg = inv_transform @ self.to_dense("minimal") @ transform + arg = inv_transform @ self.to_dense("minimal") @ transform return _ot.diamonddist(arg, other_op.to_dense("minimal")) def transform_inplace(self, s): From 54f6b5509bab762a379bcef6ccdfa85c37203f2c Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Wed, 22 Oct 2025 11:25:40 -0700 Subject: [PATCH 3/9] input validation and doc changes for DirectSumUnitaryGroup --- pygsti/models/gaugegroup.py | 33 ++++++++++++++++++++------------- 1 file changed, 20 insertions(+), 13 deletions(-) diff --git a/pygsti/models/gaugegroup.py b/pygsti/models/gaugegroup.py index 2df052807..2c5d4a2f8 100644 --- a/pygsti/models/gaugegroup.py +++ b/pygsti/models/gaugegroup.py @@ -1129,25 +1129,32 @@ def _from_nice_serialization(cls, state: dict): # memo holds already de-seriali class DirectSumUnitaryGroup(GaugeGroup): """ - A subgroup of the unitary group, where the unitary operators in the group all have a - shared block-diagonal structure. + A subgroup of the unitary group on a Hilbert space H, where H has a direct sum structure. - Example setting where this is useful: - The system's Hilbert space is naturally expressed as a direct sum, H = U ⨁ V, - and we want gauge optimization to preserve the natural separation between U and V. + Unitaries in this group are block diagonal and preserve the direct sum structure of H. """ def __init__(self, subgroups: Tuple[Union[UnitaryGaugeGroup, TrivialGaugeGroup], ...], basis, name="Direct sum gauge group"): - self.subgroups = subgroups - if isinstance(basis, _Basis): - self.basis = basis - elif isinstance(basis, str): - udim = 0 - for sg in subgroups: - udim += sg.state_space.udim + + # Step 1. Get the size of the direct sum Hilbert space + udim = 0 + for sg in subgroups: + udim += sg.state_space.udim + + # Step 2. Infer and validate the provided basis. + if isinstance(basis, str): ss = _ExplicitStateSpace(['all'], [udim]) - self.basis = _BuiltinBasis("std", ss) + basis = _BuiltinBasis("std", ss) + if basis.dim != udim**2: + msg = """ + Inconsistency between `basis.dim` and the dimension of the direct sum Hilbert + space implied by the `state_space` members of the provided subgroups. + """ + raise ValueError(msg) + + self.basis = basis + self.subgroups = subgroups self._param_dims = _np.array([sg.num_params for sg in subgroups]) self._num_params = _np.sum(self._param_dims) super().__init__(name) From f154b93091c737efc0758b0cc4bd8f4222625b25 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Fri, 24 Oct 2025 10:53:43 -0700 Subject: [PATCH 4/9] correct type annotations --- pygsti/models/gaugegroup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pygsti/models/gaugegroup.py b/pygsti/models/gaugegroup.py index 2c5d4a2f8..8bc24c510 100644 --- a/pygsti/models/gaugegroup.py +++ b/pygsti/models/gaugegroup.py @@ -817,7 +817,7 @@ class UnitaryGaugeGroup(OpGaugeGroupWithBasis): to specifying the value of `pygsti.evotypes.Evotype.default_evotype`. """ - def __init__(self, state_space: _StateSpace, basis: Optional[Union[_Basis, str]], + def __init__(self, state_space: _StateSpace, basis: Union[_Basis, str], evotype: Optional[Union[_Evotype, str]] = 'default'): state_space = _StateSpace.cast(state_space) evotype = _Evotype.cast(str(evotype), default_prefer_dense_reps=True) # since we use deriv_wrt_params From 30d52d2271eb1ecaeaed733a3a7795bba6ba91ae Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Fri, 24 Oct 2025 10:57:11 -0700 Subject: [PATCH 5/9] add private _direct_sum_unitary_group helper function in leakage.py. Have lagoified_gopparams_dicts apply the direct sum group by default, since apparently the target model stored in list-of-dicts representation of `stdgaugeopt` suite overrides my attempt at setting the default gauge group --- pygsti/tools/leakage.py | 32 +++++++++++++++++++------------- 1 file changed, 19 insertions(+), 13 deletions(-) diff --git a/pygsti/tools/leakage.py b/pygsti/tools/leakage.py index 3494bfa17..e9144ab0c 100644 --- a/pygsti/tools/leakage.py +++ b/pygsti/tools/leakage.py @@ -281,6 +281,20 @@ def subspace_superop_fro_dist(op_x, op_y, op_basis, n_leak=0): return np.linalg.norm(diff @ P) +def _direct_sum_unitary_group(subspace_bases, full_basis): + from pygsti.models.gaugegroup import UnitaryGaugeGroup, DirectSumUnitaryGroup, TrivialGaugeGroup + subgroups : list[Union[TrivialGaugeGroup, UnitaryGaugeGroup]] = [] + for sb in subspace_bases: + udim = int(np.round(sb.dim**0.5)) + assert udim**2 == sb.dim + if udim == 1: + g = TrivialGaugeGroup(sb.state_space) + else: + g = UnitaryGaugeGroup(sb.state_space, sb) + subgroups.append(g) + g_full = DirectSumUnitaryGroup(tuple(subgroups), full_basis) + return g_full + # MARK: model construction def leaky_qubit_model_from_pspec( @@ -363,12 +377,8 @@ def u2x2_to_9x9_superoperator(u2x2): gatekey = gatename if isinstance(gatename, Label) else Label((gatename, ql)) tm_3level.operations[gatekey] = u2x2_to_9x9_superoperator(unitary) - from pygsti.models.gaugegroup import UnitaryGaugeGroup, DirectSumUnitaryGroup, TrivialGaugeGroup, _ExplicitStateSpace - ss_comp = _ExplicitStateSpace(ps_2level.qubit_labels, [2]) - ss_leak = _ExplicitStateSpace(['L'], [1]) - g_comp = UnitaryGaugeGroup(ss_comp, 'pp') - g_leak = TrivialGaugeGroup(ss_leak) - g_full = DirectSumUnitaryGroup((g_comp, g_leak), mx_basis) + subgroup_bases = [Basis.cast('pp', 4), Basis.cast('pp', 1)] + g_full = _direct_sum_unitary_group(subgroup_bases, mx_basis) tm_3level.default_gauge_group = g_full tm_3level.sim = 'map' # can use 'matrix', if that's preferred for whatever reason. return tm_3level @@ -448,13 +458,9 @@ def lagoified_gopparams_dicts(gopparams_dicts: List[Dict]) -> List[Dict]: # our default because add_lago_models uses parameterization-preserving # gauge optimization. inner_dict = inner_dict.copy() - gg = tm.default_gauge_group - # ^ The most likely scenario is that gg is a DirectSumGaugeGroup, consisting - # of a UnitaryGaugeGroup and a TrivialGaugeGroup. Rather than hard-code - # that choice, we go with the default gauge group of the target model. - if gg is not None: - inner_dict['gauge_group'] = gg - inner_dict['_gaugeGroupEl'] = gg.compute_element(gg.initial_params) + gg = _direct_sum_unitary_group([Basis.cast('pp', 4), Basis.cast('pp', 1)], tm.basis) + inner_dict['gauge_group'] = gg + inner_dict['_gaugeGroupEl'] = gg.compute_element(gg.initial_params) inner_dict['gates_metric'] = 'frobenius squared' inner_dict['spam_metric'] = 'frobenius squared' inner_dict['item_weights'] = {'gates': 1.0, 'spam': 1.0} From 0ef2bb0e2f160cc4263951aabc352d08372598f4 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Fri, 24 Oct 2025 10:59:17 -0700 Subject: [PATCH 6/9] have leakage_automagic.ipynb demonstrate ability to identify when a specific gate causes leakage --- .../Examples/Leakage-automagic.ipynb | 202 ++++++++++++++++-- 1 file changed, 190 insertions(+), 12 deletions(-) diff --git a/jupyter_notebooks/Examples/Leakage-automagic.ipynb b/jupyter_notebooks/Examples/Leakage-automagic.ipynb index c586ff789..f818fe175 100644 --- a/jupyter_notebooks/Examples/Leakage-automagic.ipynb +++ b/jupyter_notebooks/Examples/Leakage-automagic.ipynb @@ -4,12 +4,22 @@ "cell_type": "code", "execution_count": null, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n" + ] + } + ], "source": [ - "from pygsti.modelpacks import smq1Q_XY, smq1Q_ZN\n", + "from pygsti.modelpacks import smq1Q_XYI as mp\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" + "from pygsti.protocols import StandardGST, ProtocolData\n", + "import numpy as np\n", + "import scipy.linalg as la" ] }, { @@ -23,27 +33,195 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 2, "metadata": {}, "outputs": [], "source": [ - "mp = smq1Q_XY\n", + "def with_leaky_gate(m, gate_label, strength):\n", + " rng = np.random.default_rng(0)\n", + " v = np.concatenate([[0.0], rng.standard_normal(size=(2,))])\n", + " v /= la.norm(v)\n", + " H = v.reshape((-1, 1)) @ v.reshape((1, -1))\n", + " H *= strength\n", + " U = la.expm(1j*H)\n", + " m_copy = m.copy()\n", + " G_ideal = m_copy.operations[gate_label]\n", + " from pygsti.modelmembers.operations import ComposedOp, StaticUnitaryOp\n", + " m_copy.operations[gate_label] = ComposedOp([G_ideal, StaticUnitaryOp(U, basis=m.basis)])\n", + " return m_copy, v\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "-- Std Practice: Iter 1 of 1 (CPTPLND) --: \n", + " --- Iterative GST: [--------------------------------------------------] 0.0% 92 circuits ---\r" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/rjmurr/Documents/pygsti-leakage/hrl-repo/pygsti_repo/pygsti/modelmembers/povms/__init__.py:615: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.\n", + "To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.\n", + " errgen_vec = _np.linalg.lstsq(phys_directions, soln.x)[0]\n", + "/Users/rjmurr/Documents/pygsti-leakage/hrl-repo/pygsti_repo/pygsti/objectivefns/objectivefns.py:4478: RuntimeWarning: divide by zero encountered in divide\n", + " p5over_lsvec = 0.5/lsvec\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + " --- Iterative GST: [########------------------------------------------] 16.67% 168 circuits ---\r" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "\n", + "WARNING: Treating result as *converged* after maximum iterations (100) were exceeded.\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + " --- Iterative GST: [#################---------------------------------] 33.33% 285 circuits ---\r" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "\n", + "WARNING: Treating result as *converged* after maximum iterations (100) were exceeded.\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + " --- Iterative GST: [#########################-------------------------] 50.0% 448 circuits ---\r" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "\n", + "WARNING: Treating result as *converged* after maximum iterations (100) were exceeded.\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + " --- Iterative GST: [#################################-----------------] 66.67% 616 circuits ---\r" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "\n", + "WARNING: Treating result as *converged* after maximum iterations (100) were exceeded.\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + " --- Iterative GST: [##########################################--------] 83.33% 784 circuits ---\r" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "\n", + "WARNING: Treating result as *converged* after maximum iterations (100) were exceeded.\n", + "\n", + "WARNING: Treating result as *converged* after maximum iterations (100) were exceeded.\n", + "\n", + "WARNING: Treating result as *converged* after maximum iterations (100) were exceeded.\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + " --- Iterative GST: [##################################################] 100.0% 784 circuits ---\n" + ] + } + ], + "source": [ "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", + "# ^ Target model. \"Leaky\" is a bit of a misnomer here. The returned model\n", + "# is simply a qutrit lift of the qubit model; leakage erorrs in the\n", + "# qubit model can manifest as CPTP Markovian errors in the qutrit model.\n", + "dgm3, leaking_state = with_leaky_gate(tm3, ('Gxpi2', 0), strength=0.125)\n", + "# ^ Data generating model. \n", + "ds = simulate_data(dgm3, ed.all_circuits_needing_data, num_samples=1000, seed=1997)\n", + "gst = StandardGST(\n", + " modes=('CPTPLND',), target_model=tm3, verbosity=2,\n", + " badfit_options={'actions': ['wildcard1d'], 'threshold': 0.0}\n", + ")\n", "pd = ProtocolData(ed, ds)\n", "res = gst.run(pd)" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 4, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/rjmurr/Documents/pygsti-leakage/hrl-repo/pygsti_repo/pygsti/models/model.py:143: UserWarning: Model.num_modeltest_params could not obtain number of *non-gauge* parameters - using total instead\n", + " _warnings.warn((\"Model.num_modeltest_params could not obtain number of *non-gauge* parameters\"\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Running idle tomography\n", + "Computing switchable properties\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/rjmurr/Documents/pygsti-leakage/hrl-repo/pygsti_repo/pygsti/report/plothelpers.py:367: UserWarning:\n", + "\n", + "Max-model params (92) <= model params (161)! Using k == 1.\n", + "\n", + "/Users/rjmurr/Documents/pygsti-leakage/hrl-repo/pygsti_repo/pygsti/forwardsims/mapforwardsim.py:732: UserWarning:\n", + "\n", + "Generating dense process matrix representations of circuits or gates \n", + "can be inefficient and should be avoided for the purposes of forward \n", + "simulation/calculation of circuit outcome probability distributions \n", + "when using the MapForwardSimulator.\n", + "\n", + "/Users/rjmurr/Documents/pygsti-leakage/hrl-repo/pygsti_repo/pygsti/report/plothelpers.py:367: UserWarning:\n", + "\n", + "Max-model params (92) <= model params (161)! Using k == 1.\n", + "\n" + ] + } + ], "source": [ "report_dir = 'example_files/leakage-report-automagic'\n", "report_object, updated_res = construct_leakage_report(res, title='easy leakage analysis!')\n", From 681bab9e83f39cc35bc36fdd53d0797b0918ccd9 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Fri, 24 Oct 2025 11:15:51 -0700 Subject: [PATCH 7/9] notebook: reduce GST circuit depth and clear output --- .../Examples/Leakage-automagic.ipynb | 166 +----------------- 1 file changed, 7 insertions(+), 159 deletions(-) diff --git a/jupyter_notebooks/Examples/Leakage-automagic.ipynb b/jupyter_notebooks/Examples/Leakage-automagic.ipynb index f818fe175..d487f30da 100644 --- a/jupyter_notebooks/Examples/Leakage-automagic.ipynb +++ b/jupyter_notebooks/Examples/Leakage-automagic.ipynb @@ -4,15 +4,7 @@ "cell_type": "code", "execution_count": null, "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "\n" - ] - } - ], + "outputs": [], "source": [ "from pygsti.modelpacks import smq1Q_XYI as mp\n", "from pygsti.tools.leakage import leaky_qubit_model_from_pspec, construct_leakage_report\n", @@ -33,7 +25,7 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -55,122 +47,16 @@ "cell_type": "code", "execution_count": null, "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "-- Std Practice: Iter 1 of 1 (CPTPLND) --: \n", - " --- Iterative GST: [--------------------------------------------------] 0.0% 92 circuits ---\r" - ] - }, - { - "name": "stderr", - "output_type": "stream", - "text": [ - "/Users/rjmurr/Documents/pygsti-leakage/hrl-repo/pygsti_repo/pygsti/modelmembers/povms/__init__.py:615: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.\n", - "To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.\n", - " errgen_vec = _np.linalg.lstsq(phys_directions, soln.x)[0]\n", - "/Users/rjmurr/Documents/pygsti-leakage/hrl-repo/pygsti_repo/pygsti/objectivefns/objectivefns.py:4478: RuntimeWarning: divide by zero encountered in divide\n", - " p5over_lsvec = 0.5/lsvec\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - " --- Iterative GST: [########------------------------------------------] 16.67% 168 circuits ---\r" - ] - }, - { - "name": "stderr", - "output_type": "stream", - "text": [ - "\n", - "WARNING: Treating result as *converged* after maximum iterations (100) were exceeded.\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - " --- Iterative GST: [#################---------------------------------] 33.33% 285 circuits ---\r" - ] - }, - { - "name": "stderr", - "output_type": "stream", - "text": [ - "\n", - "WARNING: Treating result as *converged* after maximum iterations (100) were exceeded.\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - " --- Iterative GST: [#########################-------------------------] 50.0% 448 circuits ---\r" - ] - }, - { - "name": "stderr", - "output_type": "stream", - "text": [ - "\n", - "WARNING: Treating result as *converged* after maximum iterations (100) were exceeded.\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - " --- Iterative GST: [#################################-----------------] 66.67% 616 circuits ---\r" - ] - }, - { - "name": "stderr", - "output_type": "stream", - "text": [ - "\n", - "WARNING: Treating result as *converged* after maximum iterations (100) were exceeded.\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - " --- Iterative GST: [##########################################--------] 83.33% 784 circuits ---\r" - ] - }, - { - "name": "stderr", - "output_type": "stream", - "text": [ - "\n", - "WARNING: Treating result as *converged* after maximum iterations (100) were exceeded.\n", - "\n", - "WARNING: Treating result as *converged* after maximum iterations (100) were exceeded.\n", - "\n", - "WARNING: Treating result as *converged* after maximum iterations (100) were exceeded.\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - " --- Iterative GST: [##################################################] 100.0% 784 circuits ---\n" - ] - } - ], + "outputs": [], "source": [ - "ed = mp.create_gst_experiment_design(max_max_length=32)\n", + "ed = mp.create_gst_experiment_design(max_max_length=16)\n", "tm3 = leaky_qubit_model_from_pspec(mp.processor_spec(), mx_basis='l2p1')\n", "# ^ Target model. \"Leaky\" is a bit of a misnomer here. The returned model\n", "# is simply a qutrit lift of the qubit model; leakage erorrs in the\n", "# qubit model can manifest as CPTP Markovian errors in the qutrit model.\n", "dgm3, leaking_state = with_leaky_gate(tm3, ('Gxpi2', 0), strength=0.125)\n", "# ^ Data generating model. \n", - "ds = simulate_data(dgm3, ed.all_circuits_needing_data, num_samples=1000, seed=1997)\n", + "ds = simulate_data(dgm3, ed.all_circuits_needing_data, num_samples=10_000, seed=1997)\n", "gst = StandardGST(\n", " modes=('CPTPLND',), target_model=tm3, verbosity=2,\n", " badfit_options={'actions': ['wildcard1d'], 'threshold': 0.0}\n", @@ -181,47 +67,9 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": null, "metadata": {}, - "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "/Users/rjmurr/Documents/pygsti-leakage/hrl-repo/pygsti_repo/pygsti/models/model.py:143: UserWarning: Model.num_modeltest_params could not obtain number of *non-gauge* parameters - using total instead\n", - " _warnings.warn((\"Model.num_modeltest_params could not obtain number of *non-gauge* parameters\"\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Running idle tomography\n", - "Computing switchable properties\n" - ] - }, - { - "name": "stderr", - "output_type": "stream", - "text": [ - "/Users/rjmurr/Documents/pygsti-leakage/hrl-repo/pygsti_repo/pygsti/report/plothelpers.py:367: UserWarning:\n", - "\n", - "Max-model params (92) <= model params (161)! Using k == 1.\n", - "\n", - "/Users/rjmurr/Documents/pygsti-leakage/hrl-repo/pygsti_repo/pygsti/forwardsims/mapforwardsim.py:732: UserWarning:\n", - "\n", - "Generating dense process matrix representations of circuits or gates \n", - "can be inefficient and should be avoided for the purposes of forward \n", - "simulation/calculation of circuit outcome probability distributions \n", - "when using the MapForwardSimulator.\n", - "\n", - "/Users/rjmurr/Documents/pygsti-leakage/hrl-repo/pygsti_repo/pygsti/report/plothelpers.py:367: UserWarning:\n", - "\n", - "Max-model params (92) <= model params (161)! Using k == 1.\n", - "\n" - ] - } - ], + "outputs": [], "source": [ "report_dir = 'example_files/leakage-report-automagic'\n", "report_object, updated_res = construct_leakage_report(res, title='easy leakage analysis!')\n", From 7018fa698c78d041a008818af889a0a1c7c08151 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Fri, 24 Oct 2025 12:26:49 -0700 Subject: [PATCH 8/9] steal infrastructure for setting objective function tolerances from another branch --- pygsti/objectivefns/objectivefns.py | 98 +++++++++++++++++------------ 1 file changed, 59 insertions(+), 39 deletions(-) diff --git a/pygsti/objectivefns/objectivefns.py b/pygsti/objectivefns/objectivefns.py index 81e40b1ca..d69141b44 100644 --- a/pygsti/objectivefns/objectivefns.py +++ b/pygsti/objectivefns/objectivefns.py @@ -10,6 +10,8 @@ # http://www.apache.org/licenses/LICENSE-2.0 or in the LICENSE file in the root pyGSTi directory. #*************************************************************************************************** +from typing import Optional + import itertools as _itertools import sys as _sys import time as _time @@ -29,6 +31,11 @@ from pygsti.models.model import OpModel as _OpModel from pygsti import SpaceT + +DEFAULT_MIN_PROB_CLIP = 1e-4 +DEFAULT_RADIUS = 1e-4 + + def _objfn(objfn_cls, model, dataset, circuits=None, regularization=None, penalties=None, op_label_aliases=None, comm=None, mem_limit=None, method_names=None, array_types=None, @@ -186,24 +193,25 @@ def create_from(cls, objective='logl', freq_weighted_chi2=False): """ if objective == "chi2": if freq_weighted_chi2: + mpc = RawFreqWeightedChi2Function.DEFAULT_MIN_PROB_CLIP_FOR_WEIGHTING builder = FreqWeightedChi2Function.builder( name='fwchi2', description="Freq-weighted sum of Chi^2", - regularization={'min_freq_clip_for_weighting': 1e-4}) + regularization={'min_freq_clip_for_weighting': mpc}) + else: + mpc = RawChi2Function.DEFAULT_MIN_PROB_CLIP_FOR_WEIGHTING builder = Chi2Function.builder( name='chi2', description="Sum of Chi^2", - regularization={'min_prob_clip_for_weighting': 1e-4}) + regularization={'min_prob_clip_for_weighting': mpc}) elif objective == "logl": builder = PoissonPicDeltaLogLFunction.builder( name='dlogl', description="2*Delta(log(L))", - regularization={'min_prob_clip': 1e-4, - 'radius': 1e-4}, - penalties={'cptp_penalty_factor': 0, - 'spam_penalty_factor': 0}) + regularization={'min_prob_clip': DEFAULT_MIN_PROB_CLIP, 'radius': DEFAULT_RADIUS}, + penalties={'cptp_penalty_factor': 0, 'spam_penalty_factor': 0}) elif objective == "tvd": builder = TVDFunction.builder( @@ -1660,6 +1668,9 @@ class RawChi2Function(RawObjectiveFunction): verbosity : int, optional Level of detail to print to stdout. """ + + DEFAULT_MIN_PROB_CLIP_FOR_WEIGHTING = DEFAULT_MIN_PROB_CLIP + def __init__(self, regularization=None, resource_alloc=None, name="chi2", description="Sum of Chi^2", verbosity=0): super().__init__(regularization, resource_alloc, name, description, verbosity) @@ -1678,7 +1689,7 @@ def chi2k_distributed_qty(self, objective_function_value): """ return objective_function_value - def set_regularization(self, min_prob_clip_for_weighting=1e-4): + def set_regularization(self, min_prob_clip_for_weighting: Optional[float]=None): """ Set regularization values. @@ -1692,7 +1703,11 @@ def set_regularization(self, min_prob_clip_for_weighting=1e-4): ------- None """ - self.min_prob_clip_for_weighting = min_prob_clip_for_weighting + if min_prob_clip_for_weighting is not None: + self.min_prob_clip_for_weighting = min_prob_clip_for_weighting + else: + self.min_prob_clip_for_weighting = RawChi2Function.DEFAULT_MIN_PROB_CLIP_FOR_WEIGHTING + return def lsvec(self, probs, counts, total_counts, freqs, intermediates=None): """ @@ -2286,7 +2301,6 @@ def _zero_freq_dterms_relaxed(self, total_counts, probs): # negative lsvec possible class RawFreqWeightedChi2Function(RawChi2Function): - """ The function `N(p-f)^2 / f` @@ -2307,6 +2321,9 @@ class RawFreqWeightedChi2Function(RawChi2Function): verbosity : int, optional Level of detail to print to stdout. """ + + DEFAULT_MIN_PROB_CLIP_FOR_WEIGHTING = RawChi2Function.DEFAULT_MIN_PROB_CLIP_FOR_WEIGHTING + def __init__(self, regularization=None, resource_alloc=None, name="fwchi2", description="Sum of freq-weighted Chi^2", verbosity=0): super().__init__(regularization, resource_alloc, name, description, verbosity) @@ -2326,7 +2343,7 @@ def chi2k_distributed_qty(self, objective_function_value): """ return objective_function_value # default is to assume the value *is* chi2_k distributed - def set_regularization(self, min_freq_clip_for_weighting=1e-4): + def set_regularization(self, min_freq_clip_for_weighting: Optional[float]=None): """ Set regularization values. @@ -2340,7 +2357,11 @@ def set_regularization(self, min_freq_clip_for_weighting=1e-4): ------- None """ - self.min_freq_clip_for_weighting = min_freq_clip_for_weighting + if min_freq_clip_for_weighting is not None: + self.min_freq_clip_for_weighting = min_freq_clip_for_weighting + else: + self.min_freq_clip_for_weighting = RawFreqWeightedChi2Function.DEFAULT_MIN_PROB_CLIP_FOR_WEIGHTING + return def _weights(self, p, f, total_counts): #Note: this could be computed once and cached? @@ -2744,8 +2765,10 @@ def chi2k_distributed_qty(self, objective_function_value): """ return 2 * objective_function_value # 2 * deltaLogL is what is chi2_k distributed - def set_regularization(self, min_prob_clip=1e-4, pfratio_stitchpt=None, pfratio_derivpt=None, - radius=1e-4, fmin=None): + def set_regularization(self, + min_prob_clip=DEFAULT_MIN_PROB_CLIP, pfratio_stitchpt=None, pfratio_derivpt=None, + radius=DEFAULT_RADIUS, fmin=None + ): """ Set regularization values. @@ -3144,7 +3167,7 @@ def chi2k_distributed_qty(self, objective_function_value): """ return 2 * objective_function_value # 2 * deltaLogL is what is chi2_k distributed - def set_regularization(self, min_prob_clip=1e-4, pfratio_stitchpt=None, pfratio_derivpt=None): + def set_regularization(self, min_prob_clip=DEFAULT_MIN_PROB_CLIP, pfratio_stitchpt=None, pfratio_derivpt=None): """ Set regularization values. @@ -4374,15 +4397,12 @@ def terms(self, paramvec=None, oob_check=False, profiler_str="TERMS OBJECTIVE"): self.model.from_vector(paramvec) terms = self.obj.view() - # Whether this rank is the "leader" of all the processors accessing the same shared self.jac and self.probs mem. - # Only leader processors should modify the contents of the shared memory, so we only apply operations *once* - # `unit_ralloc` is the group of all the procs targeting same destination into self.obj unit_ralloc = self.layout.resource_alloc('atom-processing') shared_mem_leader = unit_ralloc.is_host_leader with self.resource_alloc.temporarily_track_memory(self.nelements): # 'e' (terms) - self.model.sim.bulk_fill_probs(self.probs, self.layout) # syncs shared mem - self._clip_probs() # clips self.probs in place w/shared mem sync + self.model.sim.bulk_fill_probs(self.probs, self.layout) + self._clip_probs() if oob_check: # Only used for termgap cases if not self.model.sim.bulk_test_if_paths_are_sufficient(self.layout, self.probs, verbosity=1): @@ -4415,15 +4435,11 @@ def lsvec(self, paramvec=None, oob_check=False, raw_objfn_lsvec_signs=True): in {bad_locs}. """ raise RuntimeError(msg) - unit_ralloc = self.layout.resource_alloc('atom-processing') - shared_mem_leader = unit_ralloc.is_host_leader - if shared_mem_leader: - lsvec **= 0.5 + lsvec **= 0.5 if raw_objfn_lsvec_signs: - if shared_mem_leader: + if self.layout.resource_alloc('atom-processing').is_host_leader: raw_lsvec = self.raw_objfn.lsvec(self.probs, self.counts, self.total_counts, self.freqs) lsvec[:self.nelements][raw_lsvec < 0] *= -1 - unit_ralloc.host_comm_barrier() return lsvec def dterms(self, paramvec=None): @@ -5420,7 +5436,7 @@ def _array_types_for_method(cls, method_name, fsim): if method_name == 'dlsvec': return fsim._array_types_for_method('bulk_fill_timedep_dchi2') return super()._array_types_for_method(method_name, fsim) - def set_regularization(self, min_prob_clip_for_weighting=1e-4, radius=1e-4): + def set_regularization(self, min_prob_clip_for_weighting: Optional[float]=None, radius: float = DEFAULT_RADIUS): """ Set regularization values. @@ -5437,6 +5453,9 @@ def set_regularization(self, min_prob_clip_for_weighting=1e-4, radius=1e-4): ------- None """ + if min_prob_clip_for_weighting is None: + mpc = getattr(self.raw_objfn, 'DEFAULT_MIN_PROB_CLIP_FOR_WEIGHTING', DEFAULT_MIN_PROB_CLIP) + min_prob_clip_for_weighting = mpc self.min_prob_clip_for_weighting = min_prob_clip_for_weighting self.radius = radius # parameterizes "roundness" of f == 0 terms @@ -5579,7 +5598,7 @@ def _array_types_for_method(cls, method_name, fsim): if method_name == 'dlsvec': return fsim._array_types_for_method('bulk_fill_timedep_dloglpp') return super()._array_types_for_method(method_name, fsim) - def set_regularization(self, min_prob_clip=1e-4, radius=1e-4): + def set_regularization(self, min_prob_clip=DEFAULT_MIN_PROB_CLIP, radius=DEFAULT_RADIUS): """ Set regularization values. @@ -5835,7 +5854,7 @@ def _paramvec_norm_penalty(reg_factor, paramvec): return out -def _cptp_penalty_jac_fill(cp_penalty_vec_grad_to_fill, mdl, prefactor, op_basis, wrt_slice): +def _cptp_penalty_jac_fill(cp_penalty_vec_grad_to_fill, mdl, prefactor, op_basis, wrt_slice, error_tol=1e-4): """ Helper function - jacobian of CPTP penalty (sum of tracenorms of gates) Returns a (real) array of shape (len(mdl.operations), n_params). @@ -5849,7 +5868,7 @@ def _cptp_penalty_jac_fill(cp_penalty_vec_grad_to_fill, mdl, prefactor, op_basis #get sgn(chi-matrix) == d(|chi|_Tr)/dchi in std basis # so sgnchi == d(|chi_std|_Tr)/dchi_std chi = _tools.fast_jamiolkowski_iso_std(gate, op_basis) - assert(_np.linalg.norm(chi - chi.T.conjugate()) < 1e-4), \ + assert(_np.linalg.norm(chi - chi.T.conjugate()) < error_tol), \ "chi should be Hermitian!" # Alt#1 way to compute sgnchi (evals) - works equally well to svd below @@ -5862,7 +5881,7 @@ def _cptp_penalty_jac_fill(cp_penalty_vec_grad_to_fill, mdl, prefactor, op_basis # _spl.sqrtm(_np.matrix(_np.dot(chi.T.conjugate(),chi))))) sgnchi = _tools.matrix_sign(chi) - assert(_np.linalg.norm(sgnchi - sgnchi.T.conjugate()) < 1e-4), \ + assert(_np.linalg.norm(sgnchi - sgnchi.T.conjugate()) < error_tol), \ "sgnchi should be Hermitian!" # get d(gate)/dp in op_basis [shape == (nP,dim,dim)] @@ -5879,7 +5898,8 @@ def _cptp_penalty_jac_fill(cp_penalty_vec_grad_to_fill, mdl, prefactor, op_basis m_dgate_dp_std = _np.empty((nparams, mdl.dim, mdl.dim), 'complex') for p in range(nparams): # p indexes param m_dgate_dp_std[p] = _tools.fast_jamiolkowski_iso_std(dgate_dp[p], op_basis) # now "M(dGdp_std)" - assert(_np.linalg.norm(m_dgate_dp_std[p] - m_dgate_dp_std[p].T.conjugate()) < 1e-8) # check hermitian + assert(_np.linalg.norm(m_dgate_dp_std[p] - m_dgate_dp_std[p].T.conjugate()) < error_tol**2) # check hermitian + # Riley note: not clear to me why we're using the square of error_tol in the above. m_dgate_dp_std = _np.conjugate(m_dgate_dp_std) # so element-wise multiply # of complex number (einsum below) results in separately adding @@ -5890,7 +5910,7 @@ def _cptp_penalty_jac_fill(cp_penalty_vec_grad_to_fill, mdl, prefactor, op_basis #v = _np.einsum("ij,aij->a",sgnchi,MdGdp_std) v = _np.tensordot(sgnchi, m_dgate_dp_std, ((0, 1), (1, 2))) v *= prefactor * (0.5 / _np.sqrt(_tools.tracenorm(chi))) # add 0.5/|chi|_Tr factor - assert(_np.linalg.norm(v.imag) < 1e-4) + assert(_np.linalg.norm(v.imag) < error_tol) cp_penalty_vec_grad_to_fill[i, :] = 0.0 gate_gpinds_subset, within_wrtslice, within_gpinds = _slct.intersect_within(wrt_slice, gate.gpindices) @@ -5901,7 +5921,7 @@ def _cptp_penalty_jac_fill(cp_penalty_vec_grad_to_fill, mdl, prefactor, op_basis return len(mdl.operations) # the number of leading-dim indicies we filled in -def _spam_penalty_jac_fill(spam_penalty_vec_grad_to_fill, mdl, prefactor, op_basis, wrt_slice): +def _spam_penalty_jac_fill(spam_penalty_vec_grad_to_fill, mdl, prefactor, op_basis, wrt_slice, error_tol=1e-4): """ Helper function - jacobian of CPTP penalty (sum of tracenorms of gates) Returns a (real) array of shape ( _spam_penalty_size(mdl), n_params). @@ -5924,11 +5944,11 @@ def _spam_penalty_jac_fill(spam_penalty_vec_grad_to_fill, mdl, prefactor, op_bas #get sgn(denMx) == d(|denMx|_Tr)/d(denMx) in std basis # dmDim = denMx.shape[0] denmx = _tools.vec_to_stdmx(prepvec, op_basis) - assert(_np.linalg.norm(denmx - denmx.T.conjugate()) < 1e-4), \ + assert(_np.linalg.norm(denmx - denmx.T.conjugate()) < error_tol), \ "denMx should be Hermitian!" sgndm = _tools.matrix_sign(denmx) - assert(_np.linalg.norm(sgndm - sgndm.T.conjugate()) < 1e-4), \ + assert(_np.linalg.norm(sgndm - sgndm.T.conjugate()) < error_tol), \ "sgndm should be Hermitian!" # get d(prepvec)/dp in op_basis [shape == (nP,dim)] @@ -5943,7 +5963,7 @@ def _spam_penalty_jac_fill(spam_penalty_vec_grad_to_fill, mdl, prefactor, op_bas #v = _np.einsum("ij,aij,ab->b",sgndm,ddenMxdV,dVdp) v = _np.tensordot(_np.tensordot(sgndm, ddenmx_dv, ((0, 1), (1, 2))), dv_dp, (0, 0)) v *= prefactor * (0.5 / _np.sqrt(_tools.tracenorm(denmx))) # add 0.5/|denMx|_Tr factor - assert(_np.linalg.norm(v.imag) < 1e-4) + assert(_np.linalg.norm(v.imag) < error_tol) spam_penalty_vec_grad_to_fill[i, :] = 0.0 gpinds_subset, within_wrtslice, within_gpinds = _slct.intersect_within(wrt_slice, prepvec.gpindices) spam_penalty_vec_grad_to_fill[i, within_wrtslice] = v.real[within_gpinds] # slice or array index works! @@ -5960,11 +5980,11 @@ def _spam_penalty_jac_fill(spam_penalty_vec_grad_to_fill, mdl, prefactor, op_bas #get sgn(EMx) == d(|EMx|_Tr)/d(EMx) in std basis emx = _tools.vec_to_stdmx(effectvec, op_basis) # dmDim = EMx.shape[0] - assert(_np.linalg.norm(emx - emx.T.conjugate()) < 1e-4), \ + assert(_np.linalg.norm(emx - emx.T.conjugate()) < error_tol), \ "EMx should be Hermitian!" sgn_e = _tools.matrix_sign(emx) - assert(_np.linalg.norm(sgn_e - sgn_e.T.conjugate()) < 1e-4), \ + assert(_np.linalg.norm(sgn_e - sgn_e.T.conjugate()) < error_tol), \ "sgnE should be Hermitian!" # get d(prepvec)/dp in op_basis [shape == (nP,dim)] @@ -5979,7 +5999,7 @@ def _spam_penalty_jac_fill(spam_penalty_vec_grad_to_fill, mdl, prefactor, op_bas #v = _np.einsum("ij,aij,ab->b",sgnE,dEMxdV,dVdp) v = _np.tensordot(_np.tensordot(sgn_e, demx_dv, ((0, 1), (1, 2))), dv_dp, (0, 0)) v *= prefactor * (0.5 / _np.sqrt(_tools.tracenorm(emx))) # add 0.5/|EMx|_Tr factor - assert(_np.linalg.norm(v.imag) < 1e-4) + assert(_np.linalg.norm(v.imag) < error_tol) spam_penalty_vec_grad_to_fill[i, :] = 0.0 gpinds_subset, within_wrtslice, within_gpinds = _slct.intersect_within(wrt_slice, effectvec.gpindices) From 3c8ef6e927dbbfddb8a72ac01cea4bdf26fe0f80 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Fri, 24 Oct 2025 12:28:04 -0700 Subject: [PATCH 9/9] make it slightly cheaper to run the example notebook. Make a proper correctness test for leakage GST. --- .../Examples/Leakage-automagic.ipynb | 17 +++- test/unit/tools/test_leakage_gst_pipeline.py | 96 +++++++++++++++++-- 2 files changed, 101 insertions(+), 12 deletions(-) diff --git a/jupyter_notebooks/Examples/Leakage-automagic.ipynb b/jupyter_notebooks/Examples/Leakage-automagic.ipynb index d487f30da..c2e3701d8 100644 --- a/jupyter_notebooks/Examples/Leakage-automagic.ipynb +++ b/jupyter_notebooks/Examples/Leakage-automagic.ipynb @@ -49,14 +49,27 @@ "metadata": {}, "outputs": [], "source": [ - "ed = mp.create_gst_experiment_design(max_max_length=16)\n", + "ed = mp.create_gst_experiment_design(max_max_length=8)\n", + "# ^ The default max length is small so we don't have to wait as long \n", + "# for the GST fit (just for purposes of this notebook).\n", "tm3 = leaky_qubit_model_from_pspec(mp.processor_spec(), mx_basis='l2p1')\n", "# ^ Target model. \"Leaky\" is a bit of a misnomer here. The returned model\n", "# is simply a qutrit lift of the qubit model; leakage erorrs in the\n", "# qubit model can manifest as CPTP Markovian errors in the qutrit model.\n", "dgm3, leaking_state = with_leaky_gate(tm3, ('Gxpi2', 0), strength=0.125)\n", "# ^ Data generating model. \n", - "ds = simulate_data(dgm3, ed.all_circuits_needing_data, num_samples=10_000, seed=1997)\n", + "num_samples = 100_000\n", + "# ^ The number of samples is large to compensate for short circuit length.\n", + "# Feel free to change the number of samples to something more \"realistic\"\n", + "# if you'd like.\n", + "if num_samples > 10_000:\n", + " from pygsti.objectivefns import objectivefns\n", + " objectivefns.DEFAULT_MIN_PROB_CLIP = objectivefns.DEFAULT_RADIUS = 1e-12\n", + " # ^ There are numerical thresholding rules in objective function evaluation\n", + " # that lead to errors when the number of samples is extremely large.\n", + " # The lines above change those thresholding rules to be appropriate in\n", + " # the unusual setting that is this notebook.\n", + "ds = simulate_data(dgm3, ed.all_circuits_needing_data, num_samples=num_samples, seed=1997)\n", "gst = StandardGST(\n", " modes=('CPTPLND',), target_model=tm3, verbosity=2,\n", " badfit_options={'actions': ['wildcard1d'], 'threshold': 0.0}\n", diff --git a/test/unit/tools/test_leakage_gst_pipeline.py b/test/unit/tools/test_leakage_gst_pipeline.py index 967d0e394..d36e53a46 100644 --- a/test/unit/tools/test_leakage_gst_pipeline.py +++ b/test/unit/tools/test_leakage_gst_pipeline.py @@ -1,28 +1,104 @@ -from pygsti.modelpacks import smq1Q_XY +from pygsti.modelpacks import smq1Q_XYI 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 +import numpy as np +import scipy.linalg as la + + +def with_leaky_gate(m, gate_label, strength): + rng = np.random.default_rng(0) + v = np.concatenate([[0.0], rng.standard_normal(size=(2,))]) + v /= la.norm(v) + H = v.reshape((-1, 1)) @ v.reshape((1, -1)) + H *= strength + U = la.expm(1j*H) + m_copy = m.copy() + G_ideal = m_copy.operations[gate_label] + from pygsti.modelmembers.operations import ComposedOp, StaticUnitaryOp + m_copy.operations[gate_label] = ComposedOp([G_ideal, StaticUnitaryOp(U, basis=m.basis)]) + return m_copy, v class TestLeakageGSTPipeline(unittest.TestCase): """ - This is a system-integration smoke test for our leakage modeling abilities. + This is a simple system-integration test for our leakage modeling abilities; + it takes just under 1 minute to run on an Apple M2 Max machine. """ - def test_smoke(self): + def test_pipeline_1Q_XYI(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) + mp = smq1Q_XYI + ed = mp.create_gst_experiment_design(max_max_length=8) + # ^ The max length is small so we don't have to wait as long for the GST fit. + tm3 = leaky_qubit_model_from_pspec(mp.processor_spec()) + # ^ Target model. + dgm3, _ = with_leaky_gate(tm3, ('Gxpi2', 0), strength=0.125) + # ^ Data generating model. + num_samples = 100_000 + # ^ The number of samples is large to compensate for short circuit length. + from pygsti.objectivefns import objectivefns + objectivefns.DEFAULT_MIN_PROB_CLIP = objectivefns.DEFAULT_RADIUS = 1e-12 + # ^ The lines above change numerical thresholding rules in objective evaluation + # to be appropriate when the number of shots/circuit is extremely large. + ds = simulate_data(dgm3, ed.all_circuits_needing_data, num_samples=num_samples, seed=1997) + gst = StandardGST( + modes=('CPTPLND',), target_model=tm3, verbosity=2, + badfit_options={'actions': ['wildcard1d'], 'threshold': 0.0} + ) pd = ProtocolData(ed, ds) res = gst.run(pd) _, updated_res = construct_leakage_report(res, title='easy leakage analysis!') + # ^ we do that as a smoke test for construct_leakage_report and to get our hands + # on the results updated with leakage-aware gauge optimization. est = updated_res.estimates['CPTPLND'] - assert 'LAGO' in est.models - # ^ That's the leakage-aware version of 'stdgaugeopt' + + """ + Original results are shown below. We don't rely on the exact numbers here. What matters is + qualitative aspects of how Gxpi2 and Gypi2 deviate from their respective targets. Since our + data generating model only applied leakage to Gxpi2, a "good" result reports much more error + in Gxpi2 than Gypi2. (It's not clear to me why stdgaugeopt lacks wildcard error.) + + Leakage-aware guage optimization. + + | Gate | ent. infidelity | 1/2 trace dist | 1/2 diamond dist | Max TOP | Unmodeled error | + |---------|-----------------|----------------|------------------|----------|-----------------| + | [] | 0.000001 | 0.000522 | 0.000729 | 0.000384 | 0.000001 | + | Gxpi2:0 | 0.00207 | 0.045378 | 0.062144 | 0.048625 | 0.003824 | + | Gypi2:0 | 0.000188 | 0.013716 | 0.016257 | 0.010192 | 0.000152 | + + Standard gauge optimization + + | Gate | ent. infidelity | 1/2 trace dist | 1/2 diamond dist | Max TOP | + |---------|-----------------|----------------|------------------|----------| + | [] | 0.000006 | 0.002033 | 0.002866 | 0.002749 | + | Gxpi2:0 | 0.00061 | 0.024568 | 0.032364 | 0.024514 | + | Gypi2:0 | 0.000602 | 0.024526 | 0.033052 | 0.023098 | + + We'll run tests with subspace entanglement infidelity. + + * For LAGO, infidelity of Gxpi2 is 10x larger than that of Gypi2; + we'll test for a 5x difference. + + * For standard gauge optimization, Gxpi2 and Gypi2 have almost identical infidelities; + we'll test for a factor 1.1x there. + """ + from pygsti.tools.leakage import subspace_entanglement_fidelity as fidelity + + mdls = {lbl: est.models[lbl] for lbl in {'target', 'LAGO', 'stdgaugeopt'}} + assert mdls['target'].basis.name == mdls['stdgaugeopt'].basis.name == mdls['LAGO'].basis.name == 'l2p1' + + gates = dict() + for lbl, mdl in mdls.items(): + gates[lbl] = {g: mdl.operations[(f'G{g}pi2',0)].to_dense() for g in ['x', 'y']} + + infids = dict() + for lbl in ['LAGO', 'stdgaugeopt']: + infids[lbl] = {g: 1 - fidelity(gates[lbl][g], gates['target'][g], 'l2p1', n_leak=1) for g in ['x', 'y'] } + + self.assertGreater( infids['LAGO']['x'], 5.0 * infids['LAGO']['y'] ) + self.assertLess( infids['stdgaugeopt']['x'], 1.1 * infids['stdgaugeopt']['y'] ) return