Update germ selection and FPR to better work with implicit models

pygsti/algorithms/fiducialpairreduction.py

@@ -24,7 +24,10 @@
 from pygsti import circuits as _circuits
 
 from pygsti.circuits import circuitconstruction as _gsc
+from pygsti.models import ExplicitOpModel as _ExplicitOpModel
+from pygsti.models import ImplicitOpModel as _ImplicitOpModel
 from pygsti.modelmembers.operations import EigenvalueParamDenseOp as _EigenvalueParamDenseOp
+from pygsti.modelmembers.povms import convert as _convert_povm
 from pygsti.tools import remove_duplicates as _remove_duplicates
 from pygsti.tools import slicetools as _slct
 from pygsti.tools.legacytools import deprecate as _deprecated_fn
@@ -428,8 +431,7 @@ def find_sufficient_fiducial_pairs_per_germ(target_model, prep_fiducials, meas_f
             #Create a new model containing static target gates and a
             # special "germ" gate that is parameterized only by it's
             # eigenvalues (and relevant off-diagonal elements)
-            gsGerm = target_model.copy()
-            gsGerm.set_all_parameterizations("static")
+            gsGerm = _copy_to_static_explicitop_model(target_model)
             germMx = gsGerm.sim.product(germ)
             #give this state space labels equal to the line_labels of 
             gsGerm.operations['Ggerm'] = _EigenvalueParamDenseOp(
@@ -628,8 +630,7 @@ def find_sufficient_fiducial_pairs_per_germ_greedy(target_model, prep_fiducials,
             #Create a new model containing static target gates and a
             # special "germ" gate that is parameterized only by it's
             # eigenvalues (and relevant off-diagonal elements)
-            gsGerm = target_model.copy()
-            gsGerm.set_all_parameterizations("static")
+            gsGerm = _copy_to_static_explicitop_model(target_model)
             germMx = gsGerm.sim.product(germ)
             gsGerm.operations["Ggerm"] = _EigenvalueParamDenseOp(
                 germMx, True, constrain_to_tp)
@@ -804,8 +805,7 @@ def find_sufficient_fiducial_pairs_per_germ_power(target_model, prep_fiducials,
             #Create a new model containing static target gates and a
             # special "germ" gate that is parameterized only by it's
             # eigenvalues (and relevant off-diagonal elements)
-            gsGerm = target_model.copy()
-            gsGerm.set_all_parameterizations("static")
+            gsGerm = _copy_to_static_explicitop_model(target_model)
             germMx = gsGerm.sim.product(germ)
 
             gsGerm.operations["Ggerm"] = _EigenvalueParamDenseOp(
@@ -2071,6 +2071,29 @@ def _make_spam_static(model):
     model_copy._rebuild_paramvec()
 
     return model_copy
+
+
+def _copy_to_static_explicitop_model(mdl):
+    if isinstance(mdl, _ExplicitOpModel):
+        ret = mdl.copy()
+        ret.set_all_parameterizations('static')
+    else:
+        ret = _ExplicitOpModel(mdl.state_space, mdl.basis,
+                           default_gate_type='static',
+                           default_prep_type='static',
+                           default_povm_type='static',
+                           default_instrument_type='static',
+                           simulator=mdl.sim.copy(), evotype=mdl.evotype)
+        for k, v in mdl.prep_blks['layers'].items():
+            ret.preps[k] = v.to_dense()
+        for k, v in mdl.povm_blks['layers'].items():
+            ret.povms[k] = _convert_povm(v.copy(), 'static', mdl.basis)
+        for k, v in mdl.operation_blks['layers'].items():
+            ret.operations[k] = v.to_dense()
+
+    assert ret.num_params == 0
+    return ret
+
 
 #write a helper function for precomputing the jacobian dictionaries from bulk_dprobs
 #which can then be passed into the construction of the compactEVD caches.
@@ -2147,10 +2170,18 @@ def compute_jacobian_dicts(model, germs, prep_fiducials, meas_fiducials, prep_po
 def _set_up_prep_POVM_tuples(target_model, prep_povm_tuples, return_meas_dofs= False):
 
     if prep_povm_tuples == "first":
-        firstRho = list(target_model.preps.keys())[0]
-        prep_ssl = [target_model.preps[firstRho].state_space.state_space_labels]
-        firstPOVM = list(target_model.povms.keys())[0]
-        POVM_ssl = [target_model.povms[firstPOVM].state_space.state_space_labels]
+        if isinstance(target_model, _ExplicitOpModel):        
+            firstRho = list(target_model.preps.keys())[0]
+            prep_ssl = [target_model.preps[firstRho].state_space.state_space_labels]
+            firstPOVM = list(target_model.povms.keys())[0]
+            POVM_ssl = [target_model.povms[firstPOVM].state_space.state_space_labels]
+        elif isinstance(target_model, _ImplicitOpModel):
+            firstRho = list(target_model.prep_blks['layers'].keys())[0]
+            prep_ssl = [target_model.prep_blks['layers'][firstRho].state_space.state_space_labels]
+            firstPOVM = list(target_model.povm_blks['layers'].keys())[0]
+            POVM_ssl = [target_model.povm_blks['layers'][firstPOVM].state_space.state_space_labels]
+        else:
+            raise ValueError("Model must be an ExplicitOpModel or ImplicitOpModel")
         prep_povm_tuples = [(firstRho, firstPOVM)]
         #I think using the state space labels for firstRho and firstPOVM as the
         #circuit labels should work most of the time (new stricter line_label enforcement means
@@ -2160,11 +2191,24 @@ def _set_up_prep_POVM_tuples(target_model, prep_povm_tuples, return_meas_dofs= F
     #if not we still need to extract state space labels for all of these to meet new circuit
     #label handling requirements.
     else:
-        prep_ssl = [target_model.preps[lbl_tup[0]].state_space.state_space_labels for lbl_tup in prep_povm_tuples]
-        POVM_ssl = [target_model.povms[lbl_tup[1]].state_space.state_space_labels for lbl_tup in prep_povm_tuples]
+        if isinstance(target_model, _ExplicitOpModel):
+            prep_ssl = [target_model.preps[lbl_tup[0]].state_space.state_space_labels for lbl_tup in prep_povm_tuples]
+            POVM_ssl = [target_model.povms[lbl_tup[1]].state_space.state_space_labels for lbl_tup in prep_povm_tuples]
+        elif isinstance(target_model, _ImplicitOpModel):
+            prep_ssl = [target_model.prep_blks['layers'][lbl_tup[0]].state_space.state_space_labels for lbl_tup in prep_povm_tuples]
+            POVM_ssl = [target_model.povm_blks['layers'][lbl_tup[1]].state_space.state_space_labels for lbl_tup in prep_povm_tuples]
+        else:
+            raise ValueError("Model must be an ExplicitOpModel or ImplicitOpModel")
 
     #brief intercession to calculate the number of degrees of freedom for the povm.
-    num_effects= len(list(target_model.povms[prep_povm_tuples[0][1]].keys()))
+    if isinstance(target_model, _ExplicitOpModel):
+        num_effects= len(list(target_model.povms[prep_povm_tuples[0][1]].keys()))
+    elif isinstance(target_model, _ImplicitOpModel):
+        num_effects= len(list(target_model.povm_blks['layers'][prep_povm_tuples[0][1]].keys()))
+    else:
+        raise ValueError("Model must be an ExplicitOpModel or ImplicitOpModel")
+
+    #subtract 1 for the POVM constraint that effects sum to identity.
     dof_per_povm= num_effects-1
 
     prep_povm_tuples = [(_circuits.Circuit([prepLbl], line_labels=prep_ssl[i]), 

pygsti/algorithms/germselection.py

@@ -25,6 +25,7 @@
 from pygsti import baseobjs as _baseobjs
 from pygsti.tools import mpitools as _mpit
 from pygsti.models import ExplicitOpModel as _ExplicitOpModel
+from pygsti.models import ImplicitOpModel as _ImplicitOpModel
 from pygsti.forwardsims import MatrixForwardSimulator as _MatrixForwardSimulator
 
 FLOATSIZE = 8  # in bytes: TODO: a better way
@@ -193,7 +194,13 @@ def find_germs(target_model, randomize=True, randomization_strength=1e-2,
     modelList = _setup_model_list(target_model, randomize,
                                   randomization_strength, num_gs_copies, seed)
 
-    gates = list(target_model.operations.keys())
+    if isinstance(target_model, _ExplicitOpModel):
+        gates = list(target_model.operations.keys())
+    elif isinstance(target_model, _ImplicitOpModel):
+        gates = list(target_model.operation_blks['layers'].keys())
+    else:
+        raise ValueError("Model must be an ExplicitOpModel or ImplicitOpModel")
+
     availableGermsList = []
     if candidate_germ_counts is None: candidate_germ_counts = {6: 'all upto'}
     for germLength, count in candidate_germ_counts.items():
@@ -1020,17 +1027,20 @@ def _remove_spam_vectors(model):
     Model
     """
     reducedModel = model.copy()
-    try:
+    if isinstance(reducedModel, _ExplicitOpModel):
         for prepLabel in list(reducedModel.preps.keys()):
             del reducedModel.preps[prepLabel]
         for povmLabel in list(reducedModel.povms.keys()):
             del reducedModel.povms[povmLabel]
-    except AttributeError:
+    elif isinstance(reducedModel, _ImplicitOpModel):
         # Implicit model instead
         for prepLabel in list(reducedModel.prep_blks.keys()):
             del reducedModel.prep_blks[prepLabel]
         for povmLabel in list(reducedModel.povm_blks.keys()):
             del reducedModel.povm_blks[povmLabel]
+    else:
+        raise ValueError("Don't know how to remove SPAM vectors from a "
+                         f"{type(reducedModel)}")
 
     reducedModel._mark_for_rebuild()
     return reducedModel
@@ -1303,10 +1313,17 @@ def _bulk_twirled_deriv(model, circuits, eps=1e-6, check=False, comm=None, float
     numpy array
         An array of shape (num_simplified_circuits, op_dim^2, num_model_params)
     """
-    if len(model.preps) > 0 or len(model.povms) > 0:
-        model = _remove_spam_vectors(model)
-        # This function assumes model has no spam elements so `lookup` below
-        #  gives indexes into products computed by evalTree.
+
+    # This function assumes model has no spam elements so `lookup` below
+    #  gives indexes into products computed by evalTree.
+    if isinstance(model, _ExplicitOpModel):
+        if len(model.preps) > 0 or len(model.povms) > 0:
+            model = _remove_spam_vectors(model)            
+    elif isinstance(model, _ImplicitOpModel):
+        if len(model.prep_blks) > 0 or len(model.povm_blks) > 0:
+            model = _remove_spam_vectors(model)
+    else:
+        raise ValueError(f"Unknown model type: {type(model)}")
 
     resource_alloc = _baseobjs.ResourceAllocation(comm=comm)
     dProds, prods = model.sim.bulk_dproduct(circuits, flat=True, return_prods=True, resource_alloc=resource_alloc)
@@ -1786,11 +1803,15 @@ def find_germs_breadthfirst(model_list, germs_list, randomize=True,
     #assert(all([(mdl.num_params == Np) for mdl in model_list])), \
     #    "All models must have the same number of parameters!"
 
-    (_, numGaugeParams,
-     numNonGaugeParams, _) = _get_model_params(model_list)
     if num_nongauge_params is not None:
-        numGaugeParams = numGaugeParams + numNonGaugeParams - num_nongauge_params
+        reducedModelList = list(map(_remove_spam_vectors, model_list))
+        numParamsSet = {reducedModel.num_params for reducedModel in reducedModelList}
+        if len(numParamsSet) != 1:
+            raise ValueError("When specifying num_nongauge_params, all models must have the same number of parameters.")
+        numGaugeParams = next(iter(numParamsSet)) - num_nongauge_params
         numNonGaugeParams = num_nongauge_params
+    else:
+        (_, numGaugeParams, numNonGaugeParams, _) = _get_model_params(model_list)
 
     germLengths = _np.array([len(germ) for germ in germs_list], _np.int64)
 
@@ -3037,9 +3058,16 @@ def _compute_bulk_twirled_ddd_compact(model, germs_list, eps,
     #The representations of the germ process matrices are clearly independent 
     #of the spam parameters. (I say that, but I only realized I had forgotten this like
     #6 months later...)
-    if len(model.preps) > 0 or len(model.povms) > 0:
-        model = _remove_spam_vectors(model)
-        # This function assumes model has no spam elements so `lookup` below
+    if isinstance(model, _ExplicitOpModel):
+        if len(model.preps) > 0 or len(model.povms) > 0:
+            model = _remove_spam_vectors(model)
+    elif isinstance(model, _ImplicitOpModel):
+        if len(model.prep_blks['layers']) > 0 or \
+           len(model.povm_blks['layers']) > 0:
+            model = _remove_spam_vectors(model)
+    else:
+        raise ValueError(f"Unknown model type: {type(model)}.")
+    # This function assumes model has no spam elements so `lookup` below
 
     if printer is not None:
         printer.log('Generating compact EVD Cache',1)