diff --git a/pygsti/protocols/gst.py b/pygsti/protocols/gst.py index 64c015635..dba21e8e6 100644 --- a/pygsti/protocols/gst.py +++ b/pygsti/protocols/gst.py @@ -2294,43 +2294,41 @@ def _add_badfit_estimates(results, base_estimate_label, badfit_options, def _compute_wildcard_budget_1d_model(estimate, objfn_cache, mdc_objfn, parameters, badfit_options, 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 by Tim and Stefan in the RCSGST paper. - TODO: docstring (update) + using the model introduced in https://doi.org/10.1038/s41534-023-00764-y. Parameters ---------- - model : Model - The model to add a wildcard budget to. + estimate : Estimate + The estimate object containing the model and data to be used. - ds : DataSet - The data the model predictions are being compared with. + objfn_cache : ObjectiveFunctionCache + A cache for objective function evaluations. - circuits_to_use : list - The circuits whose data are compared. + 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 + a fit is unsatisfactory. Contains detailed parameters for wildcard budget creation. - comm : mpi4py.MPI.Comm, optional - An MPI communicator used to run this computation in parallel. - - mem_limit : int, optional - A rough per-processor memory limit in bytes. - verbosity : int, optional Level of detail printed to stdout. - + gaugeopt_suite : GSTGaugeOptSuite, optional (default None) - The 1-D wildcard models can be based on gauge-variant functions such as the diamond distance. When gauge optimization is applied this will ensure the gauge optimized models used are consistent in the wildcard model and the parent estimates. + The 1-D wildcard models can be based on gauge-variant functions such as the diamond distance. + When gauge optimization is applied, this will ensure the gauge optimized models used are consistent + in the wildcard model and the parent estimates. Returns ------- - PrimitiveOpsWildcardBudget + PrimitiveOpsWildcardBudget or dict + The computed wildcard budget. If gauge optimization is applied, a dictionary of + budgets keyed by gauge optimization labels is returned. + """ printer = _baseobjs.VerbosityPrinter.create_printer(verbosity, mdc_objfn.resource_alloc) badfit_options = GSTBadFitOptions.cast(badfit_options) @@ -2387,9 +2385,32 @@ def _compute_wildcard_budget_1d_model(estimate, objfn_cache, mdc_objfn, paramete def _compute_1d_reference_values_and_name(estimate, badfit_options, gaugeopt_suite=None): - ''' - DOCSTRING: TODO - ''' + """ + Compute the reference values and name for 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 + creation. + + gaugeopt_suite : GSTGaugeOptSuite, optional (default None) + The gauge optimization suite used by the estimate. If None, final iteration + estimate will be used. + + 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. + """ 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'] @@ -2398,10 +2419,13 @@ def _compute_1d_reference_values_and_name(estimate, badfit_options, gaugeopt_sui 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 + 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 + else: # Local/cloud noise models don't have default_gauge_group attribute and can't be gauge # optimized - at least not easily. @@ -2419,6 +2443,14 @@ def _compute_1d_reference_values_and_name(estimate, badfit_options, gaugeopt_sui _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 = .5* _tools.instrument_diamonddist(op, targetinsts_dict[key]) + if inst_dd < 0: # indicates that instrument_diamonddist failed + _warnings.warn(("Diamond distance failed to compute %s reference value for 1D wildcard budget!" + "No fallback presently available for instruments, so skipping.") % str(key)) + else: + dd[key] = inst_dd spamdd = {} for key, op in preps_dict.items(): @@ -2444,15 +2476,22 @@ def _compute_1d_reference_values_and_name(estimate, badfit_options, gaugeopt_sui " Falling back to trace distance.") % str(key)) dd[lbl][key] = _tools.jtracedist(op.to_dense(), target_model.operations[key].to_dense()) - spamdd = {} - for key, op in gaugeopt_model.preps.items(): - spamdd[key] = _tools.tracedist(_tools.vec_to_stdmx(op.to_dense(), 'pp'), - _tools.vec_to_stdmx(target_model.preps[key].to_dense(), 'pp')) + for key, op in gaugeopt_model.instruments.items(): + inst_dd = .5* _tools.instrument_diamonddist(op, target_model.instruments[key], gaugeopt_model.basis) + if inst_dd < 0: # indicates that instrument_diamonddist failed + _warnings.warn(("Diamond distance failed to compute %s reference value for 1D wildcard budget!" + "No fallback presently available for instruments, so skipping.") % str(key)) + else: + dd[lbl][key] = inst_dd + spamdd = {} + for key, op in gaugeopt_model.preps.items(): + spamdd[key] = _tools.tracedist(_tools.vec_to_stdmx(op.to_dense(), 'pp'), + _tools.vec_to_stdmx(target_model.preps[key].to_dense(), 'pp')) - for key in gaugeopt_model.povms.keys(): - spamdd[key] = 0.5 * _tools.optools.povm_diamonddist(gaugeopt_model, target_model, key) + for key in gaugeopt_model.povms.keys(): + spamdd[key] = 0.5 * _tools.optools.povm_diamonddist(gaugeopt_model, target_model, key) - dd[lbl]['SPAM'] = sum(spamdd.values()) + dd[lbl]['SPAM'] = sum(spamdd.values()) return dd, 'diamond distance' else: raise ValueError("Invalid wildcard1d_reference value (%s) in bad-fit options!" diff --git a/pygsti/tools/optools.py b/pygsti/tools/optools.py index 8d68a73ff..b8d8995dd 100644 --- a/pygsti/tools/optools.py +++ b/pygsti/tools/optools.py @@ -969,6 +969,64 @@ def povm_diamonddist(model, target_model, povmlbl): target_povm_mx = compute_povm_map(target_model, povmlbl) return diamonddist(povm_mx, target_povm_mx, target_model.basis) +def instrument_infidelity(a, b, mx_basis): + """ + Infidelity between instruments a and b + + Parameters + ---------- + a : Instrument + The first instrument. + + b : Instrument + The second instrument. + + mx_basis : Basis or {'pp', 'gm', 'std'} + the basis that `a` and `b` are in. + + Returns + ------- + float + """ + sqrt_component_fidelities = [_np.sqrt(entanglement_fidelity(a[l], b[l], mx_basis)) + for l in a.keys()] + return 1 - sum(sqrt_component_fidelities)**2 + + +def instrument_diamonddist(a, b, mx_basis): + """ + The diamond distance between instruments a and b. + + Parameters + ---------- + a : Instrument + The first instrument. + + b : Instrument + The second instrument. + + mx_basis : Basis or {'pp', 'gm', 'std'} + the basis that `a` and `b` are in. + + Returns + ------- + float + """ + #Turn instrument into a CPTP map on qubit + classical space. + adim = a.state_space.dim + mx_basis = _Basis.cast(mx_basis, dim=adim) + nComps = len(a.keys()) + sumbasis = _DirectSumBasis([mx_basis] * nComps) + composite_op = _np.zeros((adim * nComps, adim * nComps), 'd') + composite_top = _np.zeros((adim * nComps, adim * nComps), 'd') + for i, clbl in enumerate(a.keys()): + aa, bb = i * adim, (i + 1) * adim + for j in range(nComps): + cc, dd = j * adim, (j + 1) * adim + composite_op[aa:bb, cc:dd] = a[clbl].to_dense(on_space='HilbertSchmidt') + composite_top[aa:bb, cc:dd] = b[clbl].to_dense(on_space='HilbertSchmidt') + return diamonddist(composite_op, composite_top, sumbasis) + #decompose operation matrix into axis of rotation, etc def decompose_gate_matrix(operation_mx):