From 0865b7f3537181588a0ca679083f49ce7a89607c Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Mon, 29 Sep 2025 16:03:59 -0700 Subject: [PATCH 01/18] initial simplification of wildcard budget code in gst.py --- pygsti/protocols/gst.py | 202 ++++++++++++++++------------------------ 1 file changed, 80 insertions(+), 122 deletions(-) diff --git a/pygsti/protocols/gst.py b/pygsti/protocols/gst.py index 83680912f..8558ccb8c 100644 --- a/pygsti/protocols/gst.py +++ b/pygsti/protocols/gst.py @@ -2243,9 +2243,7 @@ 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 @@ -2316,7 +2314,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=None): """ 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. @@ -2356,7 +2354,6 @@ def _compute_wildcard_budget_1d_model(estimate, objfn_cache, mdc_objfn, paramete """ 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 +2375,49 @@ 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, + # Getting the target model is a little annoying. + target = getattr(gaugeopt_suite, 'gaugeopt_target', None) + if target is None: + target = estimate.models['target'] + + # Getting the dict of gauge-optimized models is very annoying. + gop_names = getattr(gaugeopt_suite, 'gaugeopt_suite_names', None) + if gop_names is None: + gop_names = ['stdgaugeopt'] if 'stdgaugeopt' in estimate.models else [] + 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) + 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): """ - 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 +2429,59 @@ 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 - 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 + target_ops, target_preps, target_povms, target_insts = _memberdicts(target_model) - 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: {} for lbl in gopped_models} + for lbl, gaugeopt_model in gopped_models.items(): - for key in povmops_dict.keys(): - spamdd[key] = 0.5 * _tools.optools.povm_diamonddist(gaugeopt_model, target_model, key) + ops, preps, _, insts = _memberdicts(gaugeopt_model) - dd['SPAM'] = sum(spamdd.values()) + for key, op in ops.items(): + dd[lbl][key] = 0.5 * _tools.diamonddist(op.to_dense(), target_ops[key].to_dense()) # 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 = 0.5* _tools.instrument_diamonddist(op, target_insts[key]) # 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: + dd[lbl][key] = inst_dd - 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')) + spamdd = {} + for key, op in preps.items(): + rho1 = _tools.vec_to_stdmx(op.to_dense(), 'pp') + rho2 = _tools.vec_to_stdmx(target_preps[key].to_dense(), 'pp') + spamdd[key] = _tools.tracedist(rho1, rho2) - for key in gaugeopt_model.povms.keys(): - spamdd[key] = 0.5 * _tools.optools.povm_diamonddist(gaugeopt_model, target_model, key) + for key in target_povms: + spamdd[key] = 0.5 * _tools.optools.povm_diamonddist(gaugeopt_model, target_model, key) # type: ignore - 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): From 58eff0787b8a1bc39048bd60e16b2b71e7449cdc Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Mon, 29 Sep 2025 16:12:00 -0700 Subject: [PATCH 02/18] remove unused code --- pygsti/protocols/gst.py | 38 ++------------------------------------ 1 file changed, 2 insertions(+), 36 deletions(-) diff --git a/pygsti/protocols/gst.py b/pygsti/protocols/gst.py index 8558ccb8c..71c6c588f 100644 --- a/pygsti/protocols/gst.py +++ b/pygsti/protocols/gst.py @@ -1511,6 +1511,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 +1972,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 +2004,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() @@ -2324,20 +2299,12 @@ def _compute_wildcard_budget_1d_model(estimate, mdc_objfn, verbosity, gaugeopt_s 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. @@ -2349,8 +2316,7 @@ def _compute_wildcard_budget_1d_model(estimate, mdc_objfn, verbosity, gaugeopt_s 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) From 4a8b7fb741df0425cc4da7ba8de71e94034c1cd6 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Tue, 30 Sep 2025 16:04:03 -0700 Subject: [PATCH 03/18] Clean up style in reportables.py. Move some leakage metric definitions from reportables.py into pygsti/tools/leakage.py. Bugfix for qubit label assumption in leaky_qubit_model_from_pspec. --- pygsti/report/reportables.py | 61 +++++++++++++++++++----------------- pygsti/tools/leakage.py | 53 +++++++++++++++++++++++++++---- 2 files changed, 79 insertions(+), 35 deletions(-) diff --git a/pygsti/report/reportables.py b/pygsti/report/reportables.py index f92dc8d36..a40def411 100644 --- a/pygsti/report/reportables.py +++ b/pygsti/report/reportables.py @@ -1016,14 +1016,17 @@ 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( @@ -1037,8 +1040,10 @@ def diamonddist_to_leakfree_cptp(op, ignore, mx_basis): 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( @@ -1052,49 +1057,41 @@ def subspace_diamonddist_to_leakfree_cptp(op, ignore, mx_basis): 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 + from pygsti.tools.leakage import subspace_diamonddist as _subspace_diamonddist + return _subspace_diamonddist(op_a, op_b, basis) + 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) + from pygsti.tools.leakage import leakage_profile + evals, _ = leakage_profile(op, mx_basis) + return evals[0] -def pergate_leakrate_min(op, ignore, mx_basis): - return pergate_leakrate_reduction(op, ignore, mx_basis, min) PerGateLeakRateMax = _modf.opsfn_factory(pergate_leakrate_max) + + +def pergate_leakrate_min(op, ignore, mx_basis): + from pygsti.tools.leakage import leakage_profile + evals, _ = leakage_profile(op, mx_basis) + return evals[-1] + + 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 + from pygsti.tools.leakage import seepage_profile + seeprate, _ = seepage_profile(op, mx_basis) + return seeprate[0] + PerGateSeepRate = _modf.opsfn_factory(pergate_seeprate) @@ -1169,10 +1166,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) @@ -1201,9 +1200,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) @@ -1305,10 +1306,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) diff --git a/pygsti/tools/leakage.py b/pygsti/tools/leakage.py index f748f6a2b..676cacef0 100644 --- a/pygsti/tools/leakage.py +++ b/pygsti/tools/leakage.py @@ -144,7 +144,7 @@ def apply_tensorized_to_teststate(op_x: np.ndarray, op_y, op_basis: np.ndarray, @lru_cache -def leading_dxd_submatrix_basis_vectors(d: int, n: int, current_basis: Basis): +def leading_dxd_submatrix_basis_vectors(d: int, n: int, current_basis: Basis) -> np.ndarray: """ 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 @@ -237,7 +237,7 @@ def leading_dxd_submatrix_basis_vectors(d: int, n: int, current_basis: Basis): @set_docstring(CHOI_INDUCED_METRIC % 'entanglement fidelity') -def subspace_entanglement_fidelity(op_x, op_y, op_basis, n_leak=0): +def subspace_entanglement_fidelity(op_x: np.ndarray, op_y: np.ndarray, op_basis, n_leak=0) -> float: 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) @@ -246,7 +246,7 @@ def subspace_entanglement_fidelity(op_x, op_y, op_basis, n_leak=0): @set_docstring(CHOI_INDUCED_METRIC % 'jamiolkowski trace distance') -def subspace_jtracedist(op_x, op_y, op_basis, n_leak=0): +def subspace_jtracedist(op_x: np.ndarray, op_y: np.ndarray, op_basis, n_leak=0) -> float: 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) @@ -264,7 +264,7 @@ def subspace_jtracedist(op_x, op_y, op_basis, n_leak=0): *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): +def subspace_superop_fro_dist(op_x: np.ndarray, op_y: np.ndarray, op_basis, n_leak=0) -> float: diff = op_x - op_y if n_leak == 0: return np.linalg.norm(diff, 'fro') @@ -276,6 +276,45 @@ def subspace_superop_fro_dist(op_x, op_y, op_basis, n_leak=0): return np.linalg.norm(diff @ P) +def subspace_diamonddist(op_x: np.ndarray, op_y: np.ndarray, op_basis) -> float: + dim_mixed = op_x.shape[0] + dim_pure = int(dim_mixed**0.5) + dim_pure_compsub = dim_pure - 1 + U = leading_dxd_submatrix_basis_vectors(dim_pure_compsub, dim_pure, op_basis) + P = U @ U.T.conj() + assert np.linalg.norm(P - P.real) < 1e-10 + P = P.real + from pygsti.tools.optools import diamonddist + val : float = diamonddist(op_x @ P, op_y @ P, op_basis, return_x=False) / 2 # type: ignore + return val + + +def leakage_profile(op, mx_basis) -> tuple[np.ndarray, list[np.ndarray]]: + assert op.shape == (9, 9) + lfb = BuiltinBasis('l2p1', 9) + op_lfb = pgbt.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 = pgbt.vec_to_stdmx(leakage_effect_superket, 'pp') + leakage_rates, states = np.linalg.eigh(leakage_effect) + ind = np.argsort(leakage_rates)[::-1] + rates = leakage_rates[ind] + states = [np.concatenate([s,[0.0]]) for s in states.T[ind]] + return rates, states + + +def seepage_profile(op, mx_basis) -> tuple[np.ndarray, list[np.ndarray]]: + assert op.shape == (9, 9) + lfb = BuiltinBasis('l2p1', 9) + op_lfb = pgbt.change_basis(op, mx_basis, lfb) + elinds = lfb.elindlookup + seeprate = np.atleast_1d(op_lfb[elinds['I'], elinds['L']]) + state = [np.array([0.0, 0.0, 1.0], dtype=np.complex128)] + return seeprate, state + + + # MARK: model construction def leaky_qubit_model_from_pspec(ps_2level: QubitProcessorSpec, mx_basis: Union[str, Basis]='l2p1', levels_readout_zero=(0,)) -> ExplicitOpModel: @@ -326,6 +365,8 @@ def leaky_qubit_model_from_pspec(ps_2level: QubitProcessorSpec, mx_basis: Union[ 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 +374,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,7 +389,7 @@ 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, ql) if gatename != '{idle}' else ('Gi', ql) tm_3level.operations[gatekey] = u2x2_to_9x9_superoperator(unitary) tm_3level.sim = 'map' # can use 'matrix', if that's preferred for whatever reason. From f6f31a0836e155a9449c4b8ba1467f308f90069d Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Mon, 6 Oct 2025 13:21:56 -0700 Subject: [PATCH 04/18] fix mathematical bug in definition of subspace-restricted Frobenius norm --- pygsti/modelmembers/operations/linearop.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/pygsti/modelmembers/operations/linearop.py b/pygsti/modelmembers/operations/linearop.py index 7a323145c..a7097ba07 100644 --- a/pygsti/modelmembers/operations/linearop.py +++ b/pygsti/modelmembers/operations/linearop.py @@ -14,7 +14,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 optools as _ot, matrixtools as _mt from pygsti import SpaceT from typing import Any @@ -416,11 +416,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): From d6faeb66bc915b12bea788b67f7c7baa87ce1d45 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Mon, 6 Oct 2025 13:54:19 -0700 Subject: [PATCH 05/18] apply bugfix to last commit (needed to include POVM effects, whereas last commit just applied a fix to linear operators). Also, make it possible to use different weights on state prep vs `spam` writ-large in gauge optimization --- pygsti/modelmembers/povms/effect.py | 6 +++++- pygsti/models/explicitcalc.py | 3 ++- 2 files changed, 7 insertions(+), 2 deletions(-) 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/models/explicitcalc.py b/pygsti/models/explicitcalc.py index 6717e1f66..2c1120a0f 100644 --- a/pygsti/models/explicitcalc.py +++ b/pygsti/models/explicitcalc.py @@ -172,6 +172,7 @@ def frobeniusdist(self, other_calc, transform_mx: Union[None, _np.ndarray, Trans if item_weights is None: item_weights = {} opWeight = item_weights.get('gates', 1.0) spamWeight = item_weights.get('spam', 1.0) + prepWeight = item_weights.get('prep', spamWeight) for opLabel, gate in self.operations.items(): wt = item_weights.get(opLabel, opWeight) @@ -179,7 +180,7 @@ def frobeniusdist(self, other_calc, transform_mx: Union[None, _np.ndarray, Trans nSummands += wt * (gate.dim)**2 for lbl, rhoV in self.preps.items(): - wt = item_weights.get(lbl, spamWeight) + wt = item_weights.get(lbl, prepWeight) d += wt * rhoV.frobeniusdist_squared(other_calc.preps[lbl], P, invP) nSummands += wt * rhoV.dim From 36e1e36ba902d2af3bd3d7e19e7fe6847a3e7375 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Mon, 6 Oct 2025 14:04:52 -0700 Subject: [PATCH 06/18] make note of possible changes to eigenvalue_entanglement_infidelity --- pygsti/report/reportables.py | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/pygsti/report/reportables.py b/pygsti/report/reportables.py index a40def411..b145c5d53 100644 --- a/pygsti/report/reportables.py +++ b/pygsti/report/reportables.py @@ -1600,11 +1600,17 @@ def eigenvalue_entanglement_infidelity(a, b, mx_basis): float """ d2 = a.shape[0] - evA = _np.linalg.eigvals(a) - evB = _np.linalg.eigvals(b) - _, pairs = _tools.minweight_match(evA, evB, lambda x, y: abs(x - y), + evalsA, evecsA = _np.linalg.eig(a) + evalsB, evecsB = _np.linalg.eig(b) + """ + IDEA: check if (a,b) are normal operators. If they are, then compute + proceed by computing a Schur decomposition to get orthonormal eigenbases. + From there, match them based on a dissimilarity kernel of abs(1 - ). + That could be used + """ + _, pairs = _tools.minweight_match(evalsA, evalsB, lambda x, y: abs(x - y), return_pairs=True) # just to get pairing - mlPl = abs(_np.sum([_np.conjugate(evB[j]) * evA[i] for i, j in pairs])) + mlPl = abs(_np.sum([_np.conjugate(evalsB[j]) * evalsA[i] for i, j in pairs])) return 1.0 - mlPl / float(d2) From 57fea370c4b6adb7b1b27cf8664545b4ff31764e Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Mon, 6 Oct 2025 14:11:11 -0700 Subject: [PATCH 07/18] make add_lago_models robust to estimates that are missing `stdgaugeopt` from their gopparameters member. Update lagoified_gopparams_dict so the gauge optimization step that `stdgaugeopt` does over the unitary group is REPLACED by gauge optimization over the target models default_gauge_group (when that default exists). Update leaky_qubit_model_from_pspec to handle idle gates in a more robust way and to set the default gauge group to a direct sum unitary group that preserves separation between computational subspace and leakage space. --- pygsti/tools/leakage.py | 30 +++++++++++++++++++++++++----- 1 file changed, 25 insertions(+), 5 deletions(-) diff --git a/pygsti/tools/leakage.py b/pygsti/tools/leakage.py index 676cacef0..9bdf19b6c 100644 --- a/pygsti/tools/leakage.py +++ b/pygsti/tools/leakage.py @@ -15,6 +15,7 @@ from pygsti.tools import optools as pgot from pygsti.tools import basistools as pgbt 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 warnings @@ -317,7 +318,10 @@ def seepage_profile(op, mx_basis) -> tuple[np.ndarray, list[np.ndarray]]: # 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: @@ -361,6 +365,8 @@ 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) @@ -389,9 +395,16 @@ def u2x2_to_9x9_superoperator(u2x2): return superop for gatename, unitary in Us.items(): - gatekey = (gatename, ql) if gatename != '{idle}' else ('Gi', ql) + 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 @@ -407,6 +420,10 @@ def lagoified_gopparams_dicts(gopparams_dicts: List[Dict]) -> List[Dict]: 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: + gg = inner_dict['target_model'].default_gauge_group + if gg is not None: + inner_dict['gauge_group'] = gg + inner_dict['_gaugeGroupEl'] = gg.compute_element(gg.initial_params) 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. @@ -461,9 +478,12 @@ def add_lago_models(results: ModelEstimateResults, est_key: Optional[str] = None 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} + 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) if isinstance(est_key, str): _add_param_preserving_gauge_opt(results, est_key, gos, verbosity) From 1d5df521485447d0ddd452c2e42ef944e9f6b68f Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Mon, 6 Oct 2025 17:30:54 -0700 Subject: [PATCH 08/18] seems to fix diamond-distance wildcard with leakage models. Still need to write tests. --- pygsti/protocols/gst.py | 62 ++++++++++++++++++++++++++++++----------- 1 file changed, 46 insertions(+), 16 deletions(-) diff --git a/pygsti/protocols/gst.py b/pygsti/protocols/gst.py index 71c6c588f..3a3382283 100644 --- a/pygsti/protocols/gst.py +++ b/pygsti/protocols/gst.py @@ -2224,11 +2224,10 @@ def _add_badfit_estimates(results, base_estimate_label, badfit_options, #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 @@ -2289,7 +2288,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, mdc_objfn, 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. @@ -2341,15 +2340,20 @@ def _compute_wildcard_budget_1d_model(estimate, mdc_objfn, verbosity, gaugeopt_s two_dlogl_threshold = _chi2.ppf(1 - percentile, max(ds_dof - nparams, 1)) redbox_threshold = _chi2.ppf(1 - percentile / nboxes, 1) - # Getting the target model is a little annoying. - target = getattr(gaugeopt_suite, 'gaugeopt_target', None) + target = gaugeopt_suite.gaugeopt_target if target is None: target = estimate.models['target'] # Getting the dict of gauge-optimized models is very annoying. - gop_names = getattr(gaugeopt_suite, 'gaugeopt_suite_names', None) + 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 = ['stdgaugeopt'] if 'stdgaugeopt' in estimate.models else [] + 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'] @@ -2360,7 +2364,7 @@ def _compute_wildcard_budget_1d_model(estimate, mdc_objfn, verbosity, gaugeopt_s 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) + 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()) @@ -2378,7 +2382,7 @@ def _compute_wildcard_budget_1d_model(estimate, mdc_objfn, verbosity, gaugeopt_s return wcm -def _compute_1d_reference_values_and_name(target_model, gopped_models): +def _compute_1d_reference_values_and_name(target_model, gopped_models, gaugeopt_suite): """ Compute the reference values the 1D wildcard budget. @@ -2416,20 +2420,41 @@ def _memberdicts(_m) -> tuple[dict,dict,dict,dict]: target_ops, target_preps, target_povms, target_insts = _memberdicts(target_model) - dd = {lbl: {} for lbl in gopped_models} + 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 = 0 + if isinstance(argdicts, list) and len(argdicts) > 0: + n_leak = argdicts[-1].get('n_leak', n_leak) + + if n_leak > 0: + dim = gaugeopt_model.basis.state_space.udim + assert n_leak == 1 + assert dim == 3 + B = _tools.leading_dxd_submatrix_basis_vectors(2, 3, gaugeopt_model.basis) + P = B @ B.T.conj() + if _np.linalg.norm(P.imag) > 1e-12: + msg = f"Attempting to run leakage-aware gauge optimization with basis {gaugeopt_model.basis}\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) + P = P.real + else: + P = _tools.matrixtools.IdentityOperator() + ops, preps, _, insts = _memberdicts(gaugeopt_model) + mx_basis = gaugeopt_model.basis for key, op in ops.items(): - dd[lbl][key] = 0.5 * _tools.diamonddist(op.to_dense(), target_ops[key].to_dense()) # type: ignore + 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 = 0.5* _tools.instrument_diamonddist(op, target_insts[key]) # type: ignore + 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)) @@ -2438,8 +2463,8 @@ def _memberdicts(_m) -> tuple[dict,dict,dict,dict]: spamdd = {} for key, op in preps.items(): - rho1 = _tools.vec_to_stdmx(op.to_dense(), 'pp') - rho2 = _tools.vec_to_stdmx(target_preps[key].to_dense(), 'pp') + 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) for key in target_povms: @@ -3225,6 +3250,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 From b74c4a05f667abe9471bbc1a2eceafc0f385aadc Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Tue, 7 Oct 2025 10:36:26 -0700 Subject: [PATCH 09/18] define proper constants (either module-level or class-level) that hold the values of what are currently magic numbers hard-coded to 1e-4. --- pygsti/objectivefns/objectivefns.py | 83 +++++++++++++++++++---------- 1 file changed, 55 insertions(+), 28 deletions(-) diff --git a/pygsti/objectivefns/objectivefns.py b/pygsti/objectivefns/objectivefns.py index a88975f16..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. @@ -5413,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. @@ -5430,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 @@ -5572,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. @@ -5828,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). @@ -5842,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 @@ -5855,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)] @@ -5872,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 @@ -5883,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) @@ -5894,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). @@ -5917,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)] @@ -5936,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! @@ -5953,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)] @@ -5972,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 e0a8fbeca9c9d9090c5bb77c7be3fec54161a8bf Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Wed, 8 Oct 2025 09:54:05 -0700 Subject: [PATCH 10/18] add a second stage to leakage-aware gauge optimization --- pygsti/algorithms/gaugeopt.py | 77 +++++++++++++++++++++++++++++------ pygsti/protocols/gst.py | 2 +- pygsti/tools/leakage.py | 31 ++++++++++---- 3 files changed, 88 insertions(+), 22 deletions(-) diff --git a/pygsti/algorithms/gaugeopt.py b/pygsti/algorithms/gaugeopt.py index c332553fa..c68549e21 100644 --- a/pygsti/algorithms/gaugeopt.py +++ b/pygsti/algorithms/gaugeopt.py @@ -372,6 +372,7 @@ def _create_objective_fn(model, target_model, item_weights: Optional[dict[str,fl if item_weights is None: item_weights = {} opWeight = item_weights.get('gates', 1.0) spamWeight = item_weights.get('spam', 1.0) + prepWeight = item_weights.get('prep', spamWeight) mxBasis = model.basis #Use the target model's basis if model's is unknown @@ -623,6 +624,44 @@ def _mock_objective_fn(v): else: transform_mx_arg = None # ^ It would be equivalent to set this to a pair of IdentityOperator objects. + from pygsti.report.reportables import eigenvalue_entanglement_infidelity, vec_fidelity + + + """ + We say that a (model, target) pair admit a _perfect gauge_ if "model" has no relational errors + when compared to "target." If "model" is considered in the perfect gauge, then + + (1) its fidelities with the target gates match the eigenvalue fidelities with the target gates; + + (2) the fidelities between rho_e and {E_t : E_t in M_target} will match those between + rho_e and {E_e : E_e in M_estimate}; and + + (3) the fidelities between rho_t and {E_e : E_e in M_estimate} will match those between + rho_t and {E_t : E_t in M_target}. + + In view of this, taking fidelity as the gauge-optimization objective tries to minimize an + aggregration of any mismatch between the various fidelities above. + """ + + gate_fidelity_targets : dict[ Union[str, _baseobjs.Label], Union[float, _np.floating] ] = dict() + if gates_metric == 'fidelity': + for lbl in target_model.operations: + G_target = target_model.operations[lbl] + G_curest = model.operations[lbl] + t = 1 - eigenvalue_entanglement_infidelity(G_curest.to_dense(), G_target.to_dense(), model.basis) + gate_fidelity_targets[lbl] = min(t, 1.0) + + spam_fidelity_targets : dict[ + tuple[Union[str, _baseobjs.Label], Union[str, _baseobjs.Label]], + dict[Union[str, _baseobjs.Label], Union[float, _np.floating]] + ] = dict() + if spam_metric == 'fidelity': + for preplbl in target_model.preps: + for povmlbl in target_model.povms: + rho_curest = model.preps[preplbl].to_dense() + M_curest = model.povms[povmlbl] + t = {elbl: vec_fidelity(rho_curest, e.to_dense(), mxBasis).item() for (elbl, e) in M_curest.items() } + spam_fidelity_targets[(preplbl, povmlbl)] = t def _objective_fn(gauge_group_el, oob_check): mdl = _transform_with_oob_check(model, gauge_group_el, oob_check) @@ -656,12 +695,15 @@ def _objective_fn(gauge_group_el, oob_check): ret += val elif gates_metric == "fidelity": - # If n_leak==0, then subspace_entanglement_fidelity is just entanglement_fidelity + # Leakage-aware metrics NOT available 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 + t = gate_fidelity_targets[opLbl] + v = _tools.entanglement_fidelity(top, mop, mxBasis) + z = _np.abs(t - v) + ret += wt * z elif gates_metric == "tracedist": # If n_leak==0, then subspace_jtracedist is just jtracedist. @@ -691,17 +733,26 @@ def _objective_fn(gauge_group_el, oob_check): 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 + val = 0.0 + for preplbl in target_model.preps: + wt_prep = item_weights.get(preplbl, prepWeight) + for povmlbl in target_model.povms: + wt_povm = item_weights.get(povmlbl, spamWeight) + rho_curest = model.preps[preplbl].to_dense() + rho_target = target_model.preps[preplbl].to_dense() + M_curest = model.povms[povmlbl] + M_target = target_model.povms[povmlbl] + vs_prep = {elbl: vec_fidelity(rho_curest, e.to_dense(), mxBasis) for (elbl, e) in M_target.items() } + vs_povm = {elbl: vec_fidelity(rho_target, e.to_dense(), mxBasis) for (elbl, e) in M_curest.items() } + t_dict = spam_fidelity_targets[(preplbl, povmlbl)] + val1 = 0.0 + for lbl, f in vs_prep.items(): + val1 += _np.abs(t_dict[lbl] - f) + val2 = 0.0 + for lbl, f in vs_povm.items(): + val2 += _np.abs(t_dict[lbl] - f) + val += (wt_prep * val1 + wt_povm * val2) + ret += val elif spam_metric == "tracedist": # Leakage-aware metrics NOT available. diff --git a/pygsti/protocols/gst.py b/pygsti/protocols/gst.py index 3a3382283..66aef4162 100644 --- a/pygsti/protocols/gst.py +++ b/pygsti/protocols/gst.py @@ -2426,7 +2426,7 @@ def _memberdicts(_m) -> tuple[dict,dict,dict,dict]: argdicts = gaugeopt_suite.gaugeopt_argument_dicts.get(lbl, dict()) n_leak = 0 if isinstance(argdicts, list) and len(argdicts) > 0: - n_leak = argdicts[-1].get('n_leak', n_leak) + n_leak = argdicts[0].get('n_leak', n_leak) if n_leak > 0: dim = gaugeopt_model.basis.state_space.udim diff --git a/pygsti/tools/leakage.py b/pygsti/tools/leakage.py index 9bdf19b6c..bbed30d76 100644 --- a/pygsti/tools/leakage.py +++ b/pygsti/tools/leakage.py @@ -417,13 +417,11 @@ def lagoified_gopparams_dicts(gopparams_dicts: List[Dict]) -> List[Dict]: strip out gauge optimization over the TPSpam gauge group and apply a few other changes appropriate for leakage modeling. """ + from pygsti.models.gaugegroup import UnitaryGaugeGroup gopparams_dicts = [gp for gp in gopparams_dicts if 'TPSpam' not in str(type(gp['_gaugeGroupEl']))] gopparams_dicts = copy.deepcopy(gopparams_dicts) + tm = gopparams_dicts[0]['target_model'] for inner_dict in gopparams_dicts: - gg = inner_dict['target_model'].default_gauge_group - if gg is not None: - inner_dict['gauge_group'] = gg - inner_dict['_gaugeGroupEl'] = gg.compute_element(gg.initial_params) 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. @@ -432,18 +430,35 @@ 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['n_leak'] = 1 + inner_dict['item_weights'] = {'gates': 1.0, 'spam': 1.0} + gopparams_dicts.append(inner_dict) return gopparams_dicts From d6953d76a0e1705ac7c564c52319473ddf835139 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Wed, 8 Oct 2025 13:07:22 -0700 Subject: [PATCH 11/18] make note of idea without making half-baked changes to code --- pygsti/report/reportables.py | 22 ++++++++++++++-------- 1 file changed, 14 insertions(+), 8 deletions(-) diff --git a/pygsti/report/reportables.py b/pygsti/report/reportables.py index b145c5d53..9dd7a32c2 100644 --- a/pygsti/report/reportables.py +++ b/pygsti/report/reportables.py @@ -1600,17 +1600,23 @@ def eigenvalue_entanglement_infidelity(a, b, mx_basis): float """ d2 = a.shape[0] - evalsA, evecsA = _np.linalg.eig(a) - evalsB, evecsB = _np.linalg.eig(b) + evA = _np.linalg.eigvals(a) + evB = _np.linalg.eigvals(b) """ - IDEA: check if (a,b) are normal operators. If they are, then compute - proceed by computing a Schur decomposition to get orthonormal eigenbases. - From there, match them based on a dissimilarity kernel of abs(1 - ). - That could be used + IDEA: check if (a,b) are normal operators. If they are, then proceed by + computing a Schur decomposition to get orthonormal eigenbases. From there, + match them based on a dissimilarity kernel of abs(1 - ), rather than + the dissimilarity kernel abs(ev_a - ev_b). + + This would introduce some gauge dependence to eigenvalue entanglement + fidelity. However, there's already ambiguity in the definition of this + metric since it requires an ordering of the eigenvalues. It just so + happens that this function's current implementation chooses the ordering + in a gauge-invariant way. """ - _, pairs = _tools.minweight_match(evalsA, evalsB, lambda x, y: abs(x - y), + _, pairs = _tools.minweight_match(evA, evB, lambda x, y: abs(x - y), return_pairs=True) # just to get pairing - mlPl = abs(_np.sum([_np.conjugate(evalsB[j]) * evalsA[i] for i, j in pairs])) + mlPl = abs(_np.sum([_np.conjugate(evB[j]) * evA[i] for i, j in pairs])) return 1.0 - mlPl / float(d2) From 2de11625fd14b860574a5f5491ab0f2e5829cd86 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Wed, 8 Oct 2025 16:37:51 -0700 Subject: [PATCH 12/18] =?UTF-8?q?changes=20to=20leakage.py=20*=20Lots=20of?= =?UTF-8?q?=20docs.=20*=20New=20function=20called=20=E2=80=9Csubspace=5Fpr?= =?UTF-8?q?ojector=E2=80=9D=20*=20rename=20functions=20only=20introduced?= =?UTF-8?q?=20in=20this=20PR=20(leakage=5Fprofile,=20seepage=5Fprofile=20?= =?UTF-8?q?=E2=80=94>=20gate=5Fleakage=5Fprofile,=20gate=5Fseepage=5Fprofi?= =?UTF-8?q?le).=20*=20rename=20std=5Flago=5Fgopsuite=20to=20std=5Flago=5Fg?= =?UTF-8?q?augeopt=5Fparams.=20This=20makes=20sense=20because=20each=20dic?= =?UTF-8?q?t=20in=20a=20list-of-dicts=20representation=20of=20a=20gauge=20?= =?UTF-8?q?optimization=20suite=20is=20supposed=20to=20be=20passable=20to?= =?UTF-8?q?=20gaugeopt=5Fto=5Ftarget=20entirely=20as=20kwargs.=20One=20cou?= =?UTF-8?q?ld=20argue=20that=20this=20is=20a=20breaking=20change,=20but=20?= =?UTF-8?q?it's=20for=20a=20function=20that=20should=20have=20had=20low-ex?= =?UTF-8?q?posure=20in=200.9.14.?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit updates to gst.py and gaugeopt.py to incorporate the subspace_projector function. updates to reportables.py to reflect the renames of leakage_profile / seepage_profile. --- pygsti/algorithms/gaugeopt.py | 10 +- pygsti/protocols/gst.py | 9 +- pygsti/report/reportables.py | 12 +- pygsti/tools/leakage.py | 274 +++++++++++++++++++++------------- 4 files changed, 175 insertions(+), 130 deletions(-) diff --git a/pygsti/algorithms/gaugeopt.py b/pygsti/algorithms/gaugeopt.py index c68549e21..75077bcfa 100644 --- a/pygsti/algorithms/gaugeopt.py +++ b/pygsti/algorithms/gaugeopt.py @@ -609,15 +609,7 @@ def _mock_objective_fn(v): 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 + P = _tools.subspace_projector(dim - n_leak, dim, mxBasis) transform_mx_arg = (P, _tools.matrixtools.IdentityOperator()) # ^ The semantics of this tuple are defined by the frobeniusdist function # in the ExplicitOpModelCalc class. diff --git a/pygsti/protocols/gst.py b/pygsti/protocols/gst.py index 66aef4162..a8dfb4344 100644 --- a/pygsti/protocols/gst.py +++ b/pygsti/protocols/gst.py @@ -2432,14 +2432,7 @@ def _memberdicts(_m) -> tuple[dict,dict,dict,dict]: dim = gaugeopt_model.basis.state_space.udim assert n_leak == 1 assert dim == 3 - B = _tools.leading_dxd_submatrix_basis_vectors(2, 3, gaugeopt_model.basis) - P = B @ B.T.conj() - if _np.linalg.norm(P.imag) > 1e-12: - msg = f"Attempting to run leakage-aware gauge optimization with basis {gaugeopt_model.basis}\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) - P = P.real + P = _tools.subspace_projector(2, 3, gaugeopt_model.basis) else: P = _tools.matrixtools.IdentityOperator() diff --git a/pygsti/report/reportables.py b/pygsti/report/reportables.py index 9dd7a32c2..949379388 100644 --- a/pygsti/report/reportables.py +++ b/pygsti/report/reportables.py @@ -1070,8 +1070,8 @@ def subspace_diamonddist(op_a, op_b, basis): def pergate_leakrate_max(op, ignore, mx_basis): - from pygsti.tools.leakage import leakage_profile - evals, _ = leakage_profile(op, mx_basis) + from pygsti.tools.leakage import gate_leakage_profile + evals, _ = gate_leakage_profile(op, mx_basis) return evals[0] @@ -1079,8 +1079,8 @@ def pergate_leakrate_max(op, ignore, mx_basis): def pergate_leakrate_min(op, ignore, mx_basis): - from pygsti.tools.leakage import leakage_profile - evals, _ = leakage_profile(op, mx_basis) + from pygsti.tools.leakage import gate_leakage_profile + evals, _ = gate_leakage_profile(op, mx_basis) return evals[-1] @@ -1088,8 +1088,8 @@ def pergate_leakrate_min(op, ignore, mx_basis): def pergate_seeprate(op, ignore, mx_basis): - from pygsti.tools.leakage import seepage_profile - seeprate, _ = seepage_profile(op, mx_basis) + from pygsti.tools.leakage import gate_seepage_profile + seeprate, _ = gate_seepage_profile(op, mx_basis) return seeprate[0] diff --git a/pygsti/tools/leakage.py b/pygsti/tools/leakage.py index bbed30d76..14c1e44dc 100644 --- a/pygsti/tools/leakage.py +++ b/pygsti/tools/leakage.py @@ -27,8 +27,6 @@ from pygsti.processors import QubitProcessorSpec - -# Question: include the parenthetical in the heading below? NOTATION = \ """ Default notation (deferential to text above) @@ -145,71 +143,22 @@ def apply_tensorized_to_teststate(op_x: np.ndarray, op_y, op_basis: np.ndarray, @lru_cache -def leading_dxd_submatrix_basis_vectors(d: int, n: int, current_basis: Basis) -> np.ndarray: - """ - 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: +@set_docstring( +""" +Here, H has dimension n and C ⊂ H has dimension d ≤ n. - 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 +This function returns a column-unitary matrix B where P = B B^{\\dagger} is the +orthogonal projector from M[H] to M[C] with respect to current_basis. - 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. - """ +If you only care about P, then you can call subspace_projector instead. +""" + NOTATION) +def leading_dxd_submatrix_basis_vectors(d: int, n: int, current_basis: Basis) -> np.ndarray: assert d <= n + if d == n: + return np.eye(n**2) 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 @@ -223,6 +172,41 @@ def leading_dxd_submatrix_basis_vectors(d: int, n: int, current_basis: Basis) -> return submatrix_basis_vectors +@lru_cache +@set_docstring( +""" +This function returns the superoperator in S[H] that projects orthogonally from +M[H] to M[C], where H is n-dimensional and C ⊂ H is d-dimensional (d ≤ n). + +The action of this operator is easy to understand when M[H] and M[C] are viewed +as spaces of n-by-n matrices rather than spaces of length-n^2 vectors. + +For v in M[H] and u = P v, we have + + mat(v) = [x11, x12] and mat(u) = [x11, 0] + [x21, x22] [ 0, 0]. + +This characterization makes two facts about P apparent. First, P is positive +(i.e., it takes Hermitian psd operators to Hermitian psd operators). Second, +P is trace-non-increasing. +""" + NOTATION) +def subspace_projector(d: int, n: int, basis: Basis, force_real=True) -> np.ndarray: + assert d <= n + if d == n: + return np.eye(n**2) + U = leading_dxd_submatrix_basis_vectors(d, n, basis) # type: ignore + P = U @ U.T.conj() + if force_real: + if np.linalg.norm(P.imag) > 1e-12: + msg = "The orthogonal projector onto the computational subspace in the basis\n" + msg += f"{basis} is not real-valued. Since we were passed force_real=True we're\n" + msg += "raising a ValueError. Try again with a basis like 'l2p1' or 'gm', or\n" + msg += "use force_real=False." + raise ValueError(msg) + P = P.real + return P + + CHOI_INDUCED_METRIC = \ """ The pair (op_x, op_y) represent some superoperators (X, Y) in S[H], using op_basis. @@ -243,7 +227,7 @@ def subspace_entanglement_fidelity(op_x: np.ndarray, op_y: np.ndarray, op_basis, 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 + return ent_fid # type: ignore @set_docstring(CHOI_INDUCED_METRIC % 'jamiolkowski trace distance') @@ -252,45 +236,91 @@ def subspace_jtracedist(op_x: np.ndarray, op_y: np.ndarray, op_basis, n_leak=0) 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 + return j_dist # type: ignore -@set_docstring( + +PROJECTION_INDUCED_METRIC = \ """ 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. +We return the %s 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) +""" + NOTATION + + +@set_docstring(PROJECTION_INDUCED_METRIC % 'Frobenius distance') def subspace_superop_fro_dist(op_x: np.ndarray, op_y: np.ndarray, op_basis, n_leak=0) -> float: diff = op_x - op_y if n_leak == 0: - return np.linalg.norm(diff, 'fro') + return np.linalg.norm(diff, 'fro') # type: ignore + n = int(np.sqrt(op_x.shape[0])) + assert op_x.shape == op_y.shape == (n**2, n**2) + return np.linalg.norm(diff @ P) # type: ignore + + +@set_docstring(PROJECTION_INDUCED_METRIC % 'diamond distance') +def subspace_diamonddist(op_x: np.ndarray, op_y: np.ndarray, op_basis, n_leak=0) -> float: + """ + Here we give a brief motivating derivation for defining the subspace diamond norm in + the way that we have. This derivation won't convince a skeptic that our definition + is the best-possible. + + Suppose we canonically measure the distance between two superoperators (X, Y) by + + D(X, Y; H) = max || (X - Y) v || + v is in M[H], (Eq. 1) + tr(v) = 1, + v is positive + + for some norm || * ||. + + We arrive at a natural analog of this metric when (X, Y) are restricted to M[C] + simply by replacing "H" in (Eq. 1) with "C". - 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) + Using P to denote the orthogonal projector onto M[C], we claim that + D(X, Y; C) = D(X P, Y P; H). (Eq. 2) -def subspace_diamonddist(op_x: np.ndarray, op_y: np.ndarray, op_basis) -> float: + Here's a proof of that claim: + + | It's easy to show that P is a positive trace-non-increasing map. In particular, + | if u = P v, then the matrix representations of u and v are + | + | mat(v) = [v11, v12] and mat(u) = [v11, 0] + | [v21, v22] [ 0, 0], + | + | where v11 and v22 are psd if v is positive. From here the claim follows once + | you've convinced yourself that the pair of problems below have the same optimal + | objective value + | + | max || (X - Y) P v || and max || (X - Y) P v || + | mat(v) = [v11, v12] mat(v) = [v11, v12] + | [v21, v22] [v21, v22] + | mat(v) is PSD v11 is PSD + | tr(v11) + tr(v22) = 1 tr(v11) <= 1. + + This can be taken a little further. The proof's argument goes through unchanged if, + instead of starting with the objective || (X - Y) v ||, we started with f((X - Y) v), + where f satisfies the property that f(c v) >= f(v) whenever c is a scalar >= 1. + """ + from pygsti.tools.optools import diamonddist dim_mixed = op_x.shape[0] dim_pure = int(dim_mixed**0.5) - dim_pure_compsub = dim_pure - 1 - U = leading_dxd_submatrix_basis_vectors(dim_pure_compsub, dim_pure, op_basis) - P = U @ U.T.conj() - assert np.linalg.norm(P - P.real) < 1e-10 - P = P.real - from pygsti.tools.optools import diamonddist + assert n_leak <= 1 + dim_pure_compsub = dim_pure - n_leak + P = subspace_projector(dim_pure_compsub, dim_pure, op_basis) val : float = diamonddist(op_x @ P, op_y @ P, op_basis, return_x=False) / 2 # type: ignore return val -def leakage_profile(op, mx_basis) -> tuple[np.ndarray, list[np.ndarray]]: +def gate_leakage_profile(op, mx_basis) -> tuple[np.ndarray, list[np.ndarray]]: + """ + A + """ assert op.shape == (9, 9) lfb = BuiltinBasis('l2p1', 9) op_lfb = pgbt.change_basis(op, mx_basis, lfb) @@ -305,7 +335,7 @@ def leakage_profile(op, mx_basis) -> tuple[np.ndarray, list[np.ndarray]]: return rates, states -def seepage_profile(op, mx_basis) -> tuple[np.ndarray, list[np.ndarray]]: +def gate_seepage_profile(op, mx_basis) -> tuple[np.ndarray, list[np.ndarray]]: assert op.shape == (9, 9) lfb = BuiltinBasis('l2p1', 9) op_lfb = pgbt.change_basis(op, mx_basis, lfb) @@ -413,14 +443,48 @@ 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 `estimates.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) - tm = gopparams_dicts[0]['target_model'] for inner_dict in gopparams_dicts: inner_dict['n_leak'] = 1 # ^ This function could accept n_leak as an argument instead. However, @@ -462,15 +526,11 @@ def lagoified_gopparams_dicts(gopparams_dicts: List[Dict]) -> List[Dict]: return gopparams_dicts -def std_lago_gopsuite(model: ExplicitOpModel) -> dict[str, list[dict]]: +def std_lago_gaugeopt_params(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. + Return a dictionary of the form {'LAGO': v}, where v is a list-of-dicts + representation of a gauge optimization suite suitable for leakage modeling, + obtained by modiftying the 'stdgaugeopt' suite induced by `model`. """ from pygsti.protocols.gst import GSTGaugeOptSuite std_gop_suite = GSTGaugeOptSuite(gaugeopt_suite_names=('stdgaugeopt',)) @@ -484,23 +544,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] - 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) 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_gaugeopt_params(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 e4098f70d31cfc4d95440490b29eb54441da6ce4 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Thu, 9 Oct 2025 18:02:52 -0700 Subject: [PATCH 13/18] implement fully-general gate_leakage_profile and gate_seepage_profile functions. Rename subspace_projector to superop_subspace_projector. --- pygsti/algorithms/gaugeopt.py | 2 +- pygsti/baseobjs/basisconstructors.py | 4 +- pygsti/protocols/gst.py | 2 +- pygsti/tools/leakage.py | 182 +++++++++++++++++++++------ pygsti/tools/matrixtools.py | 21 +++- 5 files changed, 172 insertions(+), 39 deletions(-) diff --git a/pygsti/algorithms/gaugeopt.py b/pygsti/algorithms/gaugeopt.py index 75077bcfa..d0af57f3a 100644 --- a/pygsti/algorithms/gaugeopt.py +++ b/pygsti/algorithms/gaugeopt.py @@ -609,7 +609,7 @@ def _mock_objective_fn(v): dim = int(_np.sqrt(mxBasis.dim)) if n_leak > 0: - P = _tools.subspace_projector(dim - n_leak, dim, mxBasis) + P = _tools.superop_subspace_projector(dim - n_leak, dim, mxBasis) transform_mx_arg = (P, _tools.matrixtools.IdentityOperator()) # ^ The semantics of this tuple are defined by the frobeniusdist function # in the ExplicitOpModelCalc class. diff --git a/pygsti/baseobjs/basisconstructors.py b/pygsti/baseobjs/basisconstructors.py index 5ba4e40bd..7132c612b 100644 --- a/pygsti/baseobjs/basisconstructors.py +++ b/pygsti/baseobjs/basisconstructors.py @@ -10,6 +10,7 @@ # http://www.apache.org/licenses/LICENSE-2.0 or in the LICENSE file in the root pyGSTi directory. #*************************************************************************************************** import itertools as _itertools +import re as _re import numbers as _numbers import numpy as _np import scipy.sparse as _sps @@ -697,6 +698,7 @@ 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() @@ -718,6 +720,7 @@ def lf_labels(matrix_dim: int) -> tuple[str,...]: ) 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 @@ -743,7 +746,6 @@ def lf_matrices(matrix_dim: int) -> list[_np.ndarray]: return leakage_basis_mxs - def qsim_matrices(matrix_dim): """ Get the elements of the QuantumSim basis with matrix dimension `matrix_dim`. diff --git a/pygsti/protocols/gst.py b/pygsti/protocols/gst.py index a8dfb4344..ea169e1a7 100644 --- a/pygsti/protocols/gst.py +++ b/pygsti/protocols/gst.py @@ -2432,7 +2432,7 @@ def _memberdicts(_m) -> tuple[dict,dict,dict,dict]: dim = gaugeopt_model.basis.state_space.udim assert n_leak == 1 assert dim == 3 - P = _tools.subspace_projector(2, 3, gaugeopt_model.basis) + P = _tools.superop_subspace_projector(2, 3, gaugeopt_model.basis) else: P = _tools.matrixtools.IdentityOperator() diff --git a/pygsti/tools/leakage.py b/pygsti/tools/leakage.py index 14c1e44dc..f0822a4b6 100644 --- a/pygsti/tools/leakage.py +++ b/pygsti/tools/leakage.py @@ -14,10 +14,13 @@ 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 +from pygsti.baseobjs.basisconstructors import is_standard_leakage_basis_name import numpy as np +import scipy.linalg as la import warnings from typing import Union, Dict, Optional, List, Any, TYPE_CHECKING @@ -29,16 +32,19 @@ NOTATION = \ """ -Default notation (deferential to text above) --------------------------------------------- +Default notation (possibly superceded by 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. + standard basis vectors of H. The orthogonal complement of C in H is H \\ C. * 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. + We can interpret A ∈ M[C] as an element of M[H] by requiring that ker(A) and + ker(A^{†}) contain H \\ C. + * Given a space of linear operators, L, we write S[L] for the set of linear transformations ("superoperators") from L to L. @@ -150,7 +156,7 @@ def apply_tensorized_to_teststate(op_x: np.ndarray, op_y, op_basis: np.ndarray, This function returns a column-unitary matrix B where P = B B^{\\dagger} is the orthogonal projector from M[H] to M[C] with respect to current_basis. -If you only care about P, then you can call subspace_projector instead. +If you only care about P, then you can call superop_subspace_projector instead. """ + NOTATION) def leading_dxd_submatrix_basis_vectors(d: int, n: int, current_basis: Basis) -> np.ndarray: assert d <= n @@ -190,7 +196,7 @@ def leading_dxd_submatrix_basis_vectors(d: int, n: int, current_basis: Basis) -> (i.e., it takes Hermitian psd operators to Hermitian psd operators). Second, P is trace-non-increasing. """ + NOTATION) -def subspace_projector(d: int, n: int, basis: Basis, force_real=True) -> np.ndarray: +def superop_subspace_projector(d: int, n: int, basis: Basis, force_real=True) -> np.ndarray: assert d <= n if d == n: return np.eye(n**2) @@ -207,7 +213,7 @@ def subspace_projector(d: int, n: int, basis: Basis, force_real=True) -> np.ndar return P -CHOI_INDUCED_METRIC = \ +CHOI_INDUCED_METRIC_TEMPLATE = \ """ The pair (op_x, op_y) represent some superoperators (X, Y) in S[H], using op_basis. @@ -221,7 +227,7 @@ def subspace_projector(d: int, n: int, basis: Basis, force_real=True) -> np.ndar """ + NOTATION -@set_docstring(CHOI_INDUCED_METRIC % 'entanglement fidelity') +@set_docstring(CHOI_INDUCED_METRIC_TEMPLATE % 'entanglement fidelity') def subspace_entanglement_fidelity(op_x: np.ndarray, op_y: np.ndarray, op_basis, n_leak=0) -> float: 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) @@ -230,7 +236,7 @@ def subspace_entanglement_fidelity(op_x: np.ndarray, op_y: np.ndarray, op_basis, return ent_fid # type: ignore -@set_docstring(CHOI_INDUCED_METRIC % 'jamiolkowski trace distance') +@set_docstring(CHOI_INDUCED_METRIC_TEMPLATE % 'jamiolkowski trace distance') def subspace_jtracedist(op_x: np.ndarray, op_y: np.ndarray, op_basis, n_leak=0) -> float: 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) @@ -240,7 +246,7 @@ def subspace_jtracedist(op_x: np.ndarray, op_y: np.ndarray, op_basis, n_leak=0) -PROJECTION_INDUCED_METRIC = \ +PROJECTION_INDUCED_METRIC_TEMPLATE = \ """ The pair (op_x, op_y) represent some superoperators (X, Y) in S[H], using op_basis. @@ -252,17 +258,18 @@ def subspace_jtracedist(op_x: np.ndarray, op_y: np.ndarray, op_basis, n_leak=0) """ + NOTATION -@set_docstring(PROJECTION_INDUCED_METRIC % 'Frobenius distance') +@set_docstring(PROJECTION_INDUCED_METRIC_TEMPLATE % 'Frobenius distance') def subspace_superop_fro_dist(op_x: np.ndarray, op_y: np.ndarray, op_basis, n_leak=0) -> float: diff = op_x - op_y if n_leak == 0: return np.linalg.norm(diff, 'fro') # type: ignore n = int(np.sqrt(op_x.shape[0])) assert op_x.shape == op_y.shape == (n**2, n**2) + P = superop_subspace_projector(n - n_leak, n, op_basis, force_leak=False) return np.linalg.norm(diff @ P) # type: ignore -@set_docstring(PROJECTION_INDUCED_METRIC % 'diamond distance') +@set_docstring(PROJECTION_INDUCED_METRIC_TEMPLATE % 'diamond distance') def subspace_diamonddist(op_x: np.ndarray, op_y: np.ndarray, op_basis, n_leak=0) -> float: """ Here we give a brief motivating derivation for defining the subspace diamond norm in @@ -312,38 +319,143 @@ def subspace_diamonddist(op_x: np.ndarray, op_y: np.ndarray, op_basis, n_leak=0) dim_pure = int(dim_mixed**0.5) assert n_leak <= 1 dim_pure_compsub = dim_pure - n_leak - P = subspace_projector(dim_pure_compsub, dim_pure, op_basis) + P = superop_subspace_projector(dim_pure_compsub, dim_pure, op_basis) val : float = diamonddist(op_x @ P, op_y @ P, op_basis, return_x=False) / 2 # type: ignore return val -def gate_leakage_profile(op, mx_basis) -> tuple[np.ndarray, list[np.ndarray]]: - """ - A - """ - assert op.shape == (9, 9) - lfb = BuiltinBasis('l2p1', 9) - op_lfb = pgbt.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 = pgbt.vec_to_stdmx(leakage_effect_superket, 'pp') - leakage_rates, states = np.linalg.eigh(leakage_effect) - ind = np.argsort(leakage_rates)[::-1] - rates = leakage_rates[ind] - states = [np.concatenate([s,[0.0]]) for s in states.T[ind]] +TRANSPORT_PROFILES_TEMPLATE = \ +""" +This function returns a description of how a gate (in some process matrix representation) +transports population from SUB to the orthogonal complement of SUB in H. + +Underlying mathematics +---------------------- +Our subspace of interest is represented by the unique Hermitian operator E_SUB satisfying +⟨ E_SUB, ρ ⟩ = 1 and ρ = E_SUB ρ E_SUB for all densities ρ ∈ M[SUB]. %s + +The process matrix `op` is the mx_basis representation of a CPTP map G : M[H] ➞ M[H]. + +Consider an experimental protocol where a system is prepared in state ρ, evolved to G(ρ), +and then measured with the 2-element POVM {E_SUB, 1 - E_SUB}. Runs of this protocol that +result in the `1 - E_SUB` measurement outcome are called transport events. + +The population transport profile of (G,SUB) is a full specification of the probabilties +of transport events, considering all possible choices of ρ ∈ M[SUB]. This profile is +encoded in the Hermitian operator + + E_transport := (E_SUB) G^{†}(1 - E_SUB) (E_SUB) ∈ M[SUB] ⊂ M[H]. + +This is a valid representation, since + + Pr{ G transports ρ | ρ ∈ M[SUB] } + = ⟨ 1 - E_SUB , G(ρ) ⟩ // definition of a transport event + = ⟨ G^{†}(1 - E_SUB) , ρ ⟩ // definition of the adjoint + = ⟨ E_transport , ρ ⟩. // using ρ = E_SUB ρ E_SUB + + +Function output and interpretation +---------------------------------- +This function returns a tuple with eigenvalues and eigenvectors of E_transport. + +Say this tuple is `(rates, states)`. If G is CPTP, then + + rates[-1] ≤ Pr{ G transports ρ | ρ ∈ M[SUB] } ≤ rates[0]. + +The upper bound is G's maximum transport of population (Max TOP) out of SUB, and +is acheived by applying G to the pure state `states[0]` ∈ SUB. + +The lower bound is G's minimum transport of population (Min TOP) out of SUB (assuming +population starts in SUB) and is acheived by applying G to `states[-1]` ∈ SUB. + +The length of `rates` and `states` is equal to the rank of E_SUB. + +""" + + +@set_docstring( +TRANSPORT_PROFILES_TEMPLATE.replace('SUB', 'sub') % '' + """ +Failure modes +------------- +We raise a ValueError if pygsti.tools.is_projector(E_sub, E_sub_tol) returns False. +""" + NOTATION ) +def pop_transport_profile(E_sub: np.ndarray, op: np.ndarray, mx_basis: Basis, E_sub_tol=1e-14) -> tuple[np.ndarray,list[np.ndarray]]: + n = int(np.sqrt(E_sub.size)) + if not pgmt.is_projector(E_sub, E_sub_tol): + msg = \ + f""" + The argument E_sub must be an orthogonal projector. The provided value was + + E_sub = {E_sub}; + + this failed the tests in pygsti.tools.is_projector at tolerance={E_sub_tol}. + """ + raise ValueError(msg) + + E_sub_perp_mat = np.eye(n) - E_sub + E_sub_perp_vec = pgbt.stdmx_to_vec(E_sub_perp_mat, mx_basis) + transport_E_vec = op.T @ E_sub_perp_vec + transport_E_mat = pgbt.vec_to_stdmx(transport_E_vec, mx_basis, keep_complex=True) + transport_E_mat = E_sub @ transport_E_mat @ E_sub + + rates, states = np.linalg.eigh(transport_E_mat) + dim_proj = int(np.trace(E_sub).real) + ind = np.argsort(np.abs(rates))[::-1][:dim_proj] + rates = rates[ind] + states = [s for s in states.T[ind]] + + return rates, states + + +@set_docstring( +TRANSPORT_PROFILES_TEMPLATE.replace('SUB', 'C') % """This operator is +proportional to the element of mx_basis labeled 'I', since that basis element is assumed +be proportional to the orthogonal projector onto C. +""" + NOTATION ) +def gate_leakage_profile(op: np.ndarray, mx_basis: Basis) -> tuple[np.ndarray, list[np.ndarray]]: + mx_basis = Basis.cast(mx_basis, dim=int(np.sqrt(op.size))) + E_comp_mat = mx_basis.ellookup['I'] # type: ignore + n = int(np.sqrt(E_comp_mat.size)) + assert E_comp_mat.shape == (n, n) + dim_comp = pgmt.matrix_rank(E_comp_mat, assume_hermitian=True) + E_comp_mat = E_comp_mat * (dim_comp / np.trace(E_comp_mat)) + if dim_comp == n: + msg = \ + """ + The provided basis suggests that the computational subspace is equal to + the entire system Hilbert space. Returning with an empty leakage profile! + """ + warnings.warn(msg) + return np.empty((0,)), [] + + rates, states = pop_transport_profile(E_comp_mat, op, mx_basis) return rates, states +@set_docstring( +TRANSPORT_PROFILES_TEMPLATE.replace('SUB', '[H \\ C]') % """This operator is +inferred from the element of mx_basis labeled 'I', since that basis element is assumed +be proportional to the orthogonal projector onto C. +""" + NOTATION ) def gate_seepage_profile(op, mx_basis) -> tuple[np.ndarray, list[np.ndarray]]: - assert op.shape == (9, 9) - lfb = BuiltinBasis('l2p1', 9) - op_lfb = pgbt.change_basis(op, mx_basis, lfb) - elinds = lfb.elindlookup - seeprate = np.atleast_1d(op_lfb[elinds['I'], elinds['L']]) - state = [np.array([0.0, 0.0, 1.0], dtype=np.complex128)] - return seeprate, state + mx_basis = Basis.cast(mx_basis, dim=int(np.sqrt(op.size))) + E_comp_mat = mx_basis.ellookup['I'] # type: ignore + n = int(np.sqrt(E_comp_mat.size)) + assert E_comp_mat.shape == (n, n) + dim_comp = pgmt.matrix_rank(E_comp_mat, assume_hermitian=True) + E_comp_mat = E_comp_mat * (dim_comp / np.trace(E_comp_mat)) + E_leak_mat = np.eye(n) - E_comp_mat + if dim_comp == n: + msg = \ + """ + The provided basis suggests that the computational subspace is equal to + the entire system Hilbert space. Returning with an empty seepage profile! + """ + warnings.warn(msg) + return np.empty((0,)), [] + rates, states = pop_transport_profile(E_leak_mat, op, mx_basis) + return rates, states # MARK: model construction @@ -478,7 +590,7 @@ def lagoified_gopparams_dicts(gopparams_dicts: List[Dict]) -> List[Dict]: # add_lago_models(results, 'CPTPLND', gos) - After those lines execute, the `estimates.models` dict will have two new + 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 diff --git a/pygsti/tools/matrixtools.py b/pygsti/tools/matrixtools.py index 210ef23cd..26fd62bf0 100644 --- a/pygsti/tools/matrixtools.py +++ b/pygsti/tools/matrixtools.py @@ -139,6 +139,18 @@ def is_valid_density_mx(mx, tol=1e-9): return abs(_np.trace(mx) - 1.0) < tol and is_pos_def(mx, tol) +def is_projector(mx, tol=1e-14): + # We use a stringent default tolerance since projectors tend to be + # computed to high accuracy. + if tol == _np.inf: + return True + if not is_hermitian(mx, tol): + return False + mx2 = mx @ mx + atol = _np.sqrt(mx.size) * tol + return _np.allclose(mx, mx2, atol=atol) + + def nullspace(m, tol=1e-7): """ Compute the nullspace of a matrix. @@ -473,6 +485,14 @@ def pinv_of_matrix_with_orthogonal_columns(m): return m_with_scaled_cols.T +def matrix_rank(m : _np.ndarray, tol=1e-14, assume_hermitian=False): + if assume_hermitian: + s = _np.abs(_spl.eigvalsh(m)) + else: + s = _spl.svdvals(m) + return _np.count_nonzero(s > tol*_np.max(s)) + + def matrix_sign(m): """ Compute the matrix s = sign(m). The eigenvectors of s are the same as those of m. @@ -533,7 +553,6 @@ def matrix_sign(m): return sign - def print_mx(mx, width=9, prec=4, withbrackets=False): """ Print matrix in pretty format. From 9b90bae3728e7da71137ae360179c8181474bbb6 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Fri, 10 Oct 2025 13:31:06 -0700 Subject: [PATCH 14/18] test for matrix_tools.py:is_projector. Use np.linalg.matrix_rank in leakage.py, rather than home-brewing an equivalent function. --- pygsti/baseobjs/basisconstructors.py | 1 - pygsti/tools/leakage.py | 7 ++----- pygsti/tools/matrixtools.py | 12 ++---------- test/unit/tools/test_matrixtools.py | 20 ++++++++++++++++++++ 4 files changed, 24 insertions(+), 16 deletions(-) diff --git a/pygsti/baseobjs/basisconstructors.py b/pygsti/baseobjs/basisconstructors.py index 7132c612b..f37473770 100644 --- a/pygsti/baseobjs/basisconstructors.py +++ b/pygsti/baseobjs/basisconstructors.py @@ -10,7 +10,6 @@ # http://www.apache.org/licenses/LICENSE-2.0 or in the LICENSE file in the root pyGSTi directory. #*************************************************************************************************** import itertools as _itertools -import re as _re import numbers as _numbers import numpy as _np import scipy.sparse as _sps diff --git a/pygsti/tools/leakage.py b/pygsti/tools/leakage.py index f0822a4b6..ead7f0dce 100644 --- a/pygsti/tools/leakage.py +++ b/pygsti/tools/leakage.py @@ -18,7 +18,6 @@ from pygsti.tools.basistools import stdmx_to_vec from pygsti.baseobjs import Label from pygsti.baseobjs.basis import TensorProdBasis, Basis, BuiltinBasis -from pygsti.baseobjs.basisconstructors import is_standard_leakage_basis_name import numpy as np import scipy.linalg as la import warnings @@ -417,7 +416,7 @@ def gate_leakage_profile(op: np.ndarray, mx_basis: Basis) -> tuple[np.ndarray, l E_comp_mat = mx_basis.ellookup['I'] # type: ignore n = int(np.sqrt(E_comp_mat.size)) assert E_comp_mat.shape == (n, n) - dim_comp = pgmt.matrix_rank(E_comp_mat, assume_hermitian=True) + dim_comp = np.linalg.matrix_rank(E_comp_mat, hermitian=True) E_comp_mat = E_comp_mat * (dim_comp / np.trace(E_comp_mat)) if dim_comp == n: msg = \ @@ -442,7 +441,7 @@ def gate_seepage_profile(op, mx_basis) -> tuple[np.ndarray, list[np.ndarray]]: E_comp_mat = mx_basis.ellookup['I'] # type: ignore n = int(np.sqrt(E_comp_mat.size)) assert E_comp_mat.shape == (n, n) - dim_comp = pgmt.matrix_rank(E_comp_mat, assume_hermitian=True) + dim_comp = np.linalg.matrix_rank(E_comp_mat, hermitian=True) E_comp_mat = E_comp_mat * (dim_comp / np.trace(E_comp_mat)) E_leak_mat = np.eye(n) - E_comp_mat if dim_comp == n: @@ -726,5 +725,3 @@ def construct_leakage_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 26fd62bf0..8be082733 100644 --- a/pygsti/tools/matrixtools.py +++ b/pygsti/tools/matrixtools.py @@ -147,8 +147,8 @@ def is_projector(mx, tol=1e-14): if not is_hermitian(mx, tol): return False mx2 = mx @ mx - atol = _np.sqrt(mx.size) * tol - return _np.allclose(mx, mx2, atol=atol) + tol = _np.sqrt(mx.size) * tol + return _np.allclose(mx, mx2, atol=tol, rtol=tol) def nullspace(m, tol=1e-7): @@ -485,14 +485,6 @@ def pinv_of_matrix_with_orthogonal_columns(m): return m_with_scaled_cols.T -def matrix_rank(m : _np.ndarray, tol=1e-14, assume_hermitian=False): - if assume_hermitian: - s = _np.abs(_spl.eigvalsh(m)) - else: - s = _spl.svdvals(m) - return _np.count_nonzero(s > tol*_np.max(s)) - - def matrix_sign(m): """ Compute the matrix s = sign(m). The eigenvectors of s are the same as those of m. diff --git a/test/unit/tools/test_matrixtools.py b/test/unit/tools/test_matrixtools.py index 3b2b20a75..0eb31c209 100644 --- a/test/unit/tools/test_matrixtools.py +++ b/test/unit/tools/test_matrixtools.py @@ -24,6 +24,25 @@ def test_is_pos_def(self): self.assertTrue(mt.is_pos_def(pos_mx)) self.assertFalse(mt.is_pos_def(non_pos_mx)) + def test_is_projector(self): + self.assertTrue(mt.is_projector(np.zeros((3, 3)))) + I2 = np.eye(2) + self.assertTrue(mt.is_projector(I2)) + I2[1,1] = 0.9999998 + self.assertFalse(mt.is_projector(I2)) + self.assertTrue( mt.is_projector(I2, tol=np.inf)) + A = np.array([ + [-0.653829-0.313923j, -0.129614+1.458021j, 0.783975+1.960258j], + [ 1.493431+1.801635j, -1.259066+1.315104j, 1.513924+0.35738j ], + [ 1.345875-1.208319j, 0.781311-0.004454j, 0.264456+0.656475j] + ]) + Q = np.linalg.qr(A)[0] + for i in range(1, 4): + U = Q[:,:i].reshape((3, i)) + P = U @ U.T.conj() + self.assertTrue(mt.is_projector(P)) + return + def test_mx_to_string(self): mx = np.array([[ 1, 1+2j], [1-2j, 3]], 'complex') @@ -192,3 +211,4 @@ def test_prime_factors(self): self.assertEqual(mt.prime_factors(7), [7]) self.assertEqual(mt.prime_factors(10), [2, 5]) self.assertEqual(mt.prime_factors(12), [2, 2, 3]) + From 06c3c36982e11fb70dd794ef359027be16d71cf4 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Fri, 10 Oct 2025 14:51:37 -0700 Subject: [PATCH 15/18] correct errors in recent changes to fidelity-based gauge optimization --- pygsti/algorithms/gaugeopt.py | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/pygsti/algorithms/gaugeopt.py b/pygsti/algorithms/gaugeopt.py index d0af57f3a..fab15cce6 100644 --- a/pygsti/algorithms/gaugeopt.py +++ b/pygsti/algorithms/gaugeopt.py @@ -625,11 +625,12 @@ def _mock_objective_fn(v): (1) its fidelities with the target gates match the eigenvalue fidelities with the target gates; - (2) the fidelities between rho_e and {E_t : E_t in M_target} will match those between - rho_e and {E_e : E_e in M_estimate}; and + (2) the (rho_e, M_t) SPAM probs will match the (rho_e, M_e) SPAM probs + + (3) the (rho_t, M_e) SPAM probs will match the (rho_e, M_e) SPAM probs + + ^ (3) originally said "fidelities ... match those between rho_e and {E_t : E_t in M_target}", which is WRONG. - (3) the fidelities between rho_t and {E_e : E_e in M_estimate} will match those between - rho_t and {E_t : E_t in M_target}. In view of this, taking fidelity as the gauge-optimization objective tries to minimize an aggregration of any mismatch between the various fidelities above. @@ -652,7 +653,7 @@ def _mock_objective_fn(v): for povmlbl in target_model.povms: rho_curest = model.preps[preplbl].to_dense() M_curest = model.povms[povmlbl] - t = {elbl: vec_fidelity(rho_curest, e.to_dense(), mxBasis).item() for (elbl, e) in M_curest.items() } + t = {elbl: _np.vdot(rho_curest, e.to_dense()) for (elbl, e) in M_curest.items() } spam_fidelity_targets[(preplbl, povmlbl)] = t def _objective_fn(gauge_group_el, oob_check): @@ -734,8 +735,8 @@ def _objective_fn(gauge_group_el, oob_check): rho_target = target_model.preps[preplbl].to_dense() M_curest = model.povms[povmlbl] M_target = target_model.povms[povmlbl] - vs_prep = {elbl: vec_fidelity(rho_curest, e.to_dense(), mxBasis) for (elbl, e) in M_target.items() } - vs_povm = {elbl: vec_fidelity(rho_target, e.to_dense(), mxBasis) for (elbl, e) in M_curest.items() } + vs_prep = {elbl: _np.vdot(rho_curest, e_target.to_dense()) for (elbl, e_target) in M_target.items() } + vs_povm = {elbl: _np.vdot(rho_target, e_curest.to_dense()) for (elbl, e_curest) in M_curest.items() } t_dict = spam_fidelity_targets[(preplbl, povmlbl)] val1 = 0.0 for lbl, f in vs_prep.items(): From a7683287356752df987754a26beae4abb4c27877 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Wed, 15 Oct 2025 08:24:19 -0700 Subject: [PATCH 16/18] heavy refactor of gaugeopt_to_target as entry point for gauge optimization --- pygsti/algorithms/gaugeopt.py | 213 ++++++++++++++++++++++++++-------- 1 file changed, 163 insertions(+), 50 deletions(-) diff --git a/pygsti/algorithms/gaugeopt.py b/pygsti/algorithms/gaugeopt.py index fab15cce6..851aed71c 100644 --- a/pygsti/algorithms/gaugeopt.py +++ b/pygsti/algorithms/gaugeopt.py @@ -20,21 +20,108 @@ from pygsti import tools as _tools from pygsti.tools import mpitools as _mpit from pygsti.tools import slicetools as _slct +from pygsti.models import ( + ExplicitOpModel as _ExplicitOpModel +) +from pygsti.models.model import OpModel as _OpModel from pygsti.models.gaugegroup import ( TrivialGaugeGroupElement as _TrivialGaugeGroupElement, GaugeGroupElement as _GaugeGroupElement ) -from typing import Callable, Union, Optional +from typing import Callable, Union, Optional, Any + + +class GaugeoptToTargetArgs: + # This class is a glorified namespace. It was introduced to make gaugeopt_to_target(...) easier to understand. + + @staticmethod + def parsed_model(model: Union[_OpModel, _ExplicitOpModel], convert_model_to: Optional[Any]) -> _OpModel: + if convert_model_to is None: + return model + + if not isinstance(model, _ExplicitOpModel): + raise ValueError('Gauge optimization only supports model conversion for ExplicitOpModels.') + conversion_args = convert_model_to if isinstance(convert_model_to, (list, tuple)) else (convert_model_to,) + model_out = model.copy() # don't alter the original model's parameterization (this would be unexpected) + for args in conversion_args: + if isinstance(args, str): + model_out.convert_members_inplace(args, set_default_gauge_group=True) + elif isinstance(args, dict): + model_out.convert_members_inplace(**args) + else: + raise ValueError("Invalid `convert_model_to` arguments: %s" % str(args)) + return model_out + + @staticmethod + def parsed_method(target_model: Optional[_OpModel], method: str, gates_metric: str, spam_metric: str, n_leak: int) -> str: + ls_mode_allowed = bool(target_model is not None + and gates_metric.startswith("frobenius") + and spam_metric.startswith("frobenius") + and n_leak == 0) + # and model.dim < 64: # least squares optimization seems uneffective if more than 3 qubits + # -- observed by Lucas - should try to debug why 3 qubits seemed to cause trouble... + if method == "ls" and not ls_mode_allowed: + raise ValueError("Least-squares method is not allowed! Target" + " model must be non-None and frobenius metrics" + " must be used.") + elif method == "auto": + return 'ls' if ls_mode_allowed else 'L-BFGS-B' + else: + return method + + # The static dicts of default values are substituted in gaugeopt_to_target's kwargs. + # This is a safe thing to do because no invocation of "gaugeopt_to_target" within pyGSTi + # used positional arguments past target_model. + + create_objective_passthrough_kwargs : dict[str,Any] = dict( + item_weights=None, + cptp_penalty_factor=0, + spam_penalty_factor=0, + check_jac=False, + n_leak=0, + ) + + gaugeopt_custom_passthrough_kwargs : dict[str, Any] = dict( + maxiter=100000, + maxfev=None, + tol=1e-8, + oob_check_interval=0, + verbosity=0, + ) + + other_kwargs : dict[str, Any] = dict( + gates_metric="frobenius", # defines ls_mode_allowed, then passed to _create_objective_fn, + spam_metric="frobenius", # defines ls_mode_allowed, then passed to _create_objective_fn, + gauge_group=None, # used iff convert_model_to is not None, then passed to _create_objective_fn and gaugeopt_custom + method='auto', # validated, then passed to _create_objective_fn and gaugeopt_custom + return_all=False, # used in the function body (only for branching after passing to gaugeopt_custom) + convert_model_to=None, # passthrough to converted_model. + comm=None, # passthrough to _create_objective_fn and gaugeopt_custom + ) + + default_kwargs : dict[str, Any] = dict( + **create_objective_passthrough_kwargs, + **gaugeopt_custom_passthrough_kwargs, + **other_kwargs + ) + + @staticmethod + def create_full_kwargs(kwargs: dict[str, Any]): + full_kwargs = GaugeoptToTargetArgs.default_kwargs.copy() + full_kwargs.update(kwargs) + return full_kwargs + + old_positional_args = ( + 'item_weights', 'cptp_penalty_factor', 'spam_penalty_factor', + 'gates_metric', 'spam_metric', 'gauge_group', 'method', + 'maxiter', 'maxfev', 'tol', 'oob_check_interval', + 'convert_model_to', 'return_all', 'comm', 'verbosity', + 'check_jac' + ) -def gaugeopt_to_target(model, target_model, item_weights=None, - cptp_penalty_factor=0, spam_penalty_factor=0, - gates_metric="frobenius", spam_metric="frobenius", - 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, n_leak=0): +def gaugeopt_to_target(model, target_model, *args, **kwargs): """ Optimize the gauge degrees of freedom of a model to that of a target. @@ -145,43 +232,68 @@ def gaugeopt_to_target(model, target_model, item_weights=None, found, gaugeMx is the gauge matrix used to transform the model, and model is the final gauge-transformed model. """ - - if item_weights is None: item_weights = {} - - ls_mode_allowed = bool(target_model is not None - and gates_metric.startswith("frobenius") - and spam_metric.startswith("frobenius")) - #and model.dim < 64: # least squares optimization seems uneffective if more than 3 qubits - # -- observed by Lucas - should try to debug why 3 qubits seemed to cause trouble... + if extra_posargs := len(args) > 0: + msg = \ + f""" + Recieved {extra_posargs} positional arguments past `model` and `target_model`. + These should be passed as appropriate keyword arguments instead. This version + of pyGSTi will infer intended keyword arguments based on the legacy argument + positions. Future versions of pyGSTi will raise an error. + """ + _warnings.warn(msg) + kwargs.update({ k: v for (k, v) in zip(GaugeoptToTargetArgs.old_positional_args, args) }) + out = gaugeopt_to_target(model, target_model, **kwargs) + return out - if method == "ls" and not ls_mode_allowed: - raise ValueError("Least-squares method is not allowed! Target" - " model must be non-None and frobenius metrics" - " must be used.") - if method == "auto": - method = 'ls' if ls_mode_allowed else 'L-BFGS-B' + full_kwargs = GaugeoptToTargetArgs.create_full_kwargs(kwargs) + """ + This function handles a strange situation where `target_model` can be None. - if convert_model_to is not None: - conversion_args = convert_model_to if isinstance(convert_model_to, (list, tuple)) else (convert_model_to,) - model = model.copy() # don't alter the original model's parameterization (this would be unexpected) - for args in conversion_args: - if isinstance(args, str): - model.convert_members_inplace(args, set_default_gauge_group=True) - elif isinstance(args, dict): - model.convert_members_inplace(**args) - else: - raise ValueError("Invalid `convert_model_to` arguments: %s" % str(args)) + In this case, the objective function will only depend on `cptp_penality_factor` + and `spam_penalty_factor`; it's forbidden to use method == 'ls'; and the + returned OpModel might not have its basis set. + """ + # arg parsing: validating `method` + gates_metric = full_kwargs['gates_metric'] + spam_metric = full_kwargs['spam_metric'] + n_leak = full_kwargs['n_leak'] + method = full_kwargs['method'] + method = GaugeoptToTargetArgs.parsed_method( + target_model, method, gates_metric, spam_metric, n_leak + ) + + # arg parsing: (possibly) converting `model` + convert_model_to = full_kwargs['convert_model_to'] + model = GaugeoptToTargetArgs.parsed_model(model, convert_model_to) + + # actual work: constructing objective_fn and jacobian_fn + item_weights = full_kwargs['item_weights'] + if item_weights is None: + item_weights = {} + cptp_penalty = full_kwargs['cptp_penalty_factor'] + spam_penalty = full_kwargs['spam_penalty_factor'] + comm = full_kwargs['comm'] + check_jac = full_kwargs['check_jac'] objective_fn, jacobian_fn = _create_objective_fn( - model, target_model, item_weights, - cptp_penalty_factor, spam_penalty_factor, - gates_metric, spam_metric, method, comm, check_jac, n_leak) - - result = gaugeopt_custom(model, objective_fn, gauge_group, method, - maxiter, maxfev, tol, oob_check_interval, - return_all, jacobian_fn, comm, verbosity) - - #If we've gauge optimized to a target model, declare that the + model, target_model, item_weights, cptp_penalty, spam_penalty, + gates_metric, spam_metric, method, comm, check_jac, n_leak + ) + + # actual work: calling the (wrapper of the wrapper of the ...) optimizer + gauge_group = full_kwargs['gauge_group'] + maxiter = full_kwargs['maxiter'] + maxfev = full_kwargs['maxfev'] + tol = full_kwargs['tol'] + oob_check = full_kwargs['oob_check_interval'] + return_all = full_kwargs['return_all'] + verbosity = full_kwargs['verbosity'] + result = gaugeopt_custom( + model, objective_fn, gauge_group, method, maxiter, maxfev, + tol, oob_check, return_all, jacobian_fn, comm, verbosity + ) + + # If we've gauge optimized to a target model, declare that the # resulting model is now in the same basis as the target. if target_model is not None: newModel = result[-1] if return_all else result @@ -190,9 +302,14 @@ def gaugeopt_to_target(model, target_model, item_weights=None, return result -def gaugeopt_custom(model, objective_fn, gauge_group=None, +GGElObjective = Callable[[_GaugeGroupElement, bool], Union[float, _np.ndarray]] + +GGElJacobian = Union[None, Callable[[_GaugeGroupElement], _np.ndarray]] + + +def gaugeopt_custom(model, objective_fn: GGElObjective, gauge_group=None, method='L-BFGS-B', maxiter=100000, maxfev=None, tol=1e-8, - oob_check_interval=0, return_all=False, jacobian_fn=None, + oob_check_interval=0, return_all=False, jacobian_fn: Optional[GGElJacobian]=None, comm=None, verbosity=0): """ Optimize the gauge of a model using a custom objective function. @@ -203,8 +320,9 @@ def gaugeopt_custom(model, objective_fn, gauge_group=None, The model to gauge-optimize objective_fn : function - The function to be minimized. The function must take a single `Model` - argument and return a float. + The function to be minimized. The function must take a GaugeGroupElement + and a bool. If method == 'ls' then objective_fn must return an ndarray; if + method != 'ls' then objective_fn must return a float. gauge_group : GaugeGroup, optional The gauge group which defines which gauge trasformations are optimized @@ -356,11 +474,6 @@ def _call_jacobian_fn(gauge_group_el_vec): return newModel -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", From 276c2541c231722da6c3dba51e0a37803d2650a7 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Wed, 15 Oct 2025 08:51:56 -0700 Subject: [PATCH 17/18] tweaks in last commit`s refactor of gaugeopt_to_target --- pygsti/algorithms/gaugeopt.py | 183 +++++++++++++++++----------------- 1 file changed, 92 insertions(+), 91 deletions(-) diff --git a/pygsti/algorithms/gaugeopt.py b/pygsti/algorithms/gaugeopt.py index 851aed71c..80dd62ccd 100644 --- a/pygsti/algorithms/gaugeopt.py +++ b/pygsti/algorithms/gaugeopt.py @@ -33,7 +33,18 @@ class GaugeoptToTargetArgs: - # This class is a glorified namespace. It was introduced to make gaugeopt_to_target(...) easier to understand. + """ + This class is basically a namespace. It was introduced to strip out tons of complexity + in gaugeopt_to_target(...) without breaking old code that might call gaugeopt_to_target. + """ + + old_trailing_positional_args = ( + 'item_weights', 'cptp_penalty_factor', 'spam_penalty_factor', + 'gates_metric', 'spam_metric', 'gauge_group', 'method', + 'maxiter', 'maxfev', 'tol', 'oob_check_interval', + 'convert_model_to', 'return_all', 'comm', 'verbosity', + 'check_jac' + ) @staticmethod def parsed_model(model: Union[_OpModel, _ExplicitOpModel], convert_model_to: Optional[Any]) -> _OpModel: @@ -79,8 +90,30 @@ def parsed_method(target_model: Optional[_OpModel], method: str, gates_metric: s cptp_penalty_factor=0, spam_penalty_factor=0, check_jac=False, - n_leak=0, ) + """ + item_weights : dict + Dictionary of weighting factors for gates and spam operators. Keys can + be gate, state preparation, or POVM effect, as well as the special values + "spam" or "gates" which apply the given weighting to *all* spam operators + or gates respectively. Values are floating point numbers. Values given + for specific gates or spam operators take precedence over "gates" and + "spam" values. The precise use of these weights depends on the model + metric(s) being used. + + cptp_penalty_factor : float + If greater than zero, the objective function also contains CPTP penalty + terms which penalize non-CPTP-ness of the gates being optimized. This factor + multiplies these CPTP penalty terms. + + spam_penalty_factor : float + If greater than zero, the objective function also contains SPAM penalty + terms which penalize non-positive-ness of the state preps being optimized. This + factor multiplies these SPAM penalty terms. + + check_jac : bool + When True, check least squares analytic jacobian against finite differences. + """ gaugeopt_custom_passthrough_kwargs : dict[str, Any] = dict( maxiter=100000, @@ -89,6 +122,23 @@ def parsed_method(target_model: Optional[_OpModel], method: str, gates_metric: s oob_check_interval=0, verbosity=0, ) + """ + maxiter : int + Maximum number of iterations for the gauge optimization. + + maxfev : int + Maximum number of function evaluations for the gauge optimization. + Defaults to maxiter. + + tol : float + The tolerance for the gauge optimization. + + oob_check_interval : int + If greater than zero, gauge transformations are allowed to fail (by raising + any exception) to indicate an out-of-bounds condition that the gauge optimizer + will avoid. If zero, then any gauge-transform failures just terminate the + optimization. + """ other_kwargs : dict[str, Any] = dict( gates_metric="frobenius", # defines ls_mode_allowed, then passed to _create_objective_fn, @@ -96,63 +146,11 @@ def parsed_method(target_model: Optional[_OpModel], method: str, gates_metric: s gauge_group=None, # used iff convert_model_to is not None, then passed to _create_objective_fn and gaugeopt_custom method='auto', # validated, then passed to _create_objective_fn and gaugeopt_custom return_all=False, # used in the function body (only for branching after passing to gaugeopt_custom) - convert_model_to=None, # passthrough to converted_model. + convert_model_to=None, # passthrough to parsed_model. comm=None, # passthrough to _create_objective_fn and gaugeopt_custom + n_leak=0 # passthrough to parsed_objective and _create_objective_fn ) - - default_kwargs : dict[str, Any] = dict( - **create_objective_passthrough_kwargs, - **gaugeopt_custom_passthrough_kwargs, - **other_kwargs - ) - - @staticmethod - def create_full_kwargs(kwargs: dict[str, Any]): - full_kwargs = GaugeoptToTargetArgs.default_kwargs.copy() - full_kwargs.update(kwargs) - return full_kwargs - - old_positional_args = ( - 'item_weights', 'cptp_penalty_factor', 'spam_penalty_factor', - 'gates_metric', 'spam_metric', 'gauge_group', 'method', - 'maxiter', 'maxfev', 'tol', 'oob_check_interval', - 'convert_model_to', 'return_all', 'comm', 'verbosity', - 'check_jac' - ) - - -def gaugeopt_to_target(model, target_model, *args, **kwargs): """ - Optimize the gauge degrees of freedom of a model to that of a target. - - Parameters - ---------- - model : Model - The model to gauge-optimize - - target_model : Model - The model to optimize to. The metric used for comparing models - is given by `gates_metric` and `spam_metric`. - - item_weights : dict, optional - Dictionary of weighting factors for gates and spam operators. Keys can - be gate, state preparation, or POVM effect, as well as the special values - "spam" or "gates" which apply the given weighting to *all* spam operators - or gates respectively. Values are floating point numbers. Values given - for specific gates or spam operators take precedence over "gates" and - "spam" values. The precise use of these weights depends on the model - metric(s) being used. - - cptp_penalty_factor : float, optional - If greater than zero, the objective function also contains CPTP penalty - terms which penalize non-CPTP-ness of the gates being optimized. This factor - multiplies these CPTP penalty terms. - - spam_penalty_factor : float, optional - If greater than zero, the objective function also contains SPAM penalty - terms which penalize non-positive-ness of the state preps being optimized. This - factor multiplies these SPAM penalty terms. - gates_metric : {"frobenius", "fidelity", "tracedist"}, optional The metric used to compare gates within models. "frobenius" computes the normalized sqrt(sum-of-squared-differences), with weights @@ -185,22 +183,6 @@ def gaugeopt_to_target(model, target_model, *args, **kwargs): - 'evolve' -- evolutionary global optimization algorithm using DEAP - 'brute' -- Experimental: scipy.optimize.brute using 4 points along each dimensions - maxiter : int, optional - Maximum number of iterations for the gauge optimization. - - maxfev : int, optional - Maximum number of function evaluations for the gauge optimization. - Defaults to maxiter. - - tol : float, optional - The tolerance for the gauge optimization. - - oob_check_interval : int, optional - If greater than zero, gauge transformations are allowed to fail (by raising - any exception) to indicate an out-of-bounds condition that the gauge optimizer - will avoid. If zero, then any gauge-transform failures just terminate the - optimization. - convert_model_to : str, dict, list, optional For use when `model` is an `ExplicitOpModel`. When not `None`, calls `model.convert_members_inplace(convert_model_to, set_default_gauge_group=False)` if @@ -217,35 +199,54 @@ def gaugeopt_to_target(model, target_model, *args, **kwargs): When not None, an MPI communicator for distributing the computation across multiple processors. - verbosity : int, optional - How much detail to send to stdout. + n_leak : int + Used in leakage modeling. If positive, this specifies how modify defintiions + of gate and SPAM metrics to reflect the fact that target gates do not have + well-defined actions outside the computational subspace. + """ - check_jac : bool - When True, check least squares analytic jacobian against finite differences + @staticmethod + def create_full_kwargs(args: tuple[Any,...], kwargs: dict[str, Any]): + full_kwargs = GaugeoptToTargetArgs.create_objective_passthrough_kwargs.copy() + full_kwargs.update(GaugeoptToTargetArgs.gaugeopt_custom_passthrough_kwargs) + full_kwargs.update(GaugeoptToTargetArgs.other_kwargs) + full_kwargs.update(kwargs) + + if extra_posargs := len(args) > 0: + msg = \ + f""" + Recieved {extra_posargs} positional arguments past `model` and `target_model`. + These should be passed as appropriate keyword arguments instead. This version + of pyGSTi will infer intended keyword arguments based on the legacy argument + positions. Future versions of pyGSTi will raise an error. + """ + _warnings.warn(msg) + for k,v in zip(GaugeoptToTargetArgs.old_trailing_positional_args, args): + full_kwargs[k] = v + + return full_kwargs + + +def gaugeopt_to_target(model, target_model, *args, **kwargs): + """ + Legacy function to optimize the gauge degrees of freedom of a model to that of a target. + + Use of more than two positional arguments is deprecated; keyword arguments should be used + instead. See the GaugeoptToTargetArgs class for available keyword arguments and their + default values. Returns ------- - model : if return_all == False + model : if kwargs.get('return_all', False) == False + + (goodnessMin, gaugeMx, model) : if kwargs.get('return_all', False) == True - (goodnessMin, gaugeMx, model) : if return_all == True Where goodnessMin is the minimum value of the goodness function (the best 'goodness') found, gaugeMx is the gauge matrix used to transform the model, and model is the final gauge-transformed model. """ - if extra_posargs := len(args) > 0: - msg = \ - f""" - Recieved {extra_posargs} positional arguments past `model` and `target_model`. - These should be passed as appropriate keyword arguments instead. This version - of pyGSTi will infer intended keyword arguments based on the legacy argument - positions. Future versions of pyGSTi will raise an error. - """ - _warnings.warn(msg) - kwargs.update({ k: v for (k, v) in zip(GaugeoptToTargetArgs.old_positional_args, args) }) - out = gaugeopt_to_target(model, target_model, **kwargs) - return out - full_kwargs = GaugeoptToTargetArgs.create_full_kwargs(kwargs) + full_kwargs = GaugeoptToTargetArgs.create_full_kwargs(args, kwargs) """ This function handles a strange situation where `target_model` can be None. From fe9c4ec53df881dc94a38ff6e2677f9ba86ae312 Mon Sep 17 00:00:00 2001 From: Riley Murray Date: Thu, 16 Oct 2025 21:19:34 -0700 Subject: [PATCH 18/18] extremely messy implementation of leakage metrics and gauge optimization that should work for arbitrary constructions of the computational subspace. Comitting to share code and get feedback before devising the clean-up and testing plan --- .../Examples/Leakage-automagic.ipynb | 227 ++++- pygsti/algorithms/gaugeopt.py | 799 ++++++++++-------- pygsti/baseobjs/basis.py | 40 + pygsti/models/explicitmodel.py | 6 + pygsti/models/model.py | 8 +- pygsti/protocols/gst.py | 18 +- pygsti/report/factory.py | 2 +- pygsti/report/reportables.py | 127 ++- pygsti/report/section/gauge.py | 12 +- pygsti/report/section/summary.py | 12 +- pygsti/tools/leakage.py | 213 +++-- pygsti/tools/matrixtools.py | 14 +- pygsti/tools/optools.py | 97 ++- pygsti/tools/sdptools.py | 19 +- 14 files changed, 1055 insertions(+), 539 deletions(-) diff --git a/jupyter_notebooks/Examples/Leakage-automagic.ipynb b/jupyter_notebooks/Examples/Leakage-automagic.ipynb index c586ff789..07f93c469 100644 --- a/jupyter_notebooks/Examples/Leakage-automagic.ipynb +++ b/jupyter_notebooks/Examples/Leakage-automagic.ipynb @@ -2,14 +2,21 @@ "cells": [ { "cell_type": "code", - "execution_count": null, + "execution_count": 1, "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", + "%load_ext autoreload\n", + "%autoreload 2\n", + "import pygsti\n", + "from pygsti.modelpacks import smq1Q_XYI as mp\n", + "from pygsti.tools.leakage import leaky_qubit_model_from_pspec, construct_leakage_report, add_lago_models, gate_leakage_profile\n", + "from pygsti.report import construct_standard_report\n", + "from pygsti.tools.optools import unitary_to_superop\n", "from pygsti.data import simulate_data\n", - "from pygsti.protocols import StandardGST, ProtocolData" + "from pygsti.protocols import StandardGST, ProtocolData, GateSetTomography, GSTGaugeOptSuite, GSTBadFitOptions\n", + "import numpy as np\n", + "import scipy.linalg as la" ] }, { @@ -23,30 +30,222 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 2, "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", + "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": 3, + "metadata": {}, + "outputs": [ + { + "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" + ] + } + ], + "source": [ + "ed = mp.create_gst_experiment_design(max_max_length=16)\n", + "ps = mp.processor_spec()\n", + "tm3 = leaky_qubit_model_from_pspec(ps, mx_basis='l2p1')\n", + "tm3.convert_members_inplace('CPTPLND')\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", + "se = 'multinomial'\n", + "leaky_gate = ('Gxpi2', 0)\n", + "dm3, leaking_state = with_leaky_gate(tm3, leaky_gate, strength=0.125)\n", + "if se == 'none':\n", + " import pygsti.objectivefns.objectivefns as pg_ofns\n", + " # Forward simulation can raise errors if there are positive\n", + " # probabilities < pg_ofns.DEFAULT_MIN_PROB_CLIP.\n", + " # We change that default here so forward simulation works with\n", + " # the ideal model.\n", + " pg_ofns.DEFAULT_MIN_PROB_CLIP = 1e-12\n", + "ds = simulate_data(dm3, ed.all_circuits_needing_data, num_samples=10_000, sample_error=se, seed=1997)" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "--- Iterative GST: [--------------------------------------------------] 0.0% 92 circuits ---\r" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/rjmurr/Documents/pygsti-leakage/hrl-repo/pygsti_repo/pygsti/objectivefns/objectivefns.py:4494: RuntimeWarning: divide by zero encountered in divide\n", + " p5over_lsvec = 0.5/lsvec\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "--- Iterative GST: [##########----------------------------------------] 20.0% 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: [####################------------------------------] 40.0% 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: [##############################--------------------] 60.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: [########################################----------] 80.0% 616 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% 616 circuits ---\n" + ] + } + ], + "source": [ + "gopsuite = GSTGaugeOptSuite.cast('stdgaugeopt')\n", + "gopsuite.gaugeopt_target = tm3\n", + "bfo = {'actions': ['wildcard1d'], 'wildcard1d_reference': 'diamond distance', 'threshold': 0.0}\n", + "bfo = GSTBadFitOptions.cast(bfo)\n", + "gst = GateSetTomography(initial_model=dm3, gaugeopt_suite=gopsuite, verbosity=1, name='CPTPLND', badfit_options=bfo)\n", "pd = ProtocolData(ed, ds)\n", "res = gst.run(pd)" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 5, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/rjmurr/Documents/pygsti-leakage/hrl-repo/pygsti_repo/pygsti/models/model.py:145: 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" + ] + } + ], + "source": [ + "import copy\n", + "resc = copy.deepcopy(res)\n", + "add_lago_models(resc, 'CPTPLND', verbosity=1)\n", + "est = resc.estimates['CPTPLND']" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "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 (158)! 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 (158)! 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", + "report_dir = f'example_files/leakage-report-automagic-251016-NleakNone-{se}'\n", + "report_object = construct_standard_report(resc, advanced_options={'n_leak': None}, title=f'fidelity-gop sample_error={se}')\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", diff --git a/pygsti/algorithms/gaugeopt.py b/pygsti/algorithms/gaugeopt.py index 80dd62ccd..0dfaecf5d 100644 --- a/pygsti/algorithms/gaugeopt.py +++ b/pygsti/algorithms/gaugeopt.py @@ -65,11 +65,11 @@ def parsed_model(model: Union[_OpModel, _ExplicitOpModel], convert_model_to: Opt return model_out @staticmethod - def parsed_method(target_model: Optional[_OpModel], method: str, gates_metric: str, spam_metric: str, n_leak: int) -> str: + def parsed_method(target_model: Optional[_OpModel], method: str, gates_metric: str, spam_metric: str, n_leak: Optional[int]) -> str: ls_mode_allowed = bool(target_model is not None and gates_metric.startswith("frobenius") and spam_metric.startswith("frobenius") - and n_leak == 0) + and isinstance(n_leak, int) and n_leak == -1) # and model.dim < 64: # least squares optimization seems uneffective if more than 3 qubits # -- observed by Lucas - should try to debug why 3 qubits seemed to cause trouble... if method == "ls" and not ls_mode_allowed: @@ -86,7 +86,7 @@ def parsed_method(target_model: Optional[_OpModel], method: str, gates_metric: s # used positional arguments past target_model. create_objective_passthrough_kwargs : dict[str,Any] = dict( - item_weights=None, + item_weights=None, # gets cast to a dict. cptp_penalty_factor=0, spam_penalty_factor=0, check_jac=False, @@ -148,7 +148,7 @@ def parsed_method(target_model: Optional[_OpModel], method: str, gates_metric: s return_all=False, # used in the function body (only for branching after passing to gaugeopt_custom) convert_model_to=None, # passthrough to parsed_model. comm=None, # passthrough to _create_objective_fn and gaugeopt_custom - n_leak=0 # passthrough to parsed_objective and _create_objective_fn + n_leak=-1 # passthrough to parsed_objective and _create_objective_fn; default value is a flag. ) """ gates_metric : {"frobenius", "fidelity", "tracedist"}, optional @@ -223,7 +223,10 @@ def create_full_kwargs(args: tuple[Any,...], kwargs: dict[str, Any]): _warnings.warn(msg) for k,v in zip(GaugeoptToTargetArgs.old_trailing_positional_args, args): full_kwargs[k] = v - + + if full_kwargs['item_weights'] is None: + full_kwargs['item_weights'] = dict() + return full_kwargs @@ -270,8 +273,6 @@ def gaugeopt_to_target(model, target_model, *args, **kwargs): # actual work: constructing objective_fn and jacobian_fn item_weights = full_kwargs['item_weights'] - if item_weights is None: - item_weights = {} cptp_penalty = full_kwargs['cptp_penalty_factor'] spam_penalty = full_kwargs['spam_penalty_factor'] comm = full_kwargs['comm'] @@ -475,15 +476,59 @@ def _call_jacobian_fn(gauge_group_el_vec): return newModel -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, n_leak=0) -> tuple[GGElObjective, GGElJacobian]: - """ - Creates the objective function and jacobian (if available) - for gaugeopt_to_target - """ - if item_weights is None: item_weights = {} +def _transform_with_oob_check(mdl, gauge_group_el, oob_check): + """ Helper function that sometimes checks if mdl.transform_inplace(gauge_group_el) fails. """ + mdl = mdl.copy() + if oob_check: + try: + mdl.transform_inplace(gauge_group_el) + except Exception as e: + raise ValueError("Out of bounds: %s" % str(e)) # signals OOB condition + else: + mdl.transform_inplace(gauge_group_el) + return mdl + + +def _gate_fidelity_targets(model, target_model): + from pygsti.report.reportables import eigenvalue_entanglement_infidelity + gate_fidelity_targets : dict[ _baseobjs.Label, tuple[_np.ndarray, Union[float, _np.floating]] ] = dict() + for lbl in target_model.operations: + G_target = target_model.operations[lbl].to_dense() + G_curest = model.operations[lbl].to_dense() + t = 1 - eigenvalue_entanglement_infidelity(G_curest, G_target, model.basis) + t = _np.clip(t, a_min=0.0, a_max=1.0) + gate_fidelity_targets[lbl] = (G_target, t) + return gate_fidelity_targets + + +def _prep_fidelity_targets(model, target_model): + prep_fidelity_targets : dict[ _baseobjs.Label, tuple[_np.ndarray, Union[float, _np.floating]] ] = dict() + for preplbl in target_model.preps: + rho_target = _tools.vec_to_stdmx( target_model.preps[preplbl].to_dense(), model.basis ) + rho_curest = _tools.vec_to_stdmx( model.preps[preplbl].to_dense(), model.basis ) + t = _tools.eigenvalue_fidelity(rho_curest, rho_target) + t = _np.clip(t, a_min=0.0, a_max=1.0) + prep_fidelity_targets[preplbl] = (rho_target, t) + return prep_fidelity_targets + + +def _povm_effect_fidelity_targets(model, target_model): + povm_fidelity_targets : dict[ _baseobjs.Label, tuple[dict, dict] ] = dict() + for povmlbl in target_model.povms: + M_target = target_model.povms[povmlbl] + M_curest = model.povms[povmlbl] + Es_target = {elbl : _tools.vec_to_stdmx(e_target.to_dense(), model.basis) for (elbl, e_target) in M_target.items() } + Es_curest = {elbl : _tools.vec_to_stdmx(e_target.to_dense(), model.basis) for (elbl, e_target) in M_curest.items() } + ts = {elbl : _tools.eigenvalue_fidelity(Es_target[elbl], Es_curest[elbl]) for elbl in M_target } + ts = {elbl : _np.clip(t, a_min=0.0, a_max=1.0) for (elbl, t) in ts.items()} + povm_fidelity_targets[povmlbl] = (Es_target, ts) + return povm_fidelity_targets + + +def _legacy_create_scalar_objective(model, target_model, + item_weights: dict[str,float], cptp_penalty_factor: float, spam_penalty_factor: float, + gates_metric: str, spam_metric: str, n_leak: Optional[int] + ) -> tuple[GGElObjective, GGElJacobian]: opWeight = item_weights.get('gates', 1.0) spamWeight = item_weights.get('spam', 1.0) prepWeight = item_weights.get('prep', spamWeight) @@ -495,393 +540,407 @@ def _create_objective_fn(model, target_model, item_weights: Optional[dict[str,fl if mxBasis.name == "unknown" and target_model is not None: mxBasis = target_model.basis - def _transform_with_oob_check(mdl, gauge_group_el, oob_check): - """ Helper function that sometimes checks if mdl.transform_inplace(gauge_group_el) fails. """ - mdl = mdl.copy() - if oob_check: - try: mdl.transform_inplace(gauge_group_el) - except Exception as e: - raise ValueError("Out of bounds: %s" % str(e)) # signals OOB condition - else: - mdl.transform_inplace(gauge_group_el) - return mdl - if method == "ls": - # least-squares case where objective function returns an array of - # the before-they're-squared difference terms and there's an analytic jacobian + udim = int(_np.sqrt(mxBasis.dim)) + I = _tools.matrixtools.IdentityOperator() + if n_leak is None and mxBasis.first_element_is_identity: + P = I + elif n_leak is None: + P = _tools.superop_subspace_projector(mxBasis) + elif n_leak > 0: + P = _tools.superop_subspace_projector(udim - n_leak, udim, mxBasis) + else: + # this also handles the case n_leak==-1, which is a flag. + P = I + n_leak = 0 + transform_mx_arg = (P, I) + # ^ The semantics of this tuple are defined by the frobeniusdist function + # in the ExplicitOpModelCalc class. + + gate_fidelity_targets : dict[ _baseobjs.Label, tuple[_np.ndarray, Union[float, _np.floating]] ] = dict() + if gates_metric == 'fidelity': + gate_fidelity_targets.update(_gate_fidelity_targets(model, target_model)) + + prep_fidelity_targets : dict[ _baseobjs.Label, tuple[_np.ndarray, Union[float, _np.floating]] ] = dict() + povm_fidelity_targets : dict[ _baseobjs.Label, tuple[dict, dict] ] = dict() + if spam_metric == 'fidelity': + prep_fidelity_targets.update(_prep_fidelity_targets(model, target_model)) + povm_fidelity_targets.update(_povm_effect_fidelity_targets(model, target_model)) + + def _objective_fn(gauge_group_el: _GaugeGroupElement, oob_check: bool) -> float: + mdl = _transform_with_oob_check(model, gauge_group_el, oob_check) + ret = 0 + + if cptp_penalty_factor > 0: + mdl.basis = mxBasis # set basis for jamiolkowski iso + cpPenaltyVec = _cptp_penalty(mdl, cptp_penalty_factor, mdl.basis) + ret += _np.sum(cpPenaltyVec) + + if spam_penalty_factor > 0: + mdl.basis = mxBasis + spamPenaltyVec = _spam_penalty(mdl, spam_penalty_factor, mdl.basis) + ret += _np.sum(spamPenaltyVec) + + 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 == 'commutators': + 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 * _np.linalg.norm((mop - top) @ P)**2 + + elif gates_metric == "fidelity": + # Leakage-aware metrics NOT available + for opLbl in mdl.operations: + wt = item_weights.get(opLbl, opWeight) + mop = mdl.operations[opLbl].to_dense() + top, t = gate_fidelity_targets[opLbl] + v = _tools.entanglement_fidelity(mop, top, mxBasis) + z = _np.abs(t - v) + ret += wt * z + + 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) - assert(gates_metric.startswith("frobenius") and spam_metric.startswith("frobenius")), \ - "Only 'frobenius' and 'frobeniustt' metrics can be used when `method='ls'`!" - assert(gates_metric == spam_metric) - frobenius_transform_target = bool(gates_metric == 'frobeniustt') # tt = "transform target" + else: + raise ValueError("Invalid gates_metric: %s" % gates_metric) - if frobenius_transform_target: - full_target_model = target_model.copy() - full_target_model.convert_members_inplace("full") # so we can gauge-transform the target model. + 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 == 'commutators': + for preplabel, m_prep in mdl.preps.items(): + wt_prep = 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) + com = rhoMx1 @ rhoMx2 - rhoMx2 - rhoMx1 + ret += wt_prep * _np.linalg.norm(com)**2 + + for povmlbl in target_model.povms: + wt_povm = item_weights.get(povmlbl, spamWeight) + M_curest = model.povms[povmlbl] + M_target = target_model.povms[povmlbl] + for elbl in M_target: + E_curest = _tools.vec_to_stdmx(M_curest[elbl].to_dense(), mxBasis) + E_target = _tools.vec_to_stdmx(M_target[elbl].to_dense(), mxBasis) + com_mat = E_curest @ E_target - E_target @ E_curest + com_vec = P @ _tools.stdmx_to_vec(com_mat, mxBasis) + val = _np.linalg.norm(com_vec) + ret += wt_povm * val**2 + + elif spam_metric == "fidelity": + # Leakage-aware metrics NOT available + val_prep = 0.0 + for preplbl in target_model.preps: + wt_prep = item_weights.get(preplbl, prepWeight) + rho_curest = _tools.vec_to_stdmx(model.preps[preplbl].to_dense(), mxBasis) + rho_target, t = prep_fidelity_targets[preplbl] + v = _tools.fidelity(rho_curest, rho_target) + val_prep += wt * abs(t - v) + val_povm = 0.0 + for povmlbl in target_model.povms: + wt_povm = item_weights.get(povmlbl, spamWeight) + M_target = target_model.povms[povmlbl] + M_curest = model.povms[povmlbl] + Es_target, ts = povm_fidelity_targets[povmlbl] + for elbl in M_target: + t_e = ts[elbl] + E = _tools.vec_to_stdmx(M_curest[elbl].to_dense(), mxBasis) + v_e = _tools.fidelity(E, Es_target[elbl]) + val_povm += wt * abs(t_e - v_e) + ret += (val_prep + val_povm) + + 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: - full_target_model = None # in case it get's referenced by mistake + raise ValueError("Invalid spam_metric: %s" % spam_metric) - def _objective_fn(gauge_group_el, oob_check): + return ret - if frobenius_transform_target: - transformed = _transform_with_oob_check(full_target_model, gauge_group_el.inverse(), oob_check) - other = model - else: - transformed = _transform_with_oob_check(model, gauge_group_el, oob_check) - other = target_model + return _objective_fn, None + + +def _legacy_create_least_squares_objective(model, target_model, + item_weights: dict[str,float], cptp_penalty_factor: float, spam_penalty_factor: float, + gates_metric: str, spam_metric: str, comm: Optional[Any], check_jac: bool + ) -> tuple[GGElObjective, GGElJacobian]: - residuals, _ = transformed.residuals(other, None, item_weights) + opWeight = item_weights.get('gates', 1.0) + spamWeight = item_weights.get('spam', 1.0) + mxBasis = model.basis - # We still the non-target model to be transformed and checked for these penalties - if cptp_penalty_factor > 0 or spam_penalty_factor > 0: - if frobenius_transform_target: - transformed = _transform_with_oob_check(model, gauge_group_el, oob_check) + #Use the target model's basis if model's is unknown + # (as it can often be if it's just come from an logl opt, + # since from_vector will clear any basis info) + if mxBasis.name == "unknown" and target_model is not None: + mxBasis = target_model.basis - if cptp_penalty_factor > 0: - transformed.basis = mxBasis - cpPenaltyVec = _cptp_penalty(transformed, cptp_penalty_factor, transformed.basis) - else: cpPenaltyVec = [] # so concatenate ignores + assert(gates_metric.startswith("frobenius") and spam_metric.startswith("frobenius")), \ + "Only 'frobenius' and 'frobeniustt' metrics can be used when `method='ls'`!" + assert(gates_metric == spam_metric) + frobenius_transform_target = bool(gates_metric == 'frobeniustt') # tt = "transform target" - if spam_penalty_factor > 0: - transformed.basis = mxBasis - spamPenaltyVec = _spam_penalty(transformed, spam_penalty_factor, transformed.basis) - else: spamPenaltyVec = [] # so concatenate ignores + if frobenius_transform_target: + full_target_model = target_model.copy() + full_target_model.convert_members_inplace("full") # so we can gauge-transform the target model. + else: + full_target_model = None # in case it get's referenced by mistake - return _np.concatenate((residuals, cpPenaltyVec, spamPenaltyVec)) - else: - return residuals + def _objective_fn(gauge_group_el: _GaugeGroupElement, oob_check: bool) -> _np.ndarray: - def _jacobian_fn(gauge_group_el): + if frobenius_transform_target: + transformed = _transform_with_oob_check(full_target_model, gauge_group_el.inverse(), oob_check) + other = model + else: + transformed = _transform_with_oob_check(model, gauge_group_el, oob_check) + other = target_model - #Penalty terms below always act on the transformed non-target model. - original_gauge_group_el = gauge_group_el + residuals, _ = transformed.residuals(other, None, item_weights) + # We still the non-target model to be transformed and checked for these penalties + if cptp_penalty_factor > 0 or spam_penalty_factor > 0: if frobenius_transform_target: - gauge_group_el = gauge_group_el.inverse() - mdl_pre = full_target_model.copy() - mdl_post = mdl_pre.copy() - else: - mdl_pre = model.copy() - mdl_post = mdl_pre.copy() - mdl_post.transform_inplace(gauge_group_el) - - # Indices: Jacobian output matrix has shape (L, N) - start = 0 - d = mdl_pre.dim - N = gauge_group_el.num_params - L = mdl_pre.num_elements - - #Compute "extra" (i.e. beyond the model-element) rows of jacobian - if cptp_penalty_factor != 0: L += _cptp_penalty_size(mdl_pre) - if spam_penalty_factor != 0: L += _spam_penalty_size(mdl_pre) - - #Set basis for pentaly term calculation - if cptp_penalty_factor != 0 or spam_penalty_factor != 0: - mdl_pre.basis = mxBasis - mdl_post.basis = mxBasis - - jacMx = _np.zeros((L, N)) - - #Overview of terms: - # objective: op_term = (S_inv * gate * S - target_op) - # jac: d(op_term) = (d (S_inv) * gate * S + S_inv * gate * dS ) - # d(op_term) = (-(S_inv * dS * S_inv) * gate * S + S_inv * gate * dS ) - - # objective: rho_term = (S_inv * rho - target_rho) - # jac: d(rho_term) = d (S_inv) * rho - # d(rho_term) = -(S_inv * dS * S_inv) * rho - - # objective: ET_term = (E.T * S - target_E.T) - # jac: d(ET_term) = E.T * dS - - #Overview of terms when frobenius_transform_target == True). Note that the objective - #expressions are identical to the above except for an additional overall minus sign and S <=> S_inv. - - # objective: op_term = (gate - S * target_op * S_inv) - # jac: d(op_term) = -(dS * target_op * S_inv + S * target_op * -(S_inv * dS * S_inv) ) - # d(op_term) = (-dS * target_op * S_inv + S * target_op * (S_inv * dS * S_inv) ) - - # objective: rho_term = (rho - S * target_rho) - # jac: d(rho_term) = - dS * target_rho - - # objective: ET_term = (E.T - target_E.T * S_inv) - # jac: d(ET_term) = - target_E.T * -(S_inv * dS * S_inv) - # d(ET_term) = target_E.T * (S_inv * dS * S_inv) - - #Distribute computation across processors - allDerivColSlice = slice(0, N) - derivSlices, myDerivColSlice, derivOwners, mySubComm = \ - _mpit.distribute_slice(allDerivColSlice, comm) - if mySubComm is not None: - _warnings.warn("Note: more CPUs(%d)" % comm.Get_size() - + " than gauge-opt derivative columns(%d)!" % N) # pragma: no cover - - n = _slct.length(myDerivColSlice) - wrtIndices = _slct.indices(myDerivColSlice) if (n < N) else None - my_jacMx = jacMx[:, myDerivColSlice] # just the columns I'm responsible for - - # S, and S_inv are shape (d,d) - #S = gauge_group_el.transform_matrix - S_inv = gauge_group_el.transform_matrix_inverse - dS = gauge_group_el.deriv_wrt_params(wrtIndices) # shape (d*d),n - dS.shape = (d, d, n) # call it (d1,d2,n) - dS = _np.rollaxis(dS, 2) # shape (n, d1, d2) - assert(dS.shape == (n, d, d)) - - # --- NOTE: ordering here, with running `start` index MUST - # correspond to those in Model.residuals, which in turn - # must correspond to those in ForwardSimulator.residuals - which - # currently orders as: gates, simplified_ops, preps, effects. - - # -- LinearOperator terms - # ------------------------- - for lbl, G in mdl_pre.operations.items(): - # d(op_term) = S_inv * (-dS * S_inv * G * S + G * dS) = S_inv * (-dS * G' + G * dS) - # Note: (S_inv * G * S) is G' (transformed G) - wt = item_weights.get(lbl, opWeight) - left = -1 * _np.dot(dS, mdl_post.operations[lbl].to_dense('minimal')) # shape (n,d1,d2) - right = _np.swapaxes(_np.dot(G.to_dense('minimal'), dS), 0, 1) # shape (d1,n,d2) -> (n,d1,d2) - result = _np.swapaxes(_np.dot(S_inv, left + right), 1, 2) # shape (d1, d2, n) - result = result.reshape((d**2, n)) # must copy b/c non-contiguous - my_jacMx[start:start + d**2] = wt * result - start += d**2 + transformed = _transform_with_oob_check(model, gauge_group_el, oob_check) - # -- Instrument terms - # ------------------------- - for ilbl, Inst in mdl_pre.instruments.items(): - wt = item_weights.get(ilbl, opWeight) - for lbl, G in Inst.items(): - # same calculation as for operation terms - left = -1 * _np.dot(dS, mdl_post.instruments[ilbl][lbl].to_dense('minimal')) # (n,d1,d2) - right = _np.swapaxes(_np.dot(G.to_dense('minimal'), dS), 0, 1) # (d1,n,d2) -> (n,d1,d2) - result = _np.swapaxes(_np.dot(S_inv, left + right), 1, 2) # shape (d1, d2, n) - result = result.reshape((d**2, n)) # must copy b/c non-contiguous - my_jacMx[start:start + d**2] = wt * result - start += d**2 - - # -- prep terms - # ------------------------- - for lbl, rho in mdl_post.preps.items(): - # d(rho_term) = -(S_inv * dS * S_inv) * rho - # Note: (S_inv * rho) is transformed rho - wt = item_weights.get(lbl, spamWeight) - Sinv_dS = _np.dot(S_inv, dS) # shape (d1,n,d2) - result = -1 * _np.dot(Sinv_dS, rho.to_dense('minimal')) # shape (d,n) - my_jacMx[start:start + d] = wt * result - start += d + if cptp_penalty_factor > 0: + transformed.basis = mxBasis + cpPenaltyVec = _cptp_penalty(transformed, cptp_penalty_factor, transformed.basis) + else: cpPenaltyVec = [] # so concatenate ignores - # -- effect terms - # ------------------------- - for povmlbl, povm in mdl_pre.povms.items(): - for lbl, E in povm.items(): - # d(ET_term) = E.T * dS - wt = item_weights.get(povmlbl + "_" + lbl, spamWeight) - result = _np.dot(E.to_dense('minimal')[None, :], dS).T # shape (1,n,d2).T => (d2,n,1) - my_jacMx[start:start + d] = wt * result.squeeze(2) # (d2,n) - start += d - - # -- penalty terms -- Note: still use original gauge transform applied to `model` - # ------------------------- - if cptp_penalty_factor > 0 or spam_penalty_factor > 0: - if frobenius_transform_target: # reset back to non-target-tranform "mode" - gauge_group_el = original_gauge_group_el - mdl_pre = model.copy() - mdl_post = mdl_pre.copy() - mdl_post.transform_inplace(gauge_group_el) - - if cptp_penalty_factor > 0: - start += _cptp_penalty_jac_fill(my_jacMx[start:], mdl_pre, mdl_post, - gauge_group_el, cptp_penalty_factor, - mdl_pre.basis, wrtIndices) - - if spam_penalty_factor > 0: - start += _spam_penalty_jac_fill(my_jacMx[start:], mdl_pre, mdl_post, - gauge_group_el, spam_penalty_factor, - mdl_pre.basis, wrtIndices) - - #At this point, each proc has filled the portions (columns) of jacMx that - # it's responsible for, and so now we gather them together. - _mpit.gather_slices(derivSlices, derivOwners, jacMx, [], 1, comm) - #Note jacMx is completely filled (on all procs) - - if check_jac and (comm is None or comm.Get_rank() == 0): - def _mock_objective_fn(v): - return _objective_fn(gauge_group_el, False) - vec = gauge_group_el.to_vector() - _opt.check_jac(_mock_objective_fn, vec, jacMx, tol=1e-5, eps=1e-9, err_type='abs', - verbosity=1) - - return jacMx + if spam_penalty_factor > 0: + transformed.basis = mxBasis + spamPenaltyVec = _spam_penalty(transformed, spam_penalty_factor, transformed.basis) + else: spamPenaltyVec = [] # so concatenate ignores - else: - # 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: - P = _tools.superop_subspace_projector(dim - n_leak, dim, mxBasis) - transform_mx_arg = (P, _tools.matrixtools.IdentityOperator()) - # ^ The semantics of this tuple are defined by the frobeniusdist function - # in the ExplicitOpModelCalc class. + return _np.concatenate((residuals, cpPenaltyVec, spamPenaltyVec)) else: - transform_mx_arg = None - # ^ It would be equivalent to set this to a pair of IdentityOperator objects. - from pygsti.report.reportables import eigenvalue_entanglement_infidelity, vec_fidelity + return residuals + def _jacobian_fn(gauge_group_el: _GaugeGroupElement) -> _np.ndarray: - """ - We say that a (model, target) pair admit a _perfect gauge_ if "model" has no relational errors - when compared to "target." If "model" is considered in the perfect gauge, then - - (1) its fidelities with the target gates match the eigenvalue fidelities with the target gates; - - (2) the (rho_e, M_t) SPAM probs will match the (rho_e, M_e) SPAM probs + #Penalty terms below always act on the transformed non-target model. + original_gauge_group_el = gauge_group_el - (3) the (rho_t, M_e) SPAM probs will match the (rho_e, M_e) SPAM probs - - ^ (3) originally said "fidelities ... match those between rho_e and {E_t : E_t in M_target}", which is WRONG. - - - In view of this, taking fidelity as the gauge-optimization objective tries to minimize an - aggregration of any mismatch between the various fidelities above. - """ - - gate_fidelity_targets : dict[ Union[str, _baseobjs.Label], Union[float, _np.floating] ] = dict() - if gates_metric == 'fidelity': - for lbl in target_model.operations: - G_target = target_model.operations[lbl] - G_curest = model.operations[lbl] - t = 1 - eigenvalue_entanglement_infidelity(G_curest.to_dense(), G_target.to_dense(), model.basis) - gate_fidelity_targets[lbl] = min(t, 1.0) - - spam_fidelity_targets : dict[ - tuple[Union[str, _baseobjs.Label], Union[str, _baseobjs.Label]], - dict[Union[str, _baseobjs.Label], Union[float, _np.floating]] - ] = dict() - if spam_metric == 'fidelity': - for preplbl in target_model.preps: - for povmlbl in target_model.povms: - rho_curest = model.preps[preplbl].to_dense() - M_curest = model.povms[povmlbl] - t = {elbl: _np.vdot(rho_curest, e.to_dense()) for (elbl, e) in M_curest.items() } - spam_fidelity_targets[(preplbl, povmlbl)] = t + if frobenius_transform_target: + gauge_group_el = gauge_group_el.inverse() + mdl_pre = full_target_model.copy() + mdl_post = mdl_pre.copy() + else: + mdl_pre = model.copy() + mdl_post = mdl_pre.copy() + mdl_post.transform_inplace(gauge_group_el) + + # Indices: Jacobian output matrix has shape (L, N) + start = 0 + d = mdl_pre.dim + N = gauge_group_el.num_params + L = mdl_pre.num_elements + + #Compute "extra" (i.e. beyond the model-element) rows of jacobian + if cptp_penalty_factor != 0: L += _cptp_penalty_size(mdl_pre) + if spam_penalty_factor != 0: L += _spam_penalty_size(mdl_pre) + + #Set basis for pentaly term calculation + if cptp_penalty_factor != 0 or spam_penalty_factor != 0: + mdl_pre.basis = mxBasis + mdl_post.basis = mxBasis + + jacMx = _np.zeros((L, N)) + + #Overview of terms: + # objective: op_term = (S_inv * gate * S - target_op) + # jac: d(op_term) = (d (S_inv) * gate * S + S_inv * gate * dS ) + # d(op_term) = (-(S_inv * dS * S_inv) * gate * S + S_inv * gate * dS ) + + # objective: rho_term = (S_inv * rho - target_rho) + # jac: d(rho_term) = d (S_inv) * rho + # d(rho_term) = -(S_inv * dS * S_inv) * rho + + # objective: ET_term = (E.T * S - target_E.T) + # jac: d(ET_term) = E.T * dS + + #Overview of terms when frobenius_transform_target == True). Note that the objective + #expressions are identical to the above except for an additional overall minus sign and S <=> S_inv. + + # objective: op_term = (gate - S * target_op * S_inv) + # jac: d(op_term) = -(dS * target_op * S_inv + S * target_op * -(S_inv * dS * S_inv) ) + # d(op_term) = (-dS * target_op * S_inv + S * target_op * (S_inv * dS * S_inv) ) + + # objective: rho_term = (rho - S * target_rho) + # jac: d(rho_term) = - dS * target_rho + + # objective: ET_term = (E.T - target_E.T * S_inv) + # jac: d(ET_term) = - target_E.T * -(S_inv * dS * S_inv) + # d(ET_term) = target_E.T * (S_inv * dS * S_inv) + + #Distribute computation across processors + allDerivColSlice = slice(0, N) + derivSlices, myDerivColSlice, derivOwners, mySubComm = \ + _mpit.distribute_slice(allDerivColSlice, comm) + if mySubComm is not None: + _warnings.warn("Note: more CPUs(%d)" % comm.Get_size() + + " than gauge-opt derivative columns(%d)!" % N) # pragma: no cover + + n = _slct.length(myDerivColSlice) + wrtIndices = _slct.indices(myDerivColSlice) if (n < N) else None + my_jacMx = jacMx[:, myDerivColSlice] # just the columns I'm responsible for + + # S, and S_inv are shape (d,d) + #S = gauge_group_el.transform_matrix + S_inv = gauge_group_el.transform_matrix_inverse + dS = gauge_group_el.deriv_wrt_params(wrtIndices) # shape (d*d),n + dS.shape = (d, d, n) # call it (d1,d2,n) + dS = _np.rollaxis(dS, 2) # shape (n, d1, d2) + assert(dS.shape == (n, d, d)) + + # --- NOTE: ordering here, with running `start` index MUST + # correspond to those in Model.residuals, which in turn + # must correspond to those in ForwardSimulator.residuals - which + # currently orders as: gates, simplified_ops, preps, effects. + + # -- LinearOperator terms + # ------------------------- + for lbl, G in mdl_pre.operations.items(): + # d(op_term) = S_inv * (-dS * S_inv * G * S + G * dS) = S_inv * (-dS * G' + G * dS) + # Note: (S_inv * G * S) is G' (transformed G) + wt = item_weights.get(lbl, opWeight) + left = -1 * _np.dot(dS, mdl_post.operations[lbl].to_dense('minimal')) # shape (n,d1,d2) + right = _np.swapaxes(_np.dot(G.to_dense('minimal'), dS), 0, 1) # shape (d1,n,d2) -> (n,d1,d2) + result = _np.swapaxes(_np.dot(S_inv, left + right), 1, 2) # shape (d1, d2, n) + result = result.reshape((d**2, n)) # must copy b/c non-contiguous + my_jacMx[start:start + d**2] = wt * result + start += d**2 + + # -- Instrument terms + # ------------------------- + for ilbl, Inst in mdl_pre.instruments.items(): + wt = item_weights.get(ilbl, opWeight) + for lbl, G in Inst.items(): + # same calculation as for operation terms + left = -1 * _np.dot(dS, mdl_post.instruments[ilbl][lbl].to_dense('minimal')) # (n,d1,d2) + right = _np.swapaxes(_np.dot(G.to_dense('minimal'), dS), 0, 1) # (d1,n,d2) -> (n,d1,d2) + result = _np.swapaxes(_np.dot(S_inv, left + right), 1, 2) # shape (d1, d2, n) + result = result.reshape((d**2, n)) # must copy b/c non-contiguous + my_jacMx[start:start + d**2] = wt * result + start += d**2 - def _objective_fn(gauge_group_el, oob_check): - mdl = _transform_with_oob_check(model, gauge_group_el, oob_check) - ret = 0 + # -- prep terms + # ------------------------- + for lbl, rho in mdl_post.preps.items(): + # d(rho_term) = -(S_inv * dS * S_inv) * rho + # Note: (S_inv * rho) is transformed rho + wt = item_weights.get(lbl, spamWeight) + Sinv_dS = _np.dot(S_inv, dS) # shape (d1,n,d2) + result = -1 * _np.dot(Sinv_dS, rho.to_dense('minimal')) # shape (d,n) + my_jacMx[start:start + d] = wt * result + start += d + + # -- effect terms + # ------------------------- + for povmlbl, povm in mdl_pre.povms.items(): + for lbl, E in povm.items(): + # d(ET_term) = E.T * dS + wt = item_weights.get(povmlbl + "_" + lbl, spamWeight) + result = _np.dot(E.to_dense('minimal')[None, :], dS).T # shape (1,n,d2).T => (d2,n,1) + my_jacMx[start:start + d] = wt * result.squeeze(2) # (d2,n) + start += d + + # -- penalty terms -- Note: still use original gauge transform applied to `model` + # ------------------------- + if cptp_penalty_factor > 0 or spam_penalty_factor > 0: + if frobenius_transform_target: # reset back to non-target-tranform "mode" + gauge_group_el = original_gauge_group_el + mdl_pre = model.copy() + mdl_post = mdl_pre.copy() + mdl_post.transform_inplace(gauge_group_el) if cptp_penalty_factor > 0: - mdl.basis = mxBasis # set basis for jamiolkowski iso - cpPenaltyVec = _cptp_penalty(mdl, cptp_penalty_factor, mdl.basis) - ret += _np.sum(cpPenaltyVec) + start += _cptp_penalty_jac_fill(my_jacMx[start:], mdl_pre, mdl_post, + gauge_group_el, cptp_penalty_factor, + mdl_pre.basis, wrtIndices) if spam_penalty_factor > 0: - mdl.basis = mxBasis - spamPenaltyVec = _spam_penalty(mdl, spam_penalty_factor, mdl.basis) - ret += _np.sum(spamPenaltyVec) - - 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": - # Leakage-aware metrics NOT available - for opLbl in mdl.operations: - wt = item_weights.get(opLbl, opWeight) - top = target_model.operations[opLbl].to_dense() - mop = mdl.operations[opLbl].to_dense() - t = gate_fidelity_targets[opLbl] - v = _tools.entanglement_fidelity(top, mop, mxBasis) - z = _np.abs(t - v) - ret += wt * z - - 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) + start += _spam_penalty_jac_fill(my_jacMx[start:], mdl_pre, mdl_post, + gauge_group_el, spam_penalty_factor, + mdl_pre.basis, wrtIndices) - if "frobenius" in spam_metric and gates_metric == spam_metric: - # We already handled SPAM error in this case. Just return. - return ret + #At this point, each proc has filled the portions (columns) of jacMx that + # it's responsible for, and so now we gather them together. + _mpit.gather_slices(derivSlices, derivOwners, jacMx, [], 1, comm) + #Note jacMx is completely filled (on all procs) - 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 - val = 0.0 - for preplbl in target_model.preps: - wt_prep = item_weights.get(preplbl, prepWeight) - for povmlbl in target_model.povms: - wt_povm = item_weights.get(povmlbl, spamWeight) - rho_curest = model.preps[preplbl].to_dense() - rho_target = target_model.preps[preplbl].to_dense() - M_curest = model.povms[povmlbl] - M_target = target_model.povms[povmlbl] - vs_prep = {elbl: _np.vdot(rho_curest, e_target.to_dense()) for (elbl, e_target) in M_target.items() } - vs_povm = {elbl: _np.vdot(rho_target, e_curest.to_dense()) for (elbl, e_curest) in M_curest.items() } - t_dict = spam_fidelity_targets[(preplbl, povmlbl)] - val1 = 0.0 - for lbl, f in vs_prep.items(): - val1 += _np.abs(t_dict[lbl] - f) - val2 = 0.0 - for lbl, f in vs_povm.items(): - val2 += _np.abs(t_dict[lbl] - f) - val += (wt_prep * val1 + wt_povm * val2) - ret += val - - 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) + if check_jac and (comm is None or comm.Get_rank() == 0): + def _mock_objective_fn(v): + return _objective_fn(gauge_group_el, False) + vec = gauge_group_el.to_vector() + _opt.check_jac(_mock_objective_fn, vec, jacMx, tol=1e-5, eps=1e-9, err_type='abs', + verbosity=1) - else: - raise ValueError("Invalid spam_metric: %s" % spam_metric) + return jacMx - return ret + return _objective_fn, _jacobian_fn - _jacobian_fn = None - return _objective_fn, _jacobian_fn +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, n_leak: Optional[int]=None) -> tuple[GGElObjective, GGElJacobian]: + if item_weights is None: + item_weights = dict() + if method == 'ls': + return _legacy_create_least_squares_objective( + model, target_model, item_weights, cptp_penalty_factor, spam_penalty_factor, gates_metric, + spam_metric, comm, check_jac + ) + else: + return _legacy_create_scalar_objective( + model, target_model, item_weights, cptp_penalty_factor, spam_penalty_factor, gates_metric, + spam_metric, n_leak + ) def _cptp_penalty_size(mdl): diff --git a/pygsti/baseobjs/basis.py b/pygsti/baseobjs/basis.py index 9e5734a1e..8d13db56b 100644 --- a/pygsti/baseobjs/basis.py +++ b/pygsti/baseobjs/basis.py @@ -19,6 +19,8 @@ from typing import ( Union, Optional, + Any, + Sequence ) import numpy as _np @@ -263,6 +265,8 @@ def __init__(self, name: str, longname: str, real: bool, sparse: bool): self.longname = longname self.real = real # whether coefficients must be real (*not* whether elements are real - they're always complex) self.sparse = sparse # whether elements are stored as sparse vectors/matrices + self._is_hermitian = None + self._implies_leakage = None @property def dim(self): @@ -317,6 +321,31 @@ def elsize(self): if self.elshape is None: return 0 return int(_np.prod(self.elshape)) + # @property + # def ellookup(self) -> dict[str, Any]: + # raise NotImplementedError() + + # @property + # def elements(self) -> Sequence[Any]: + # raise NotImplementedError() + + # @property.setter + # def elements(self, val) -> None: + # self. + + # @property + # def labels(self) -> Sequence[Any]: + # raise ValueError() + + @property + def implies_leakage_modeling(self) -> bool: + if (not hasattr(self, '_implies_leakage')) or self._implies_leakage is None: + if 'I' not in self.labels: + return False + I = self.ellookup['I'] + self._implies_leakage = _np.linalg.matrix_rank(I) < I.shape[0] + return self._implies_leakage + @property def first_element_is_identity(self): """ @@ -356,6 +385,17 @@ def is_partial(self) -> bool: """ return not self.is_complete() + def is_hermitian(self) -> bool: + from pygsti.tools.matrixtools import is_hermitian as matrix_is_hermitian + if self.elndim != 2 or self.elshape[0] != self.elshape[1]: + return False + if not hasattr(self, '_is_hermitian') or self._is_hermitian is None: + tol = self.dim * _np.finfo(self.elements[0].dtype).eps + is_hermitian = all([matrix_is_hermitian(el, tol=tol) for el in self.elements]) + self._is_hermitian = is_hermitian + return self._is_hermitian + + @property def vector_elements(self): """ diff --git a/pygsti/models/explicitmodel.py b/pygsti/models/explicitmodel.py index 5d276d9e8..7e3ae7c30 100644 --- a/pygsti/models/explicitmodel.py +++ b/pygsti/models/explicitmodel.py @@ -10,6 +10,8 @@ # http://www.apache.org/licenses/LICENSE-2.0 or in the LICENSE file in the root pyGSTi directory. #*************************************************************************************************** +from __future__ import annotations + import collections as _collections import itertools as _itertools import uuid as _uuid @@ -1635,6 +1637,10 @@ def add_availability(opkey, op): availability, instrument_names=list(self.instruments.keys()), nonstd_instruments=self.instruments) + def copy(self) -> ExplicitOpModel: + c = super().copy() # <-- that is indeed an ExplicitOpModel + return c # type: ignore + def create_modelmember_graph(self): return _MMGraph({ 'preps': self.preps, diff --git a/pygsti/models/model.py b/pygsti/models/model.py index 5f6f8e940..27dba1742 100644 --- a/pygsti/models/model.py +++ b/pygsti/models/model.py @@ -10,6 +10,8 @@ # http://www.apache.org/licenses/LICENSE-2.0 or in the LICENSE file in the root pyGSTi directory. #*************************************************************************************************** +from __future__ import annotations + import bisect as _bisect import copy as _copy import itertools as _itertools @@ -2309,7 +2311,7 @@ def _post_copy(self, copy_into, memo): copy_into._reinit_opcaches() super(OpModel, self)._post_copy(copy_into, memo) - def copy(self): + def copy(self) -> OpModel: """ Copy this model. @@ -2320,13 +2322,13 @@ def copy(self): """ self._clean_paramvec() # ensure _paramvec is rebuilt if needed if OpModel._pcheck: self._check_paramvec() - ret = Model.copy(self) + ret = Model.copy(self) # <-- don't be fooled! That's an OpModel. if self._param_bounds is not None and self.parameter_labels is not None: ret._clean_paramvec() # will *always* rebuild paramvec; do now so we can preserve param bounds assert _np.all(self.parameter_labels == ret.parameter_labels) # ensure ordering is the same ret._param_bounds = self._param_bounds.copy() if OpModel._pcheck: ret._check_paramvec() - return ret + return ret # type: ignore def create_modelmember_graph(self): """ diff --git a/pygsti/protocols/gst.py b/pygsti/protocols/gst.py index ea169e1a7..283ca42b9 100644 --- a/pygsti/protocols/gst.py +++ b/pygsti/protocols/gst.py @@ -2424,17 +2424,21 @@ def _memberdicts(_m) -> tuple[dict,dict,dict,dict]: for lbl, gaugeopt_model in gopped_models.items(): argdicts = gaugeopt_suite.gaugeopt_argument_dicts.get(lbl, dict()) - n_leak = 0 + n_leak: Optional[int] = None if isinstance(argdicts, list) and len(argdicts) > 0: n_leak = argdicts[0].get('n_leak', n_leak) - if n_leak > 0: - dim = gaugeopt_model.basis.state_space.udim - assert n_leak == 1 - assert dim == 3 - P = _tools.superop_subspace_projector(2, 3, gaugeopt_model.basis) + basis = gaugeopt_model.basis + udim = int(_np.round(_np.sqrt(basis.dim))) + I = _tools.matrixtools.IdentityOperator() + if n_leak is None and basis.first_element_is_identity: + P = I + elif n_leak is None: + P = _tools.superop_subspace_projector(basis) + elif n_leak > 0: + P = _tools.superop_subspace_projector(udim - n_leak, udim, basis) else: - P = _tools.matrixtools.IdentityOperator() + P = I ops, preps, _, insts = _memberdicts(gaugeopt_model) mx_basis = gaugeopt_model.basis diff --git a/pygsti/report/factory.py b/pygsti/report/factory.py index c8781f53f..aae5827bd 100644 --- a/pygsti/report/factory.py +++ b/pygsti/report/factory.py @@ -1230,7 +1230,7 @@ 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) + n_leak = advanced_options.get('n_leak', None) # ^ 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) diff --git a/pygsti/report/reportables.py b/pygsti/report/reportables.py index 949379388..5111059ea 100644 --- a/pygsti/report/reportables.py +++ b/pygsti/report/reportables.py @@ -17,6 +17,7 @@ """ import importlib import warnings as _warnings +from typing import Union import numpy as _np import scipy.linalg as _spl @@ -40,6 +41,8 @@ FINITE_DIFF_EPS = 1e-7 +BasisLike = Union[str, _Basis] + def _null_fn(*arg): return None @@ -1020,7 +1023,13 @@ def maximum_trace_dist(gate, mx_basis): 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 + # .... ^ What does that function call do ??? + if not isinstance(mx_basis, _Basis): + mx_basis = _Basis.cast(mx_basis, dim=gate.shape[0]) + if mx_basis.implies_leakage_modeling: + n_leak = None # we'll infer a detailed leakage model from mx_basis + else: + n_leak = 1 # assume one leakage level return _tools.subspace_jtracedist(gate, closestUOpMx, mx_basis, n_leak) @@ -1168,7 +1177,12 @@ def entanglement_fidelity(a, b, mx_basis): def subspace_entanglement_fidelity(a, b, mx_basis): - n_leak = 1 + if not isinstance(mx_basis, _Basis): + mx_basis = _Basis.cast(mx_basis, dim=a.shape[0]) + if mx_basis.implies_leakage_modeling: + n_leak = None # we'll infer a detailed leakage model from mx_basis + else: + n_leak = 1 # assume one leakage level return _tools.subspace_entanglement_fidelity(a, b, mx_basis, n_leak) @@ -1274,7 +1288,12 @@ def frobenius_diff(a, b, mx_basis): # assume vary model1, model2 fixed def leaky_gate_frob_dist(a, b, mx_basis): - n_leak = 1 + if not isinstance(mx_basis, _Basis): + mx_basis = _Basis.cast(mx_basis, dim=a.shape[0]) + if mx_basis.implies_leakage_modeling: + n_leak = None # we'll infer a detailed leakage model from mx_basis + else: + n_leak = 1 # assume one leakage level return _tools.subspace_superop_fro_dist(a, b, mx_basis, n_leak) @@ -1308,7 +1327,12 @@ def jtrace_diff(a, b, mx_basis): # assume vary model1, model2 fixed def leaky_jtrace_diff(a, b, mx_basis): - n_leak = 1 + if not isinstance(mx_basis, _Basis): + mx_basis = _Basis.cast(mx_basis, dim=a.shape[0]) + if mx_basis.implies_leakage_modeling: + n_leak = None # we'll infer a detailed leakage model from mx_basis + else: + n_leak = 1 # assume one leakage level return _tools.subspace_jtracedist(a, b, mx_basis, n_leak) @@ -1580,44 +1604,74 @@ def eigenvalue_nonunitary_avg_gate_infidelity(a, b, mx_basis): # init args == (model1, model2, op_label) -def eigenvalue_entanglement_infidelity(a, b, mx_basis): +def eigenvalue_entanglement_infidelity(a: _np.ndarray, b: _np.ndarray, mx_basis: BasisLike, is_tp=None, is_unitary=None, tol=1e-8) -> _np.floating: """ - Eigenvalue entanglement infidelity between a and b + Eigenvalue entanglement infidelity is the infidelity between certain diagonal matrices that + contain [see Note 1] the eigenvalues of J(a) and J(b), where J(.) is the Jamiolkowski + isomorphism map. This is equivalent to computing infidelity of (J(a), J(b)) while pretending + that they share an eigenbasis. Parameters ---------- - a : numpy.ndarray - The first process (transfer) matrix. + is_tp : bool, optional (default None) + Flag indicating that a and b are TP in their provided basis. If None (the default), + an explicit check is performed up to numerical tolerance `tol`. If True/False, then + the check is skipped and this is used as the result of the check as though it were + performed. - b : numpy.ndarray - The second process (transfer) matrix. + is_unitary : bool, optional (default None) + Flag indicating that b is unitary. If None (the default) an explicit check is performed + up to tolerance `tol`. If True/False, then the check is skipped and this is used as the + result of the check as though it were performed. - mx_basis : Basis or {'pp', 'gm', 'std'} - the basis that `a` and `b` are in. + Notes + ----- + [1] The eigenvalues are ordered using a min-weight matching algorithm. - Returns - ------- - float - """ - d2 = a.shape[0] - evA = _np.linalg.eigvals(a) - evB = _np.linalg.eigvals(b) - """ - IDEA: check if (a,b) are normal operators. If they are, then proceed by - computing a Schur decomposition to get orthonormal eigenbases. From there, - match them based on a dissimilarity kernel of abs(1 - ), rather than - the dissimilarity kernel abs(ev_a - ev_b). + [2] If a and b are trace-preserving (TP) and b is unitary, then we compute eigenvalue + entanglement fidelity by leveraging three facts: + + 1. If (x,y) share an eigenbasis and have (consistently ordered) eigenvalues (v(x),v(y)), + then Tr(x y^{\\dagger}) is the inner product of v(x) and v(y). + + 2. J(a) and J(b) share an eigenbasis of if `a` and `b` share an eigenbasis. + + 3. If `a` and `b` are TP and `b` is unitary, then their entanglement fidelity can be + expressed as |Tr( a.b^{\\dagger} )| / d^2. + + Then together, we see that if (a,b) are TP and `b` is unitary, then their eigenvalue + entanglement fidelity is expressible as the inner product of v(a) and v(b), where + v(a) and v(b) vectors holding are suitably-ordered eigenvalues of a and b. - This would introduce some gauge dependence to eigenvalue entanglement - fidelity. However, there's already ambiguity in the definition of this - metric since it requires an ordering of the eigenvalues. It just so - happens that this function's current implementation chooses the ordering - in a gauge-invariant way. """ - _, pairs = _tools.minweight_match(evA, evB, lambda x, y: abs(x - y), - return_pairs=True) # just to get pairing - mlPl = abs(_np.sum([_np.conjugate(evB[j]) * evA[i] for i, j in pairs])) - return 1.0 - mlPl / float(d2) + d2 = a.shape[0] + + if is_unitary is None: + is_unitary = _np.allclose(_np.eye(d2), b @ b.T.conj(), atol=tol, rtol=tol) + if is_tp is None: + is_tp = _tools.is_trace_preserving(a, mx_basis, tol) and _tools.is_trace_preserving(b, mx_basis, tol) + + if is_unitary and is_tp: + evA = _np.linalg.eigvals(a) + evB = _np.linalg.eigvals(b) + _, pairs = _tools.minweight_match(evA, evB, lambda x, y: abs(x - y), + return_pairs=True) # just to get pairing + fid = abs(_np.sum([_np.conjugate(evB[j]) * evA[i] for i, j in pairs])) / d2 + else: + Ja = _tools.fast_jamiolkowski_iso_std(a, mx_basis) + Jb = _tools.fast_jamiolkowski_iso_std(b, mx_basis) + from pygsti.tools.optools import eigenvalue_fidelity + fid = eigenvalue_fidelity(Ja, Jb, gauge_invariant=True) + if fid < 0.9: + pass + # valsA, vecsA = _spl.eigh(Ja) + # valsB, vecsB = _spl.eigh(Jb) + # dissimilarity = lambda vec_x, vec_y : abs(1 - vec_x @ vec_y) + # _, pairs = _tools.minweight_match(vecsA.T.conj(), vecsB.T.conj(), dissimilarity, return_pairs=True) + # fid = abs(_np.sum([_np.conjugate(evB[j]) * evA[i] for i, j in pairs])) / d2 + # fid = _np.linalg.norm(_np.sqrt(evA) * _np.sqrt(evB), ord=1)**2 + + return 1.0 - fid Eigenvalue_entanglement_infidelity = _modf.opsfn_factory(eigenvalue_entanglement_infidelity) @@ -1658,7 +1712,8 @@ def eigenvalue_avg_gate_infidelity(a, b, mx_basis): def eigenvalue_diamondnorm(a, b, mx_basis): """ - Eigenvalue diamond distance between a and b + Eigenvalue diamond distance between `a` and `b`, assuming + that `a` and `b` are TP and that `b` is unitary. Parameters ---------- @@ -1678,8 +1733,8 @@ def eigenvalue_diamondnorm(a, b, mx_basis): d2 = a.shape[0] evA = _np.linalg.eigvals(a) evB = _np.linalg.eigvals(b) - return (d2 - 1.0) / d2 * _np.max(_tools.minweight_match(evA, evB, lambda x, y: abs(x - y), - return_pairs=False)) + temp = _tools.minweight_match(evA, evB, lambda x, y: abs(x - y), return_pairs=False) + return (d2 - 1.0) / d2 * _np.max(temp) Eigenvalue_diamondnorm = _modf.opsfn_factory(eigenvalue_diamondnorm) diff --git a/pygsti/report/section/gauge.py b/pygsti/report/section/gauge.py index 219bba41f..8a92ecf8f 100644 --- a/pygsti/report/section/gauge.py +++ b/pygsti/report/section/gauge.py @@ -167,7 +167,17 @@ 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: + + if switchboard is not None and 'mdl_target_grid' in switchboard: + basis_representative = switchboard['mdl_target_grid'][0].basis + else: + res_representative = kwargs['results'][list(kwargs['results'])[0]] + est_representative = res_representative.estimates[list(res_representative.estimates)[0]] + mdl_representative = est_representative.models[list(est_representative.models)[-1]] + basis_representative = mdl_representative.basis + n_leak_default = None if basis_representative.implies_leakage_modeling else 0 + + if kwargs.get('n_leak', n_leak_default) == 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' ) diff --git a/pygsti/report/section/summary.py b/pygsti/report/section/summary.py index 4d6764763..b841db073 100644 --- a/pygsti/report/section/summary.py +++ b/pygsti/report/section/summary.py @@ -38,7 +38,17 @@ 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): - if kwargs.get('n_leak', 0) == 0: + + if switchboard is not None and 'mdl_target_grid' in switchboard: + basis_representative = switchboard['mdl_target_grid'][0].basis + else: + res_representative = kwargs['results'][list(kwargs['results'])[0]] + est_representative = res_representative.estimates[list(res_representative.estimates)[0]] + mdl_representative = est_representative.models[list(est_representative.models)[-1]] + basis_representative = mdl_representative.basis + n_leak_default = None if basis_representative.implies_leakage_modeling else 0 + + if kwargs.get('n_leak', n_leak_default) == 0: summary_display = ('inf', 'trace', 'diamond', 'evinf', 'evdiamond') else: summary_display = ('sub-inf', 'sub-trace', 'sub-diamond', 'plf-sub-diamond', 'leak-rate-max') diff --git a/pygsti/tools/leakage.py b/pygsti/tools/leakage.py index ead7f0dce..5c4b68908 100644 --- a/pygsti/tools/leakage.py +++ b/pygsti/tools/leakage.py @@ -29,6 +29,9 @@ from pygsti.processors import QubitProcessorSpec +BasisLike = Union[Basis, str] + + NOTATION = \ """ Default notation (possibly superceded by text above) @@ -78,15 +81,28 @@ def assign(fn): *A remark.* The Choi matrix for a superoperator G in S[H] is G(rho_t), where rho_t is the output of this function when n_leak=0. """ + NOTATION) -def tensorized_teststate_density(dim: int, n_leak: int) -> 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() - # +def tensorized_teststate_density(*args) -> np.ndarray: + if len(args) == 2: + dim : int = args[0] # type: ignore + n_leak : int = args[1] # type: ignore + args_ok = isinstance(dim, int) and isinstance(n_leak, int) and n_leak < dim + if not args_ok: + raise ValueError() + E = np.eye(dim, dtype=np.complex128) + if n_leak > 0: + E[-n_leak:,-n_leak:] = 0.0 + elif len(args) == 1: + basis : Basis = args[0] # type: ignore + args_ok = isinstance(basis, Basis) and basis.is_hermitian() and ('I' in basis.labels) + if not args_ok: + raise ValueError() + E = basis.ellookup['I'].copy() + else: + raise ValueError() + psi = pgbt.stdmx_to_stdvec(E).ravel() + psi /= la.norm(psi) + # In a standard leakage basis we have # |psi> = (|00> + |11> + ... + |dim - n_leak - 1>) / sqrt(dim - n_leak). - # rho_test = np.outer(psi, psi) return rho_test @@ -106,18 +122,20 @@ def tensorized_teststate_density(dim: int, n_leak: int) -> np.ndarray: *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) +def apply_tensorized_to_teststate(op_x: np.ndarray, op_y: np.ndarray, op_basis: BasisLike, n_leak: Optional[int]=None) -> tuple[TensorProdBasis, np.ndarray, np.ndarray]: + udim = int(np.sqrt(op_x.shape[0])) + dim = udim**2 + assert op_x.shape == (dim, dim) + assert op_y.shape == (dim, dim) # 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 }. + # U = {rho in M[H] : = 0 for all i >= udim - 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_basis = Basis.cast(op_basis, dim=dim) + std_basis = BuiltinBasis('std', dim) op_x_std = pgbt.change_basis(op_x, op_basis, std_basis) op_y_std = pgbt.change_basis(op_y, op_basis, std_basis) @@ -129,12 +147,15 @@ def apply_tensorized_to_teststate(op_x: np.ndarray, op_y, op_basis: np.ndarray, # # 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) + idle_gate = np.eye(dim, 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) + if n_leak is not None: + rho_test = tensorized_teststate_density(udim, n_leak) + else: + rho_test = tensorized_teststate_density(op_basis) # type: ignore # 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]. @@ -158,6 +179,9 @@ def apply_tensorized_to_teststate(op_x: np.ndarray, op_y, op_basis: np.ndarray, If you only care about P, then you can call superop_subspace_projector instead. """ + NOTATION) def leading_dxd_submatrix_basis_vectors(d: int, n: int, current_basis: Basis) -> np.ndarray: + """ + TODO: deprecate this function. + """ assert d <= n if d == n: return np.eye(n**2) @@ -177,38 +201,86 @@ def leading_dxd_submatrix_basis_vectors(d: int, n: int, current_basis: Basis) -> return submatrix_basis_vectors +# TODO: document me +def projective_measurement_subspace_basis_vectors(E: np.ndarray, basis: Basis) -> np.ndarray: + """ + Our desired subspace is + { rho : ker(rho) contains ker(E) } . + Since E is (proportional to) a projector, this is equivalent to + { rho : E rho E = rho }. + """ + if not basis.is_hermitian(): + raise ValueError() + E = E.copy() + k = np.linalg.matrix_rank(E) + E *= (k/np.trace(E)) + if not pgmt.is_projector(E): + raise ValueError() + proj_elements = [ E @ B @ E for B in basis.elements ] + subspace_frame = np.column_stack([ pgbt.stdmx_to_vec(pB, basis) for pB in proj_elements ]) + # ^ a "frame" is an overcomplete basis. + subspace_frame = subspace_frame.real + # ^ Since basis.elements are Hermitian and E is Hermitian, we have that proj_elements are + # Hermitian and that their vectorizations are real. We cast to real here just to + # eliminate possible rounding errors. + U_full : np.ndarray = la.qr(subspace_frame, pivoting=True)[0] # type: ignore + U = U_full[:, :k**2] + return U + + @lru_cache @set_docstring( -""" -This function returns the superoperator in S[H] that projects orthogonally from -M[H] to M[C], where H is n-dimensional and C ⊂ H is d-dimensional (d ≤ n). +NOTATION) +def superop_subspace_projector(*args) -> np.ndarray: -The action of this operator is easy to understand when M[H] and M[C] are viewed -as spaces of n-by-n matrices rather than spaces of length-n^2 vectors. + num_positional_args = len(args) + if num_positional_args not in {3, 1}: + raise ValueError() -For v in M[H] and u = P v, we have + if num_positional_args == 1: + basis = args[0] + if not isinstance(basis, Basis): + raise ValueError() + if not hasattr(basis, 'ellookup') or ('I' not in basis.ellookup): + raise ValueError() + dim = basis.dim + if basis.first_element_is_identity: + return np.eye(dim) + E = basis.ellookup['I'] + U = projective_measurement_subspace_basis_vectors(E, basis) - mat(v) = [x11, x12] and mat(u) = [x11, 0] - [x21, x22] [ 0, 0]. + else: + """ + This function returns the superoperator in S[H] that projects orthogonally from + M[H] to M[C], where H is n-dimensional and C ⊂ H is d-dimensional (d ≤ n). -This characterization makes two facts about P apparent. First, P is positive -(i.e., it takes Hermitian psd operators to Hermitian psd operators). Second, -P is trace-non-increasing. -""" + NOTATION) -def superop_subspace_projector(d: int, n: int, basis: Basis, force_real=True) -> np.ndarray: - assert d <= n - if d == n: - return np.eye(n**2) - U = leading_dxd_submatrix_basis_vectors(d, n, basis) # type: ignore - P = U @ U.T.conj() - if force_real: - if np.linalg.norm(P.imag) > 1e-12: - msg = "The orthogonal projector onto the computational subspace in the basis\n" - msg += f"{basis} is not real-valued. Since we were passed force_real=True we're\n" - msg += "raising a ValueError. Try again with a basis like 'l2p1' or 'gm', or\n" - msg += "use force_real=False." - raise ValueError(msg) - P = P.real + The action of this operator is easy to understand when M[H] and M[C] are viewed + as spaces of n-by-n matrices rather than spaces of length-n^2 vectors. + + For v in M[H] and u = P v, we have + + mat(v) = [x11, x12] and mat(u) = [x11, 0] + [x21, x22] [ 0, 0]. + + This characterization makes two facts about P apparent. First, P is positive + (i.e., it takes Hermitian psd operators to Hermitian psd operators). Second, + P is trace-non-increasing. + """ + d, n, basis = args + types_ok = isinstance(d, int) and isinstance(n, int) and isinstance(basis, Basis) + if not types_ok: + raise ValueError() + if d > n: + raise ValueError() + if n**2 != basis.elsize: + raise ValueError() + if d == n: + return np.eye(n**2) + E = np.eye(n) + E[d:, d:] = 0.0 + U = projective_measurement_subspace_basis_vectors(E, basis) + + P = U @ U.T return P @@ -227,7 +299,7 @@ def superop_subspace_projector(d: int, n: int, basis: Basis, force_real=True) -> @set_docstring(CHOI_INDUCED_METRIC_TEMPLATE % 'entanglement fidelity') -def subspace_entanglement_fidelity(op_x: np.ndarray, op_y: np.ndarray, op_basis, n_leak=0) -> float: +def subspace_entanglement_fidelity(op_x: np.ndarray, op_y: np.ndarray, op_basis, n_leak: Optional[int] = None) -> float: 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) @@ -236,7 +308,7 @@ def subspace_entanglement_fidelity(op_x: np.ndarray, op_y: np.ndarray, op_basis, @set_docstring(CHOI_INDUCED_METRIC_TEMPLATE % 'jamiolkowski trace distance') -def subspace_jtracedist(op_x: np.ndarray, op_y: np.ndarray, op_basis, n_leak=0) -> float: +def subspace_jtracedist(op_x: np.ndarray, op_y: np.ndarray, op_basis, n_leak: Optional[int] = None) -> float: 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) @@ -258,18 +330,21 @@ def subspace_jtracedist(op_x: np.ndarray, op_y: np.ndarray, op_basis, n_leak=0) @set_docstring(PROJECTION_INDUCED_METRIC_TEMPLATE % 'Frobenius distance') -def subspace_superop_fro_dist(op_x: np.ndarray, op_y: np.ndarray, op_basis, n_leak=0) -> float: +def subspace_superop_fro_dist(op_x: np.ndarray, op_y: np.ndarray, op_basis, n_leak:Optional[int]=None) -> float: diff = op_x - op_y if n_leak == 0: - return np.linalg.norm(diff, 'fro') # type: ignore + return la.norm(diff, 'fro') # type: ignore n = int(np.sqrt(op_x.shape[0])) assert op_x.shape == op_y.shape == (n**2, n**2) - P = superop_subspace_projector(n - n_leak, n, op_basis, force_leak=False) - return np.linalg.norm(diff @ P) # type: ignore + if n_leak is None: + P = superop_subspace_projector(op_basis) + else: + P = superop_subspace_projector(n - n_leak, n, op_basis) + return la.norm(diff @ P) # type: ignore @set_docstring(PROJECTION_INDUCED_METRIC_TEMPLATE % 'diamond distance') -def subspace_diamonddist(op_x: np.ndarray, op_y: np.ndarray, op_basis, n_leak=0) -> float: +def subspace_diamonddist(op_x: np.ndarray, op_y: np.ndarray, op_basis, n_leak: Optional[int]=None) -> float: """ Here we give a brief motivating derivation for defining the subspace diamond norm in the way that we have. This derivation won't convince a skeptic that our definition @@ -316,9 +391,10 @@ def subspace_diamonddist(op_x: np.ndarray, op_y: np.ndarray, op_basis, n_leak=0) from pygsti.tools.optools import diamonddist dim_mixed = op_x.shape[0] dim_pure = int(dim_mixed**0.5) - assert n_leak <= 1 - dim_pure_compsub = dim_pure - n_leak - P = superop_subspace_projector(dim_pure_compsub, dim_pure, op_basis) + if n_leak is None: + P = superop_subspace_projector(op_basis) + else: + P = superop_subspace_projector(dim_pure - n_leak, dim_pure, op_basis) val : float = diamonddist(op_x @ P, op_y @ P, op_basis, return_x=False) / 2 # type: ignore return val @@ -335,9 +411,9 @@ def subspace_diamonddist(op_x: np.ndarray, op_y: np.ndarray, op_basis, n_leak=0) The process matrix `op` is the mx_basis representation of a CPTP map G : M[H] ➞ M[H]. -Consider an experimental protocol where a system is prepared in state ρ, evolved to G(ρ), -and then measured with the 2-element POVM {E_SUB, 1 - E_SUB}. Runs of this protocol that -result in the `1 - E_SUB` measurement outcome are called transport events. +Consider an experimental protocol where a system is prepared in state ρ ∈ M[SUB], evolved +to G(ρ), and then measured with the 2-element POVM {E_SUB, 1 - E_SUB}. Runs of this protocol +that result in the `1 - E_SUB` measurement outcome are called transport events. The population transport profile of (G,SUB) is a full specification of the probabilties of transport events, considering all possible choices of ρ ∈ M[SUB]. This profile is @@ -397,8 +473,8 @@ def pop_transport_profile(E_sub: np.ndarray, op: np.ndarray, mx_basis: Basis, E_ transport_E_mat = pgbt.vec_to_stdmx(transport_E_vec, mx_basis, keep_complex=True) transport_E_mat = E_sub @ transport_E_mat @ E_sub - rates, states = np.linalg.eigh(transport_E_mat) - dim_proj = int(np.trace(E_sub).real) + rates, states = la.eigh(transport_E_mat) + dim_proj = int(np.round(np.trace(E_sub).real)) ind = np.argsort(np.abs(rates))[::-1][:dim_proj] rates = rates[ind] states = [s for s in states.T[ind]] @@ -594,10 +670,18 @@ def lagoified_gopparams_dicts(gopparams_dicts: List[Dict]) -> List[Dict]: """ from pygsti.models.gaugegroup import UnitaryGaugeGroup tm = gopparams_dicts[0]['target_model'] + if tm.basis.implies_leakage_modeling: + # The basis carries enough information for us to infer how leakage should + # be modeled. + n_leak_arg = None + else: + # The basis doesn't say anything about how the leakage should be modeled. + # We'll assume a 3-level system with one leakage level. + n_leak_arg = 1 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 + inner_dict['n_leak'] = n_leak_arg # ^ This function could accept n_leak as an argument instead. However, # downstream functions for gauge optimization only support n_leak=0 or 1. # @@ -605,8 +689,8 @@ 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['gates_metric'] = 'frobenius squared' - inner_dict['spam_metric'] = 'frobenius squared' + inner_dict['gates_metric'] = 'fidelity' + inner_dict['spam_metric'] = 'fidelity' inner_dict['item_weights'] = {'gates': 0.0, 'spam': 1.0} gg = UnitaryGaugeGroup(tm.basis.state_space, tm.basis) inner_dict['gauge_group'] = gg @@ -631,7 +715,9 @@ def lagoified_gopparams_dicts(gopparams_dicts: List[Dict]) -> List[Dict]: if gg is not None: inner_dict['gauge_group'] = gg inner_dict['_gaugeGroupEl'] = gg.compute_element(gg.initial_params) - inner_dict['n_leak'] = 1 + inner_dict['gates_metric'] = 'frobenius squared' + inner_dict['spam_metric'] = 'frobenius squared' + inner_dict['n_leak'] = n_leak_arg inner_dict['item_weights'] = {'gates': 1.0, 'spam': 1.0} gopparams_dicts.append(inner_dict) return gopparams_dicts @@ -720,8 +806,11 @@ def construct_leakage_report( for ek in est_key: assert isinstance(ek, str) add_lago_models(results, ek, verbosity=gaugeopt_verbosity) + est = results.estimates[ek] + n_leak_arg = None if est.models['LAGO'].basis.implies_leakage_modeling else 1 + from pygsti.report import construct_standard_report report = construct_standard_report( - results, advanced_options={'n_leak': 1}, **extra_report_kwargs + results, advanced_options={'n_leak': n_leak_arg}, **extra_report_kwargs ) return report, results diff --git a/pygsti/tools/matrixtools.py b/pygsti/tools/matrixtools.py index 8be082733..04728cba2 100644 --- a/pygsti/tools/matrixtools.py +++ b/pygsti/tools/matrixtools.py @@ -65,7 +65,7 @@ def gram_matrix(m, adjoint=False): return out -def is_hermitian(mx, tol=1e-9): +def is_hermitian(mx: _np.ndarray, tol: float = 1e-9) -> bool: """ Test whether mx is a hermitian matrix. @@ -86,7 +86,17 @@ def is_hermitian(mx, tol=1e-9): if m != n: return False else: - return _np.all(_np.abs(mx - mx.T.conj()) <= tol) + return _np.all(_np.abs(mx - mx.T.conj()) <= tol) # type: ignore + + +def assert_hermitian(mat : _np.ndarray, tol: _np.floating) -> None: + hermiticity_error = _np.abs(mat - mat.T.conj()) + if _np.any(hermiticity_error > tol): + message = f""" + Input matrix 'mat' is not Hermitian, up to tolerance {tol}. + The absolute values of entries in (mat - mat^H) are \n{hermiticity_error}. + """ + raise ValueError(message) def is_pos_def(mx, tol=1e-9, attempt_cholesky=False): diff --git a/pygsti/tools/optools.py b/pygsti/tools/optools.py index 5dab09e8a..24c23ebac 100644 --- a/pygsti/tools/optools.py +++ b/pygsti/tools/optools.py @@ -86,18 +86,9 @@ def fidelity(a, b): # ^ use for checks that have no dimensional dependence; about 1e-12 for double precision. __VECTOR_TOL__ = (a.shape[0] ** 0.5) * __SCALAR_TOL__ # ^ use for checks that do have dimensional dependence (will naturally increase for larger matrices) - - def assert_hermitian(mat): - hermiticity_error = _np.abs(mat - mat.T.conj()) - if _np.any(hermiticity_error > __SCALAR_TOL__): - message = f""" - Input matrix 'mat' is not Hermitian, up to tolerance {__SCALAR_TOL__}. - The absolute values of entries in (mat - mat^H) are \n{hermiticity_error}. - """ - raise ValueError(message) - assert_hermitian(a) - assert_hermitian(b) + _mt.assert_hermitian(a, __SCALAR_TOL__) + _mt.assert_hermitian(b, __VECTOR_TOL__) def check_rank_one_density(mat): """ @@ -184,6 +175,32 @@ def psd_square_root(mat): return f +def eigenvalue_fidelity(a, b, gauge_invariant=True) -> _np.floating: + from pygsti.tools import minweight_match + tol = _np.finfo(a.dtype).eps ** 0.75 + _mt.assert_hermitian(a, tol) + _mt.assert_hermitian(b, tol) + if gauge_invariant: + valsA = _spl.eigvalsh(a) + valsB = _spl.eigvalsh(b) + dissimilarity = lambda x, y: abs(x - y) + _, pairs = minweight_match(valsA, valsB, dissimilarity, return_pairs=True) + else: + valsA, vecsA = _spl.eigh(a) + valsB, vecsB = _spl.eigh(b) + dissimilarity = lambda vec_x, vec_y : abs(1 - vec_x @ vec_y) + _, pairs = minweight_match(vecsA.T.conj(), vecsB.T.conj(), dissimilarity, return_pairs=True) + ind_a, ind_b = zip(*pairs) + arg_a = _np.maximum(valsA[list(ind_a)], 0) + arg_b = _np.maximum(valsB[list(ind_b)], 0) + f = _np.linalg.norm(arg_a**0.5 * arg_b**0.5, ord=1) + return f + + +def eigenvalue_infidelity(a, b, gauge_invariant=True) -> _np.floating: + return 1 - eigenvalue_fidelity(a, b, gauge_invariant) + + def frobeniusdist(a, b) -> _np.floating[Any]: """ Returns the frobenius distance between arrays: ||a - b||_Fro. @@ -377,22 +394,40 @@ def jtracedist(a, b, mx_basis='pp'): # Jamiolkowski trace distance: Tr(|J(a)-J return tracedist(JA, JB) +def is_trace_preserving(a: _np.ndarray, mx_basis: BasisLike='pp', tol=1e-14) -> bool: + dim = int(_np.round( a.size**0.5 )) + udim = int(_np.round( dim**0.5 )) + basis = _Basis.cast(mx_basis, dim=dim) + if basis.first_element_is_identity: + checkone = _np.isclose(a[0,0], 1.0, rtol=tol, atol=tol) + checktwo = _np.allclose(a[1:], 0.0, atol=tol) + return checkone and checktwo + # else, check that the adjoint of a is unital. + I_mat = _np.eye(udim) + I_vec = _bt.stdmx_to_vec(I_mat, basis) + if _np.isrealobj(a): + expect_I_vec = a.T @ I_vec + else: + expect_I_vec = a.T.conj() @ I_vec + atol = tol * udim + check = float(_np.linalg.norm(I_vec - expect_I_vec)) <= atol + return check + + def entanglement_fidelity(a, b, mx_basis: BasisLike='pp', is_tp=None, is_unitary=None): """ Returns the "entanglement" process fidelity between gate matrices. - This is given by: - - `F = Tr( sqrt{ sqrt(J(a)) * J(b) * sqrt(J(a)) } )^2` - - where J(.) is the Jamiolkowski isomorphism map that maps a operation matrix - to it's corresponding Choi Matrix. + The entanglement fidelity of (a, b) is defined as the fidelity of (J(a), J(b)), + where J(.) is the Jamiolkowski isomorphism map. - When the both of the input matrices a and b are TP, and - the target matrix b is unitary then we can use a more efficient - formula: + When the both of the input matrices a and b are TP, and the target matrix b + is unitary, then we can use the formula - `F= Tr(a @ b.conjugate().T)/d^2` + `F = Tr(a @ b.conjugate().T)/d^2` + + That formula avoids forming J(a) or J(b), and it avoids the eigenvalue computations + needed to evaluate fidelity in the general case. Parameters ---------- @@ -439,26 +474,14 @@ def entanglement_fidelity(a, b, mx_basis: BasisLike='pp', is_tp=None, is_unitary #if the tp flag isn't set we'll calculate whether it is true here if is_tp is None: - def is_tp_fn(x): - return _np.isclose(x[0, 0], 1.0) and _np.allclose(x[0, 1:d2], 0) - - is_tp= (is_tp_fn(a) and is_tp_fn(b)) + is_tp = is_trace_preserving(a, mx_basis) and is_trace_preserving(b, mx_basis) #if the unitary flag isn't set we'll calculate whether it is true here if is_unitary is None: - is_unitary= _np.allclose(_np.identity(d2, 'd'), _np.dot(b, b.conjugate().T)) + is_unitary = _np.allclose(_np.identity(d2, 'd'), b @ b.T.conj()) if is_tp and is_unitary: # then assume TP-like gates & use simpler formula - #old version, slower than einsum - #TrLambda = _np.trace(_np.dot(a, b.conjugate().T)) # same as using _np.linalg.inv(b) - - #Use einsum black magic to only calculate the diagonal elements - #if the basis is either pp or gm we know the elements are real-valued, so we - #don't need to take the conjugate - if mx_basis=='pp' or mx_basis=='gm': - TrLambda = _np.einsum('ij,ji->',a, b.T) - else: - TrLambda = _np.einsum('ij,ji->',a, b.conjugate().T) + TrLambda = _np.vdot(b, a) # == Trace(b.T.conj() @ a) == Trace(a @ b.T.conj()) return TrLambda / d2 JA = _jam.jamiolkowski_iso(a, mx_basis, mx_basis) @@ -637,6 +660,7 @@ def entanglement_infidelity(a, b, mx_basis: BasisLike = 'pp', is_tp=None, is_uni """ return 1 - entanglement_fidelity(a, b, mx_basis, is_tp, is_unitary) + def generator_infidelity(a, b, mx_basis = 'pp'): """ Returns the generator infidelity between a and b, where b is the "target" operation. @@ -687,6 +711,7 @@ def generator_infidelity(a, b, mx_basis = 'pp'): return _np.real_if_close(gen_infid) + def gateset_infidelity(model, target_model, itype='EI', weights=None, mx_basis=None, is_tp=None, is_unitary=None): """ diff --git a/pygsti/tools/sdptools.py b/pygsti/tools/sdptools.py index 071195a20..11da0dcc6 100644 --- a/pygsti/tools/sdptools.py +++ b/pygsti/tools/sdptools.py @@ -14,7 +14,7 @@ import numpy as np -from typing import Union, List, Tuple, TYPE_CHECKING +from typing import Union, List, Tuple, Optional, TYPE_CHECKING if TYPE_CHECKING: import cvxpy as cp ExpressionLike = Union[cp.Expression, np.ndarray] @@ -159,23 +159,30 @@ def cptp_superop_variable(purestate_dim: int, basis: BasisLike) -> Tuple[cp.Expr 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): +def diamond_distance_projection_model(superop: np.ndarray, basis: BasisLike, leakfree: bool=False, seepfree: bool=False, n_leak: Optional[int]=ModuleNotFoundError, 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 = [] + from pygsti.tools.leakage import leading_dxd_submatrix_basis_vectors, projective_measurement_subspace_basis_vectors + + if not isinstance(basis, Basis): + basis = Basis.cast(basis, dim=dim_mixed) + if n_leak is None: + U = projective_measurement_subspace_basis_vectors(basis.ellookup['I'], basis) + else: + dim_pure_compsub = dim_pure - n_leak + U = leading_dxd_submatrix_basis_vectors(dim_pure_compsub, dim_pure, basis) # type: ignore + 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: