Merge pull request #527 from sandialabs/feature-instrument-wildcard

Wildcard Support for Instruments

Author: sserita

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!"

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):