From 106ecbe390e8e4c147619e5a2f5b8b921db9e409 Mon Sep 17 00:00:00 2001
From: Juan G Mendoza <jgonz19@sandia.gov>
Date: Tue, 24 Jun 2025 14:28:30 -0600
Subject: [PATCH 01/34] recovered mostly changes to allow for FOGI model
 serialization

---
 pygsti/baseobjs/errorgenbasis.py | 157 ++++++++++++++++++++++++-
 pygsti/baseobjs/errorgenspace.py |  30 ++++-
 pygsti/baseobjs/statespace.py    |   3 +-
 pygsti/models/explicitmodel.py   |  20 +++-
 pygsti/models/fogistore.py       | 194 +++++++++++++++++++++----------
 pygsti/models/model.py           | 193 +++++++++++++++---------------
 6 files changed, 425 insertions(+), 172 deletions(-)

diff --git a/pygsti/baseobjs/errorgenbasis.py b/pygsti/baseobjs/errorgenbasis.py
index dd0a6ea70..99333dfe8 100644
--- a/pygsti/baseobjs/errorgenbasis.py
+++ b/pygsti/baseobjs/errorgenbasis.py
@@ -16,11 +16,12 @@
 from pygsti.baseobjs import Basis as _Basis
 from pygsti.baseobjs.errorgenlabel import GlobalElementaryErrorgenLabel as _GlobalElementaryErrorgenLabel,\
 LocalElementaryErrorgenLabel as _LocalElementaryErrorgenLabel
-
+from pygsti.baseobjs.nicelyserializable import NicelySerializable as _NicelySerializable
 from pygsti.tools import optools as _ot
+from pygsti.baseobjs.statespace import StateSpace as _StateSpace
 
 
-class ElementaryErrorgenBasis(object):
+class ElementaryErrorgenBasis(_NicelySerializable):
     """
     A basis for error-generator space defined by a set of elementary error generators.
 
@@ -88,6 +89,7 @@ def __init__(self, state_space, labels, basis_1q=None):
             comprise the basis element labels for the values of the
             `ElementaryErrorgenLabels` in `labels`.
         """
+        super().__init__()
         labels = tuple(labels)
 
         #add an assertion that the labels are ElementaryErrorgenLabels and that all of the labels are the same type.
@@ -116,6 +118,18 @@ def __init__(self, state_space, labels, basis_1q=None):
         self._cached_dual_matrices = None
         self._cached_supports = None
 
+    def _to_nice_serialization(self):
+        state = super()._to_nice_serialization()
+        state.update({'state_space' : self.state_space._to_nice_serialization(),
+                      'labels' : [label.__str__() for label in self.labels],
+                      '_basis_1q' : self._basis_1q if isinstance(self._basis_1q, str) else self._basis_1q._to_nice_serialization()
+
+        }
+        )
+        return state
+    @classmethod
+    def from_nice_serialization(cls, state):
+        return cls(_StateSpace.from_nice_serialization(state['state_space']), [_GlobalElementaryErrorgenLabel.cast(label) for label in state['labels']], state['_basis_1q'] if isinstance(state['_basis_1q'], str) else _Basis.from_nice_serialization(state['_basis_1q']))
     @property
     def labels(self):
         return self._labels
@@ -259,7 +273,6 @@ def union(self, other_basis):
         #assert that these two bases have compatible label types.
         msg = 'Incompatible `ElementaryErrrogenLabel` types, the two `ElementaryErrorgenBasis` should have the same label type.'
         assert type(self._labels[0]) == type(other_basis.labels[0]), msg
-        
         #Get the union of the two bases labels.
         union_labels = set(self._labels) | set(other_basis.labels)
         union_state_space = self.state_space.union(other_basis.state_space)
@@ -741,6 +754,15 @@ def label_index(self, label, ok_if_missing=False, identity_label='I'):
         identity_label : str, optional (default 'I')
             An optional string specifying the label used to denote the identity in basis element labels.
         """
+        try:
+            return self.labels.index(label)
+        except ValueError as error:
+            
+            if ok_if_missing:
+                return None
+            else:
+                raise error
+        '''
         if isinstance(label, _LocalElementaryErrorgenLabel):
             label = _GlobalElementaryErrorgenLabel.cast(label, self.sslbls, identity_label=identity_label)
 
@@ -774,6 +796,13 @@ def label_index(self, label, ok_if_missing=False, identity_label='I'):
             raise ValueError("Invalid elementary errorgen type: %s" % str(eetype))
 
         return base + indices[label]
+        return base + indices[elemgen_label]
+        '''
+
+    #@property
+    #def sslbls(self):
+    #    """ The support of this errorgen space, e.g., the qubits where its elements may be nontrivial """
+    #    return self.sslbls
 
     def create_subbasis(self, sslbl_overlap, retain_max_weights=True):
         """
@@ -830,3 +859,125 @@ def difference(self, other_basis):
             `ElementaryErrorgenBasis` to construct the difference with.
         """
         return self.to_explicit_basis().difference(other_basis)
+
+
+#OLD - maybe not needed?
+#class LowWeightElementaryErrorgenBasis(ElementaryErrorgenBasis):
+#    """
+#    Spanned by the elementary error generators of given type(s) (e.g. "Hamiltonian" and/or "other")
+#    and with elements corresponding to a `Basis`, usually of Paulis.
+#    """
+#
+#    def __init__(self, basis_1q, state_space, other_mode, max_ham_weight=None, max_other_weight=None,
+#                 must_overlap_with_these_sslbls=None):
+#        self._basis_1q = basis_1q
+#        self._other_mode = other_mode
+#        self.state_space = state_space
+#        self._max_ham_weight = max_ham_weight
+#        self._max_other_weight = max_other_weight
+#        self._must_overlap_with_these_sslbls = must_overlap_with_these_sslbls
+#
+#        assert(self.state_space.is_entirely_qubits), "FOGI only works for models containing just qubits (so far)"
+#        sslbls = self.state_space.sole_tensor_product_block_labels  # all the model's state space labels
+#        self.sslbls = sslbls  # the "support" of this space - the qubit labels
+#
+#        self._cached_label_indices = None
+#        self._cached_labels_by_support = None
+#        self._cached_elements = None
+#
+#        #Needed?
+#        # self.dim = len(self.labels)  # TODO - update this so we don't always need to build labels
+#        # # (this defeats lazy building via property below) - we can just compute this, especially if
+#        # # not too fancy
+#
+#    @property
+#    def labels(self):
+#        if self._cached_label_indices is None:
+#
+#            def _basis_el_strs(possible_bels, wt):
+#                for els in _itertools.product(*([possible_bels] * wt)):
+#                    yield ''.join(els)
+#
+#            labels = {}
+#            all_bels = self.basis_1q.labels[1:]  # assume first element is identity
+#            nontrivial_bels = self.basis_1q.labels[1:]  # assume first element is identity
+#
+#            max_weight = self._max_ham_weight if (self._max_ham_weight is not None) else len(self.sslbls)
+#            for weight in range(1, max_weight + 1):
+#                for support in _itertools.combinations(self.sslbls, weight):
+#                    if (self._must_overlap_with_these_sslbls is not None
+#                       and len(self._must_overlap_with_these_sslbls.intersection(support)) == 0):
+#                        continue
+#                    if support not in labels: labels[support] = []  # always True?
+#                    labels[support].extend([('H', bel) for bel in _basis_el_strs(nontrivial_bels, weight)])
+#
+#            max_weight = self._max_other_weight if (self._max_other_weight is not None) else len(self.sslbls)
+#            if self._other_mode != "all":
+#                for weight in range(1, max_weight + 1):
+#                    for support in _itertools.combinations(self.sslbls, weight):
+#                        if (self._must_overlap_with_these_sslbls is not None
+#                           and len(self._must_overlap_with_these_sslbls.intersection(support)) == 0):
+#                            continue
+#                        if support not in labels: labels[support] = []
+#                        labels[support].extend([('S', bel) for bel in _basis_el_strs(nontrivial_bels, weight)])
+#            else:
+#                #This is messy code that relies on basis labels being single characters -- TODO improve(?)
+#                idle_char = self.basis_1q.labels[1:]  # assume first element is identity
+#                assert(len(idle_char) == 1 and all([len(c) == 1 for c in nontrivial_bels])), \
+#                    "All basis el labels must be single chars for other_mode=='all'!"
+#                for support in _itertools.combinations(self.sslbls, max_weight):
+#                    # Loop over all possible basis elements for this max-weight support, computing the actual support
+#                    # of each one individually and appending it to the appropriate list
+#                    for bel1 in _basis_el_strs(all_bels, max_weight):
+#                        nonidle_indices1 = [i for i in range(max_weight) if bel1[i] != idle_char]
+#                        for bel2 in _basis_el_strs(all_bels, max_weight):
+#                            nonidle_indices2 = [i for i in range(max_weight) if bel2[i] != idle_char]
+#                            nonidle_indices = list(sorted(set(nonidle_indices1) + set(nonidle_indices2)))
+#                            bel1 = ''.join([bel1[i] for i in nonidle_indices])  # trim to actual support
+#                            bel2 = ''.join([bel2[i] for i in nonidle_indices])  # trim to actual support
+#                            actual_support = tuple([support[i] for i in nonidle_indices])
+#
+#                            if (self._must_overlap_with_these_sslbls is not None
+#                               and len(self._must_overlap_with_these_sslbls.intersection(actual_support)) == 0):
+#                                continue
+#
+#                            if actual_support not in labels: labels[actual_support] = []
+#                            labels[actual_support].append(('S', bel1, bel2))
+#
+#            self._cached_labels_by_support = labels
+#            self._cached_label_indices = _collections.OrderedDict(((support_lbl, i) for i, support_lbl in enumerate(
+#                ((support, lbl) for support, lst in labels.items() for lbl in lst))))
+#
+#        return tuple(self._cached_label_indices.keys())
+#
+#    @property
+#    def element_supports_and_matrices(self):
+#        if self._cached_elements is None:
+#            self._cached_elements = tuple(
+#                ((support, _ot.lindblad_error_generator(elemgen_label, self.basis_1q, normalize=True, sparse=False))
+#                 for support, elemgen_label in self.labels))
+#        return self._cached_elements
+#
+#    def element_index(self, label):
+#        """
+#        TODO: docstring
+#        """
+#        if self._cached_label_indices is None:
+#            self.labels  # triggers building of labels
+#        return self._cached_label_indices[label]
+#
+#    @property
+#    def sslbls(self):
+#        """ The support of this errorgen space, e.g., the qubits where its elements may be nontrivial """
+#        return self.sslbls
+#
+#    def create_subbasis(self, must_overlap_with_these_sslbls, retain_max_weights=True):
+#        """
+#        Create a sub-basis of this basis by including only the elements
+#        that overlap the given support (state space labels)
+#        """
+#        #Note: can we reduce self.state_space?
+#        return CompleteErrorgenBasis(self._basis_1q, self.state_space, self._other_mode,
+#                                     self._max_ham_weight if retain_max_weights else None,
+#                                     self._max_other_weight if retain_max_weights else None,
+#                                     self._must_overlap_with_these_sslbls)
\ No newline at end of file
diff --git a/pygsti/baseobjs/errorgenspace.py b/pygsti/baseobjs/errorgenspace.py
index 0c1b062c9..407864d90 100644
--- a/pygsti/baseobjs/errorgenspace.py
+++ b/pygsti/baseobjs/errorgenspace.py
@@ -11,11 +11,11 @@
 #***************************************************************************************************
 
 import numpy as _np
-
+from pygsti.baseobjs.nicelyserializable import NicelySerializable as _NicelySerializable
 from pygsti.tools import matrixtools as _mt
+from pygsti.baseobjs.errorgenbasis import ExplicitElementaryErrorgenBasis
 
-
-class ErrorgenSpace(object):
+class ErrorgenSpace(_NicelySerializable):
     """
     A vector space of error generators, spanned by some basis.
 
@@ -24,13 +24,25 @@ class ErrorgenSpace(object):
     """
 
     def __init__(self, vectors, basis):
+        super().__init__()
         self.vectors = vectors
         self.elemgen_basis = basis
         #Question: have multiple bases or a single one?
         #self._vectors = [] if (items is None) else items  # list of (basis, vectors_mx) pairs
         # map sslbls => (vectors, basis) where basis.sslbls == sslbls
         # or basis => vectors if bases can hash well(?)
-
+    def _to_nice_serialization(self):
+        state = super()._to_nice_serialization()
+        state.update({'vectors' : self._encodemx(self.vectors),
+                      'basis': self.elemgen_basis._to_nice_serialization()
+
+        }
+
+        )
+        return state
+    @classmethod
+    def from_nice_serialization(cls, state):
+        return cls(state['vectors'], ExplicitElementaryErrorgenBasis.from_nice_serialization(state['basis']))
     def intersection(self, other_space, free_on_unspecified_space=False, use_nice_nullspace=False):
         """
         TODO: docstring
@@ -97,3 +109,13 @@ def normalize(self, norm_order=2):
         for j in range(self.vectors.shape[1]):
             sign = +1 if max(self.vectors[:, j]) >= -min(self.vectors[:, j]) else -1
             self.vectors[:, j] /= sign * _np.linalg.norm(self.vectors[:, j], ord=norm_order)
+
+
+#class LowWeightErrorgenSpace(ErrorgenSpace):
+#    """
+#    Like a SimpleErrorgenSpace but spanned by only the elementary error generators corresponding to
+#    low-weight (up to some maximum weight) basis elements
+#    (so far, only Pauli-product bases work for this, since `Basis` objects don't keep track of each
+#    element's weight (?)).
+#    """
+#    pass
diff --git a/pygsti/baseobjs/statespace.py b/pygsti/baseobjs/statespace.py
index e261daf07..873eb7f9d 100644
--- a/pygsti/baseobjs/statespace.py
+++ b/pygsti/baseobjs/statespace.py
@@ -418,7 +418,7 @@ def create_subspace(self, labels):
         StateSpace
         """
         # Default, generic, implementation constructs an explicit state space
-        labels = set(labels)
+        labels = sorted(set(labels))
         sub_tpb_labels = []
         sub_tpb_udims = []
         sub_tpb_types = []
@@ -1388,7 +1388,6 @@ def __str__(self):
             ['*'.join(["%s(%d%s)" % (lbl, self.label_dims[lbl], 'c' if (self.label_types[lbl] == 'C') else '')
                        for lbl in tpb]) for tpb in self._labels])
 
-
 def default_space_for_dim(dim):
     """
     Create a state space for a given superoperator dimension.
diff --git a/pygsti/models/explicitmodel.py b/pygsti/models/explicitmodel.py
index 01eed989b..179145c1f 100644
--- a/pygsti/models/explicitmodel.py
+++ b/pygsti/models/explicitmodel.py
@@ -23,6 +23,7 @@
 from pygsti.models.memberdict import OrderedMemberDict as _OrderedMemberDict
 from pygsti.models.layerrules import LayerRules as _LayerRules
 from pygsti.models.modelparaminterposer import ModelParamsInterposer as _ModelParamsInterposer
+from pygsti.models.fogistore import FirstOrderGaugeInvariantStore as _FirstOrderGaugeInvariantStore
 from pygsti.models.gaugegroup import GaugeGroup as _GaugeGroup
 from pygsti.forwardsims.forwardsim import ForwardSimulator as _FSim
 from pygsti.forwardsims import matrixforwardsim as _matrixfwdsim
@@ -186,6 +187,16 @@ def _excalc(self):
         return _explicitcalc.ExplicitOpModelCalc(self.state_space.dim, simplified_preps, simplified_ops,
                                                  simplified_effects, self.num_params, self._param_interposer)
 
+    #Unneeded - just use string processing & rely on effect labels *not* having underscores in them
+    #def simplify_spamtuple_to_outcome_label(self, simplified_spamTuple):
+    #    #TODO: make this more efficient (prep lbl isn't even used!)
+    #    for prep_lbl in self.preps:
+    #        for povm_lbl in self.povms:
+    #            for elbl in self.povms[povm_lbl]:
+    #                if simplified_spamTuple == (prep_lbl, povm_lbl + "_" + elbl):
+    #                    return (elbl,) # outcome "label" (a tuple)
+    #    raise ValueError("No outcome label found for simplified spam_tuple: ", simplified_spamTuple)
+
     def _embed_operation(self, op_target_labels, op_val, force=False):
         """
         Called by OrderedMemberDict._auto_embed to create an embedded-gate
@@ -906,7 +917,7 @@ def __str__(self):
         for lbl, inst in self.instruments.items():
             s += "%s = " % str(lbl) + str(inst) + "\n"
         for lbl, factory in self.factories.items():
-            s += "%s = (factory)" % str(lbl) + '\n'
+            s += "%s = (factory)" % lbl + '\n'
         s += "\n"
 
         return s
@@ -1603,7 +1614,9 @@ def _to_nice_serialization(self):
                       'default_gauge_group': (self.default_gauge_group.to_nice_serialization()
                                               if (self.default_gauge_group is not None) else None),
                       'parameter_interposer': (self._param_interposer.to_nice_serialization()
-                                               if (self._param_interposer is not None) else None)
+                                               if (self._param_interposer is not None) else None),
+                      'fogi_store': (self.fogi_store.to_nice_serialization()
+                                               if (self.fogi_store is not None) else None)
                       })
 
         mmgraph = self.create_modelmember_graph()
@@ -1619,6 +1632,8 @@ def _from_nice_serialization(cls, state):
             if (state['default_gauge_group'] is not None) else None
         param_interposer = _ModelParamsInterposer.from_nice_serialization(state['parameter_interposer']) \
             if (state['parameter_interposer'] is not None) else None
+        fogi_store = _FirstOrderGaugeInvariantStore.from_nice_serialization(state['fogi_store']) \
+            if (state['fogi_store'] is not None) else None
         param_labels = state.get('parameter_labels', None)
         param_bounds = state.get('parameter_bounds', None)
 
@@ -1637,6 +1652,7 @@ def _from_nice_serialization(cls, state):
         mdl._clean_paramvec()
         mdl.default_gauge_group = default_gauge_group
         mdl.param_interposer = param_interposer
+        mdl.fogi_store = fogi_store
 
         Np = len(mdl._paramlbls)  # _clean_paramvec sets up ._paramlbls so its length == # of params
         if param_labels and len(param_labels) == Np:
diff --git a/pygsti/models/fogistore.py b/pygsti/models/fogistore.py
index 7dc3f8381..de826fadf 100644
--- a/pygsti/models/fogistore.py
+++ b/pygsti/models/fogistore.py
@@ -21,23 +21,45 @@
 from pygsti.tools import optools as _ot
 from pygsti.tools import slicetools as _slct
 from pygsti.tools import fogitools as _fogit
-
-
-class FirstOrderGaugeInvariantStore(object):
+from pygsti.baseobjs.nicelyserializable import NicelySerializable as _NicelySerializable
+from pygsti.baseobjs.label import Label as _Label
+from pygsti.baseobjs.errorgenlabel import GlobalElementaryErrorgenLabel as _GlobalElementaryErrorgenLabel
+from pygsti.tools.slicetools import slice_hash as _slice_hash
+from pygsti.baseobjs.errorgenspace import ErrorgenSpace as _ErrorgenSpace
+class FirstOrderGaugeInvariantStore(_NicelySerializable):
     """
     An object that computes and stores the first-order-gauge-invariant quantities of a model.
 
     Currently, it is only compatible with :class:`ExplicitOpModel` objects.
     """
+    def __init__(self, primitive_op_labels, gauge_space, elem_errorgen_labels_by_op, op_errorgen_indices, fogi_directions, fogi_metadata, dependent_dir_indices, fogv_directions, allop_gauge_action, gauge_space_directions, norm_order=None, dependent_fogi_action='drop'):
+        super().__init__()
+        self.primitive_op_labels = primitive_op_labels
+        self.gauge_space = gauge_space
+        self.elem_errorgen_labels_by_op = elem_errorgen_labels_by_op
+        self.op_errorgen_indices = op_errorgen_indices
+        self.fogi_directions = fogi_directions
+        self.fogi_metadata = fogi_metadata
+        self.dependent_dir_indices = dependent_dir_indices
+        self.fogv_directions = fogv_directions
+        self.allop_gauge_action = allop_gauge_action
+        self.gauge_space_directions = gauge_space_directions
+        self.norm_order = norm_order
+        self.dependent_fogi_action = dependent_fogi_action
 
-    def __init__(self, gauge_action_matrices_by_op, gauge_action_gauge_spaces_by_op, errorgen_coefficient_labels_by_op,
+        self.errorgen_space_op_elem_labels = tuple([(op_label, elem_lbl) for op_label in self.primitive_op_labels
+                                                    for elem_lbl in self.elem_errorgen_labels_by_op[op_label]])
+        self.fogv_labels = ["%d gauge action" % i for i in range(self.fogv_directions.shape[1])]
+    
+    @classmethod
+    def compute_fogis(cls, gauge_action_matrices_by_op, gauge_action_gauge_spaces_by_op, errorgen_coefficient_labels_by_op,
                  op_label_abbrevs=None, reduce_to_model_space=True,
                  dependent_fogi_action='drop', norm_order=None):
         """
         TODO: docstring
         """
 
-        self.primitive_op_labels = tuple(gauge_action_matrices_by_op.keys())
+        primitive_op_labels = tuple(gauge_action_matrices_by_op.keys())
 
         # Construct common gauge space by special way of intersecting the gauge spaces for all the ops
         # Note: the gauge_space of each op is constructed (see `setup_fogi`) so that the gauge action is
@@ -47,6 +69,7 @@ def __init__(self, gauge_action_matrices_by_op, gauge_action_gauge_spaces_by_op,
         #  *difference* K - UKU^dag in the Lindblad mapping under gauge transform exp(K),  L -> L + K - UKU^dag)
         common_gauge_space = None
         for op_label, gauge_space in gauge_action_gauge_spaces_by_op.items():
+            
             #FOGI DEBUG print("DEBUG gauge space of ", op_label, "has dim", gauge_space.vectors.shape[1])
             if common_gauge_space is None:
                 common_gauge_space = gauge_space
@@ -58,20 +81,20 @@ def __init__(self, gauge_action_matrices_by_op, gauge_action_gauge_spaces_by_op,
         # column space of self.fogi_directions
         #FOGI DEBUG print("DEBUG common gauge space of has dim", common_gauge_space.vectors.shape[1])
         common_gauge_space.normalize()
-        self.gauge_space = common_gauge_space
+        gauge_space = common_gauge_space
 
         # row space of self.fogi_directions - "errgen-set space" lookups
         # -- maybe make this into an "ErrgenSetSpace" object in FUTURE?
-        self.elem_errorgen_labels_by_op = errorgen_coefficient_labels_by_op
+        elem_errorgen_labels_by_op = errorgen_coefficient_labels_by_op
 
-        self.op_errorgen_indices = _fogit._create_op_errgen_indices_dict(self.primitive_op_labels,
-                                                                         self.elem_errorgen_labels_by_op)
-        self.errorgen_space_op_elem_labels = tuple([(op_label, elem_lbl) for op_label in self.primitive_op_labels
-                                                    for elem_lbl in self.elem_errorgen_labels_by_op[op_label]])
+        op_errorgen_indices = _fogit._create_op_errgen_indices_dict(primitive_op_labels,
+                                                                         elem_errorgen_labels_by_op)
+        errorgen_space_op_elem_labels = tuple([(op_label, elem_lbl) for op_label in primitive_op_labels
+                                                    for elem_lbl in elem_errorgen_labels_by_op[op_label]])
         # above is same as flattened self.elem_errorgen_labels_by_op - labels final "row basis" of fogi dirs
 
-        num_elem_errgens = sum([len(labels) for labels in self.elem_errorgen_labels_by_op.values()])
-        allop_gauge_action = _sps.lil_matrix((num_elem_errgens, self.gauge_space.vectors.shape[1]), dtype=complex)
+        num_elem_errgens = sum([len(labels) for labels in elem_errorgen_labels_by_op.values()])
+        allop_gauge_action = _sps.lil_matrix((num_elem_errgens, gauge_space.vectors.shape[1]), dtype=complex)
 
         # Now update (restrict) each op's gauge_action to use common gauge space
         # - need to write the vectors of the common (final) gauge space, w_i, as linear combos of
@@ -81,6 +104,7 @@ def __init__(self, gauge_action_matrices_by_op, gauge_action_gauge_spaces_by_op,
         #  (these could be seen as staring in the union of the common gauge's and op's elemgen
         #   bases - which would just be the common gauge's elemgen basis since it's strictly larger -
         #   restricted to the op's elemben basis)
+        
         for op_label, orig_gauge_space in gauge_action_gauge_spaces_by_op.items():
             #FOGI DEBUG print("DEBUG: ", op_label, orig_gauge_space.vectors.shape, len(orig_gauge_space.elemgen_basis))
             gauge_action = gauge_action_matrices_by_op[op_label]  # a sparse matrix
@@ -93,7 +117,7 @@ def __init__(self, gauge_action_matrices_by_op, gauge_action_gauge_spaces_by_op,
 
             # update gauge action to use common gauge space
             sparse_gauge_action = gauge_action.dot(alpha)
-            allop_gauge_action[self.op_errorgen_indices[op_label], :] = sparse_gauge_action[:, :]
+            allop_gauge_action[op_errorgen_indices[op_label], :] = sparse_gauge_action[:, :]
             gauge_action_matrices_by_op[op_label] = sparse_gauge_action.toarray()  # make **DENSE** here
             # Hopefully matrices aren't too large after this reduction and dense matrices are ok,
             # otherwise we need to change downstream nullspace and pseduoinverse operations to be sparse compatible.
@@ -101,39 +125,40 @@ def __init__(self, gauge_action_matrices_by_op, gauge_action_gauge_spaces_by_op,
             #FUTURE: if update above creates zero-rows in gauge action matrix, maybe remove
             # these rows from the row basis, i.e. self.elem_errorgen_labels_by_op[op_label]
 
-        self.gauge_action_for_op = gauge_action_matrices_by_op
+
+        gauge_action_for_op = gauge_action_matrices_by_op
+
 
         (indep_fogi_directions, indep_fogi_metadata, dep_fogi_directions, dep_fogi_metadata) = \
-            _fogit.construct_fogi_quantities(self.primitive_op_labels, self.gauge_action_for_op,
-                                             self.elem_errorgen_labels_by_op, self.op_errorgen_indices,
-                                             self.gauge_space, op_label_abbrevs, dependent_fogi_action, norm_order)
-        self.fogi_directions = _sps.hstack((indep_fogi_directions, dep_fogi_directions))
-        self.fogi_metadata = indep_fogi_metadata + dep_fogi_metadata  # list concatenation
-        self.dependent_dir_indices = _np.arange(len(indep_fogi_metadata), len(self.fogi_metadata))
-        for j, meta in enumerate(self.fogi_metadata):
-            meta['raw'] = _fogit.op_elem_vec_name(self.fogi_directions[:, j], self.errorgen_space_op_elem_labels,
+            _fogit.construct_fogi_quantities(primitive_op_labels, gauge_action_for_op,
+                                             elem_errorgen_labels_by_op, op_errorgen_indices,
+                                             gauge_space, op_label_abbrevs, dependent_fogi_action, norm_order)
+        fogi_directions = _sps.hstack((indep_fogi_directions, dep_fogi_directions))
+        fogi_metadata = indep_fogi_metadata + dep_fogi_metadata  # list concatenation
+        dependent_dir_indices = _np.arange(len(indep_fogi_metadata), len(fogi_metadata))
+        for j, meta in enumerate(fogi_metadata):
+            meta['raw'] = _fogit.op_elem_vec_name(fogi_directions[:, j], errorgen_space_op_elem_labels,
                                                   op_label_abbrevs if (op_label_abbrevs is not None) else {})
 
-        assert(len(self.errorgen_space_op_elem_labels) == self.fogi_directions.shape[0])
+        assert(len(errorgen_space_op_elem_labels) == fogi_directions.shape[0])
 
         # Note: currently PUNT on doing below with sparse math, as a sparse nullspace routine is unavailable
 
         #First order gauge *variant* directions (the complement of FOGI directions in errgen set space)
         # (directions in errorgen space that correspond to gauge transformations -- to first order)
         #fogv_directions = _mt.nice_nullspace(self.fogi_directions.T)  # can be dependent!
-        self.fogv_directions = _mt.nullspace(self.fogi_directions.toarray().T)  # complement of fogi directions
-        self.fogv_directions = _sps.csc_matrix(self.fogv_directions)  # as though we used sparse math above
-        self.fogv_labels = ["%d gauge action" % i for i in range(self.fogv_directions.shape[1])]
+        fogv_directions = _mt.nullspace(fogi_directions.toarray().T)  # complement of fogi directions
+        fogv_directions = _sps.csc_matrix(fogv_directions)  # as though we used sparse math above
         #self.fogv_labels = ["%s gauge action" % nm
         #                    for nm in _fogit.elem_vec_names(gauge_space_directions, gauge_elemgen_labels)]
 
         #Get gauge-space directions corresponding to the fogv directions
         # (pinv_allop_gauge_action takes errorgen-set -> gauge-gen space)
-        self.allop_gauge_action = allop_gauge_action
-        pinv_allop_gauge_action = _np.linalg.pinv(self.allop_gauge_action.toarray(), rcond=1e-7)
+        allop_gauge_action = allop_gauge_action
+        pinv_allop_gauge_action = _np.linalg.pinv(allop_gauge_action.toarray(), rcond=1e-7)
         gauge_space_directions = _np.dot(pinv_allop_gauge_action,
-                                         self.fogv_directions.toarray())  # in gauge-generator space
-        self.gauge_space_directions = gauge_space_directions
+                                         fogv_directions.toarray())  # in gauge-generator space
+        gauge_space_directions = gauge_space_directions
 
         #Notes on error-gen vs gauge-gen space:
         # self.fogi_directions and self.fogv_directions are dual vectors in error-generator space,
@@ -146,35 +171,86 @@ def __init__(self, gauge_action_matrices_by_op, gauge_action_gauge_spaces_by_op,
         # (or sets of operations).
 
         #Store auxiliary info for later use
-        self.norm_order = norm_order
-        self._dependent_fogi_action = dependent_fogi_action
-
+        fogi_store = cls(primitive_op_labels, gauge_space, elem_errorgen_labels_by_op, op_errorgen_indices, fogi_directions, fogi_metadata, dependent_dir_indices, fogv_directions, allop_gauge_action, gauge_space_directions, norm_order, dependent_fogi_action)
+        fogi_store._check_fogi_store()
+        return fogi_store
         #Assertions to check that everything looks good
-        if True:
-            fogi_dirs = self.fogi_directions.toarray()  # don't bother with sparse math yet
-            fogv_dirs = self.fogv_directions.toarray()
-
-            # We must reduce X_gauge_action to the "in-model gauge space" before testing if the computed vecs are FOGI:
-            assert(_np.linalg.norm(_np.dot(self.allop_gauge_action.toarray().T, fogi_dirs)) < 1e-8)
-
-            #Check that pseudo-inverse was computed correctly (~ matrices are full rank)
-            # fogi_coeffs = dot(fogi_directions.T, elem_errorgen_vec), where elem_errorgen_vec is filled from model
-            #                 params, since fogi_directions columns are *dual* vectors in error-gen space.  Thus,
-            #                 to go in reverse:
-            # elem_errogen_vec = dot(pinv_fogi_dirs_T, fogi_coeffs), where dot(fogi_directions.T, pinv_fogi_dirs_T) == I
-            # (This will only be the case when fogi_vecs are linearly independent, so when dependent_indices == 'drop')
-
-            if dependent_fogi_action == 'drop':
-                #assert(_mt.columns_are_orthogonal(self.fogi_directions))  # not true unless we construct them so...
-                assert(_np.linalg.norm(_np.dot(fogi_dirs.T, _np.linalg.pinv(fogi_dirs.T))
-                                       - _np.identity(fogi_dirs.shape[1], 'd')) < 1e-6)
-
-            # A similar relationship should always hold for the gauge directions, except for these we never
-            #  keep linear dependencies
-            assert(_mt.columns_are_orthogonal(fogv_dirs))
-            assert(_np.linalg.norm(_np.dot(fogv_dirs.T, _np.linalg.pinv(fogv_dirs.T))
-                                   - _np.identity(fogv_dirs.shape[1], 'd')) < 1e-6)
-
+    def _check_fogi_store(self):
+        fogi_dirs = self.fogi_directions.toarray()  # don't bother with sparse math yet
+        fogv_dirs = self.fogv_directions.toarray()
+
+        # We must reduce X_gauge_action to the "in-model gauge space" before testing if the computed vecs are FOGI:
+        assert(_np.linalg.norm(_np.dot(self.allop_gauge_action.toarray().T, fogi_dirs)) < 1e-8)
+
+        #Check that pseudo-inverse was computed correctly (~ matrices are full rank)
+        # fogi_coeffs = dot(fogi_directions.T, elem_errorgen_vec), where elem_errorgen_vec is filled from model
+        #                 params, since fogi_directions columns are *dual* vectors in error-gen space.  Thus,
+        #                 to go in reverse:
+        # elem_errogen_vec = dot(pinv_fogi_dirs_T, fogi_coeffs), where dot(fogi_directions.T, pinv_fogi_dirs_T) == I
+        # (This will only be the case when fogi_vecs are linearly independent, so when dependent_indices == 'drop')
+
+        if self.dependent_fogi_action == 'drop':
+            #assert(_mt.columns_are_orthogonal(self.fogi_directions))  # not true unless we construct them so...
+            assert(_np.linalg.norm(_np.dot(fogi_dirs.T, _np.linalg.pinv(fogi_dirs.T))
+                                    - _np.identity(fogi_dirs.shape[1], 'd')) < 1e-6)
+
+        # A similar relationship should always hold for the gauge directions, except for these we never
+        #  keep linear dependencies
+        assert(_mt.columns_are_orthogonal(fogv_dirs))
+        assert(_np.linalg.norm(_np.dot(fogv_dirs.T, _np.linalg.pinv(fogv_dirs.T))
+                                - _np.identity(fogv_dirs.shape[1], 'd')) < 1e-6)
+    def _to_nice_serialization(self):
+        state = super()._to_nice_serialization()
+        ops_labels = self.primitive_op_labels
+        encoded_metadata = []
+        for object in self.fogi_metadata:
+            encoded_gauge_space = self._encodemx(object['gaugespace_dir'])
+            native_opset_labels = [label.to_native() for label in object['opset']]
+            object_copy = object.copy()
+            object_copy['gaugespace_dir'] = encoded_gauge_space
+            object_copy['opset'] = native_opset_labels
+            encoded_metadata.append(object_copy)
+        
+        state.update({'ops' : [op.to_native() for op in ops_labels],
+                      'elem_errorgen_labels_by_op_keys' : [op.to_native() for op in self.elem_errorgen_labels_by_op.keys()],
+                      'elem_errorgen_labels_by_op_values' : [([label.__str__() for label in errgen_labels]) for errgen_labels in self.elem_errorgen_labels_by_op.values()],
+                      'op_errorgen_indices_keys' : [op.to_native() for op in self.op_errorgen_indices.keys()],
+                      'op_errorgen_indices_values' : [_slice_hash(slice) for slice in self.op_errorgen_indices.values()],
+                      'fogv_directions' : self._encodemx(self.fogv_directions),
+                      'fogi_directions' : self._encodemx(self.fogi_directions),
+                      'fogi_metadata' : encoded_metadata,
+                      'dependent_dir_indices' : self._encodemx(self.dependent_dir_indices),
+                      'allop_gauge_action' : self._encodemx(self.allop_gauge_action),
+                      'gauge_space_directions' : self._encodemx(self.gauge_space_directions),
+                      'norm_order' : self.norm_order,
+                      'dependent_fogi_action' : self.dependent_fogi_action,
+                      'gauge_space' : self.gauge_space._to_nice_serialization()
+                      })
+        return state
+    @classmethod
+    def _from_nice_serialization(cls, state):  
+        primitive_op_labels = [_Label(op) for op in state['ops']]
+        elem_errorgen_labels_by_op_keys = [_Label(op) for op in state['elem_errorgen_labels_by_op_keys']]
+        elem_errorgen_labels_by_op_values = tuple(tuple(_GlobalElementaryErrorgenLabel.cast(label) for label in errgen_labels) for errgen_labels in state['elem_errorgen_labels_by_op_values'])
+        elem_errorgen_labels_by_op = {op_label : errgen_labels for op_label, errgen_labels in zip(elem_errorgen_labels_by_op_keys, elem_errorgen_labels_by_op_values)}
+        op_errorgen_indices_keys = [_Label(op) for op in state['op_errorgen_indices_keys']]
+        op_errorgen_indices =  {op_label :slice(*slice_args) for op_label, slice_args in zip(op_errorgen_indices_keys, state['op_errorgen_indices_values']) }
+        fogv_directions = cls._decodemx(state['fogv_directions'])
+        fogi_directions = cls._decodemx(state['fogi_directions'])
+        
+        for object in state['fogi_metadata']:
+            object['gaugespace_dir'] = cls._decodemx(object['gaugespace_dir'])
+            object['opset'] = tuple(_Label(label) for label in object['opset'])
+        fogi_metadata = state['fogi_metadata']
+        dependent_dir_indices = cls._decodemx(state['dependent_dir_indices'])
+        allop_gauge_action = cls._decodemx(state['allop_gauge_action'])
+        gauge_space_directions = cls._decodemx(state['gauge_space_directions'])
+        dependent_fogi_action = state['dependent_fogi_action']
+        gauge_space = _ErrorgenSpace.from_nice_serialization(state['gauge_space'])
+        norm_order = state['norm_order']
+
+        return cls(primitive_op_labels, gauge_space, elem_errorgen_labels_by_op, op_errorgen_indices, fogi_directions, fogi_metadata, dependent_dir_indices, fogv_directions, allop_gauge_action, gauge_space_directions, norm_order, dependent_fogi_action)
+    
     def find_nice_fogiv_directions(self):
         # BELOW: an attempt to find nice FOGV directions - but we'd like all the vecs to be
         #  orthogonal and this seems to interfere with that, so we'll just leave the fogv dirs messy for now.
diff --git a/pygsti/models/model.py b/pygsti/models/model.py
index 50e89760e..53ae2aa0c 100644
--- a/pygsti/models/model.py
+++ b/pygsti/models/model.py
@@ -657,56 +657,6 @@ def num_modeltest_params(self):
         self._clean_paramvec()
         return Model.num_modeltest_params.fget(self)
 
-    @property
-    def parameter_labels(self):
-        """
-        A list of labels, usually of the form `(op_label, string_description)` describing this model's parameters.
-        """
-        self._clean_paramvec()
-        return self._ops_paramlbls_to_model_paramlbls(self._paramlbls)
-    
-    def set_parameter_label(self, index, label):
-        """
-        Set the label of a single model parameter.
-
-        Parameters
-        ----------
-        index : int
-            The index of the paramter whose label should be set.
-
-        label : object
-            An object that serves to label this parameter.  Often a string.
-
-        Returns
-        -------
-        None
-        """
-        self._clean_paramvec()
-        self._paramlbls[index] = label
-    
-    @property
-    def parameter_bounds(self):
-        """ Upper and lower bounds on the values of each parameter, utilized by optimization routines """
-        self._clean_paramvec()
-        return self._param_bounds
-    
-    @property
-    def num_modeltest_params(self):
-        """
-        The parameter count to use when testing this model against data.
-
-        Often times, this is the same as :meth:`num_params`, but there are times
-        when it can convenient or necessary to use a parameter count different than
-        the actual number of parameters in this model.
-
-        Returns
-        -------
-        int
-            the number of model parameters.
-        """
-        self._clean_paramvec()
-        return Model.num_modeltest_params.fget(self)
-
     def _iter_parameterized_objs(self):
         raise NotImplementedError("Derived Model classes should implement _iter_parameterized_objs")
         #return # default is to have no parameterized objects
@@ -1187,21 +1137,20 @@ def _build_index_mm_map(self):
 
         #Mapping between the model index and the corresponding model members will be more complicated
         #when there is a parameter interposer, so table implementing this for that case.
-        if self.param_interposer is not None:
-            self._index_mm_map = None
-            self._index_mm_label_map = None
-        else:
-            index_mm_map = [[] for _ in range(len(self._paramvec))]
-            index_mm_label_map = [[] for _ in range(len(self._paramvec))]
-            
-            for lbl, obj in self._iter_parameterized_objs():
-                #if the gpindices are a slice then convert to a list of indices.
-                gpindices = _slct.indices(obj.gpindices) if isinstance(obj.gpindices, slice) else obj.gpindices
-                for gpidx in gpindices:
-                    index_mm_map[gpidx].append(obj)
-                    index_mm_label_map[gpidx].append(lbl)
-            self._index_mm_map = index_mm_map
-            self._index_mm_label_map = index_mm_label_map
+        
+
+        ops_param_vec = self._model_paramvec_to_ops_paramvec(self._paramvec)
+        index_mm_map = [[] for _ in range(len(ops_param_vec))]
+        index_mm_label_map = [[] for _ in range(len(ops_param_vec))]
+        
+        for lbl, obj in self._iter_parameterized_objs():
+            #if the gpindices are a slice then convert to a list of indices.
+            gpindices = _slct.indices(obj.gpindices) if isinstance(obj.gpindices, slice) else obj.gpindices
+            for gpidx in gpindices:
+                index_mm_map[gpidx].append(obj)
+                index_mm_label_map[gpidx].append(lbl)
+        self._index_mm_map = index_mm_map
+        self._index_mm_label_map = index_mm_label_map
         #Note to future selves. If we add a flag indicating the presence of collected parameters
         #then we can improve the performance of this by using a simpler structure when no collected
 
@@ -1304,53 +1253,75 @@ def set_parameter_values(self, indices, values, close=False):
         -------
         None
         """
+        orig_param_vec = self._paramvec.copy()
 
         if isinstance(indices[0], tuple):
             #parse the strings into integer indices.
             param_labels_list = self.parameter_labels.tolist()
             indices = [param_labels_list.index(lbl) for lbl in indices]
-                        
+            
+
         for idx, val in zip(indices, values):
             self._paramvec[idx] = val
 
-        if self._param_interposer is not None or self._index_mm_map is None:
-            #fall back to standard from_vector call.
+        if self._index_mm_map is None:
             self.from_vector(self._paramvec)
+            #print("optimized code was skipped")
+            return
+
+        if self._param_interposer is not None:
+            
+            original_errgen_vec = self._param_interposer.transform_matrix @ orig_param_vec
+            new_errgen_vec = self._param_interposer.transform_matrix @ self._paramvec
+            diff_vec =  original_errgen_vec -  new_errgen_vec
+            diff_vec[_np.abs(diff_vec) < 1e-14] = 0
+            non_zero_errgens = _np.nonzero(diff_vec)
+
+            indices = non_zero_errgens[0]
+            values = new_errgen_vec[indices]
+            vec_to_access = new_errgen_vec
         else:
-            #get all of the model members which need to be be updated and loop through them to update their
-            #parameters.
-            unique_mms = {lbl:val for idx in indices for lbl, val in zip(self._index_mm_label_map[idx], self._index_mm_map[idx])}
+            vec_to_access = self._paramvec.copy()
+
+        #get all of the model members which need to be be updated and loop through them to update their
+        #parameters.
+        #test_model = self.copy()
+        #test_model.from_vector(self._paramvec)
+        unique_mms = {lbl:val for idx in indices for lbl, val in zip(self._index_mm_label_map[idx], self._index_mm_map[idx])}
+        for obj in unique_mms.values():
+            obj.from_vector(vec_to_access[obj.gpindices].copy(), close, dirty_value=False)
+        
+        #go through the model members which have been updated and identify whether any of them have children
+        #which may be present in the _opcaches which have already been updated by the parents. I think the
+        #conditions under which this should be safe are: a) the layer rules are ExplicitLayerRules,
+        #b) The parent is a POVM (it should be safe to assume that POVMs update their children, 
+        #and c) the effect is a child of that POVM.
+        
+        if isinstance(self._layer_rules, _ExplicitLayerRules):
+            updated_children = []
             for obj in unique_mms.values():
-                obj.from_vector(self._paramvec[obj.gpindices].copy(), close, dirty_value=False)
-            
-            #go through the model members which have been updated and identify whether any of them have children
-            #which may be present in the _opcaches which have already been updated by the parents. I think the
-            #conditions under which this should be safe are: a) the layer rules are ExplicitLayerRules,
-            #b) The parent is a POVM (it should be safe to assume that POVMs update their children, 
-            #and c) the effect is a child of that POVM.
-            
-            if isinstance(self._layer_rules, _ExplicitLayerRules):
-                updated_children = []
-                for obj in unique_mms.values():
-                    if isinstance(obj, _POVM):
-                        updated_children.extend(obj.values())
-            else:
-                updated_children = None
+                if isinstance(obj, _POVM):
+                    updated_children.extend(obj.values())
+        else:
+            updated_children = None
 
-            # Call from_vector on elements of the cache
-            if self._call_fromvector_on_cache:
-                #print(f'{self._opcaches=}')
-                for opcache in self._opcaches.values():
-                    for obj in opcache.values():
-                        opcache_elem_gpindices = _slct.indices(obj.gpindices) if isinstance(obj.gpindices, slice) else obj.gpindices
-                        if any([idx in opcache_elem_gpindices for idx in indices]):
-                            #check whether we have already updated this object.
-                            if updated_children is not None and any([child is obj for child in updated_children]):
-                                continue
-                            obj.from_vector(self._paramvec[opcache_elem_gpindices], close, dirty_value=False)
+        # Call from_vector on elements of the cache
+        if self._call_fromvector_on_cache:
+            #print(f'{self._opcaches=}')
+            for opcache in self._opcaches.values():
+                for obj in opcache.values():
+                    opcache_elem_gpindices = _slct.indices(obj.gpindices) if isinstance(obj.gpindices, slice) else obj.gpindices
+                    if any([idx in opcache_elem_gpindices for idx in indices]):
+                        #check whether we have already updated this object.
+                        if updated_children is not None and any([child is obj for child in updated_children]):
+                            continue
+                        obj.from_vector(vec_to_access[obj.gpindices].copy(), close, dirty_value=False)
 
-            if OpModel._pcheck: self._check_paramvec()
+        if OpModel._pcheck: self._check_paramvec()
 
+        #for (_, obj), (_, obj2) in zip(test_model._iter_parameterized_objs(), self._iter_parameterized_objs()):
+        #    assert _np.allclose(obj.to_vector(), obj2.to_vector())
+        #    print('checked')
     @property
     def param_interposer(self):
         return self._param_interposer
@@ -2550,7 +2521,7 @@ def _add_reparameterization(self, primitive_op_labels, fogi_dirs, errgenset_spac
 
     def setup_fogi(self, initial_gauge_basis, create_complete_basis_fn=None,
                    op_label_abbrevs=None, reparameterize=False, reduce_to_model_space=True,
-                   dependent_fogi_action='drop', include_spam=True, primitive_op_labels=None):
+                   dependent_fogi_action='drop', include_spam=True, primitive_op_labels=None, include_fogv = False):
 
         from pygsti.baseobjs.errorgenbasis import CompleteElementaryErrorgenBasis as _CompleteElementaryErrorgenBasis
         from pygsti.baseobjs.errorgenbasis import ExplicitElementaryErrorgenBasis as _ExplicitElementaryErrorgenBasis
@@ -2703,16 +2674,34 @@ def create_complete_basis_fn(target_sslbls):
 
         norm_order = "auto"  # NOTE - should be 1 for normalizing 'S' quantities and 2 for 'H',
         # so 'auto' utilizes intelligence within FOGIStore
-        self.fogi_store = _FOGIStore(gauge_action_matrices, gauge_action_gauge_spaces,
+
+
+        self.fogi_store = _FOGIStore.compute_fogis(gauge_action_matrices, gauge_action_gauge_spaces,
                                      errorgen_coefficient_labels,  # gauge_errgen_space_labels,
                                      op_label_abbrevs, reduce_to_model_space, dependent_fogi_action,
                                      norm_order=norm_order)
-
+        #print(f'{primitive_op_labels=}')
+        #print(f'{primitive_prep_labels=}')
+        #print(f'{primitive_povm_labels=}')
+        #print(f'{self.fogi_store.fogi_directions.toarray()=}')
+        #print(f'{self.fogi_store.errorgen_space_op_elem_labels=}')
         if reparameterize:
-            self.param_interposer = self._add_reparameterization(
+            fogi_interposer = self._add_reparameterization(
                 primitive_op_labels + primitive_prep_labels + primitive_povm_labels,
                 self.fogi_store.fogi_directions.toarray(),  # DENSE now (leave sparse in FUTURE?)
                 self.fogi_store.errorgen_space_op_elem_labels)
+            
+            if include_fogv:
+                fogv_interposer = self._add_reparameterization(
+                primitive_op_labels + primitive_prep_labels + primitive_povm_labels,
+                self.fogi_store.fogv_directions.toarray(),  # DENSE now (leave sparse in FUTURE?)
+                self.fogi_store.errorgen_space_op_elem_labels)
+            
+            self.param_interposer = fogi_interposer
+
+            if include_fogv:
+                self.param_interposer.gaugefull_inv_transform_matrix = _np.vstack([fogi_interposer.inv_transform_matrix, fogv_interposer.inv_transform_matrix])
+                self.param_interposer.gaugefull_transform_matrix = _np.linalg.inv(self.param_interposer.gaugefull_inv_transform_matrix)
 
     def fogi_errorgen_component_labels(self, include_fogv=False, typ='normal'):
         labels = self.fogi_store.fogi_errorgen_direction_labels(typ)

From f13821d6a6f92e2b72ffa22de417768b42033a8c Mon Sep 17 00:00:00 2001
From: Juan G Mendoza <jgonz19@sandia.gov>
Date: Tue, 24 Jun 2025 14:50:26 -0600
Subject: [PATCH 02/34] changed back some changes from develop that were lost

---
 pygsti/baseobjs/errorgenbasis.py | 143 +------------------------------
 pygsti/baseobjs/errorgenspace.py |  12 +--
 pygsti/models/explicitmodel.py   |  14 +--
 pygsti/models/model.py           |  23 +----
 4 files changed, 8 insertions(+), 184 deletions(-)

diff --git a/pygsti/baseobjs/errorgenbasis.py b/pygsti/baseobjs/errorgenbasis.py
index 99333dfe8..24be81bf5 100644
--- a/pygsti/baseobjs/errorgenbasis.py
+++ b/pygsti/baseobjs/errorgenbasis.py
@@ -754,15 +754,7 @@ def label_index(self, label, ok_if_missing=False, identity_label='I'):
         identity_label : str, optional (default 'I')
             An optional string specifying the label used to denote the identity in basis element labels.
         """
-        try:
-            return self.labels.index(label)
-        except ValueError as error:
-            
-            if ok_if_missing:
-                return None
-            else:
-                raise error
-        '''
+
         if isinstance(label, _LocalElementaryErrorgenLabel):
             label = _GlobalElementaryErrorgenLabel.cast(label, self.sslbls, identity_label=identity_label)
 
@@ -796,14 +788,7 @@ def label_index(self, label, ok_if_missing=False, identity_label='I'):
             raise ValueError("Invalid elementary errorgen type: %s" % str(eetype))
 
         return base + indices[label]
-        return base + indices[elemgen_label]
-        '''
-
-    #@property
-    #def sslbls(self):
-    #    """ The support of this errorgen space, e.g., the qubits where its elements may be nontrivial """
-    #    return self.sslbls
-
+        
     def create_subbasis(self, sslbl_overlap, retain_max_weights=True):
         """
         Create a sub-basis of this basis by including only the elements
@@ -858,126 +843,4 @@ def difference(self, other_basis):
         other_basis : `ElementaryErrorgenBasis`
             `ElementaryErrorgenBasis` to construct the difference with.
         """
-        return self.to_explicit_basis().difference(other_basis)
-
-
-#OLD - maybe not needed?
-#class LowWeightElementaryErrorgenBasis(ElementaryErrorgenBasis):
-#    """
-#    Spanned by the elementary error generators of given type(s) (e.g. "Hamiltonian" and/or "other")
-#    and with elements corresponding to a `Basis`, usually of Paulis.
-#    """
-#
-#    def __init__(self, basis_1q, state_space, other_mode, max_ham_weight=None, max_other_weight=None,
-#                 must_overlap_with_these_sslbls=None):
-#        self._basis_1q = basis_1q
-#        self._other_mode = other_mode
-#        self.state_space = state_space
-#        self._max_ham_weight = max_ham_weight
-#        self._max_other_weight = max_other_weight
-#        self._must_overlap_with_these_sslbls = must_overlap_with_these_sslbls
-#
-#        assert(self.state_space.is_entirely_qubits), "FOGI only works for models containing just qubits (so far)"
-#        sslbls = self.state_space.sole_tensor_product_block_labels  # all the model's state space labels
-#        self.sslbls = sslbls  # the "support" of this space - the qubit labels
-#
-#        self._cached_label_indices = None
-#        self._cached_labels_by_support = None
-#        self._cached_elements = None
-#
-#        #Needed?
-#        # self.dim = len(self.labels)  # TODO - update this so we don't always need to build labels
-#        # # (this defeats lazy building via property below) - we can just compute this, especially if
-#        # # not too fancy
-#
-#    @property
-#    def labels(self):
-#        if self._cached_label_indices is None:
-#
-#            def _basis_el_strs(possible_bels, wt):
-#                for els in _itertools.product(*([possible_bels] * wt)):
-#                    yield ''.join(els)
-#
-#            labels = {}
-#            all_bels = self.basis_1q.labels[1:]  # assume first element is identity
-#            nontrivial_bels = self.basis_1q.labels[1:]  # assume first element is identity
-#
-#            max_weight = self._max_ham_weight if (self._max_ham_weight is not None) else len(self.sslbls)
-#            for weight in range(1, max_weight + 1):
-#                for support in _itertools.combinations(self.sslbls, weight):
-#                    if (self._must_overlap_with_these_sslbls is not None
-#                       and len(self._must_overlap_with_these_sslbls.intersection(support)) == 0):
-#                        continue
-#                    if support not in labels: labels[support] = []  # always True?
-#                    labels[support].extend([('H', bel) for bel in _basis_el_strs(nontrivial_bels, weight)])
-#
-#            max_weight = self._max_other_weight if (self._max_other_weight is not None) else len(self.sslbls)
-#            if self._other_mode != "all":
-#                for weight in range(1, max_weight + 1):
-#                    for support in _itertools.combinations(self.sslbls, weight):
-#                        if (self._must_overlap_with_these_sslbls is not None
-#                           and len(self._must_overlap_with_these_sslbls.intersection(support)) == 0):
-#                            continue
-#                        if support not in labels: labels[support] = []
-#                        labels[support].extend([('S', bel) for bel in _basis_el_strs(nontrivial_bels, weight)])
-#            else:
-#                #This is messy code that relies on basis labels being single characters -- TODO improve(?)
-#                idle_char = self.basis_1q.labels[1:]  # assume first element is identity
-#                assert(len(idle_char) == 1 and all([len(c) == 1 for c in nontrivial_bels])), \
-#                    "All basis el labels must be single chars for other_mode=='all'!"
-#                for support in _itertools.combinations(self.sslbls, max_weight):
-#                    # Loop over all possible basis elements for this max-weight support, computing the actual support
-#                    # of each one individually and appending it to the appropriate list
-#                    for bel1 in _basis_el_strs(all_bels, max_weight):
-#                        nonidle_indices1 = [i for i in range(max_weight) if bel1[i] != idle_char]
-#                        for bel2 in _basis_el_strs(all_bels, max_weight):
-#                            nonidle_indices2 = [i for i in range(max_weight) if bel2[i] != idle_char]
-#                            nonidle_indices = list(sorted(set(nonidle_indices1) + set(nonidle_indices2)))
-#                            bel1 = ''.join([bel1[i] for i in nonidle_indices])  # trim to actual support
-#                            bel2 = ''.join([bel2[i] for i in nonidle_indices])  # trim to actual support
-#                            actual_support = tuple([support[i] for i in nonidle_indices])
-#
-#                            if (self._must_overlap_with_these_sslbls is not None
-#                               and len(self._must_overlap_with_these_sslbls.intersection(actual_support)) == 0):
-#                                continue
-#
-#                            if actual_support not in labels: labels[actual_support] = []
-#                            labels[actual_support].append(('S', bel1, bel2))
-#
-#            self._cached_labels_by_support = labels
-#            self._cached_label_indices = _collections.OrderedDict(((support_lbl, i) for i, support_lbl in enumerate(
-#                ((support, lbl) for support, lst in labels.items() for lbl in lst))))
-#
-#        return tuple(self._cached_label_indices.keys())
-#
-#    @property
-#    def element_supports_and_matrices(self):
-#        if self._cached_elements is None:
-#            self._cached_elements = tuple(
-#                ((support, _ot.lindblad_error_generator(elemgen_label, self.basis_1q, normalize=True, sparse=False))
-#                 for support, elemgen_label in self.labels))
-#        return self._cached_elements
-#
-#    def element_index(self, label):
-#        """
-#        TODO: docstring
-#        """
-#        if self._cached_label_indices is None:
-#            self.labels  # triggers building of labels
-#        return self._cached_label_indices[label]
-#
-#    @property
-#    def sslbls(self):
-#        """ The support of this errorgen space, e.g., the qubits where its elements may be nontrivial """
-#        return self.sslbls
-#
-#    def create_subbasis(self, must_overlap_with_these_sslbls, retain_max_weights=True):
-#        """
-#        Create a sub-basis of this basis by including only the elements
-#        that overlap the given support (state space labels)
-#        """
-#        #Note: can we reduce self.state_space?
-#        return CompleteErrorgenBasis(self._basis_1q, self.state_space, self._other_mode,
-#                                     self._max_ham_weight if retain_max_weights else None,
-#                                     self._max_other_weight if retain_max_weights else None,
-#                                     self._must_overlap_with_these_sslbls)
\ No newline at end of file
+        return self.to_explicit_basis().difference(other_basis)
\ No newline at end of file
diff --git a/pygsti/baseobjs/errorgenspace.py b/pygsti/baseobjs/errorgenspace.py
index 407864d90..9ede840e4 100644
--- a/pygsti/baseobjs/errorgenspace.py
+++ b/pygsti/baseobjs/errorgenspace.py
@@ -108,14 +108,4 @@ def normalize(self, norm_order=2):
         """
         for j in range(self.vectors.shape[1]):
             sign = +1 if max(self.vectors[:, j]) >= -min(self.vectors[:, j]) else -1
-            self.vectors[:, j] /= sign * _np.linalg.norm(self.vectors[:, j], ord=norm_order)
-
-
-#class LowWeightErrorgenSpace(ErrorgenSpace):
-#    """
-#    Like a SimpleErrorgenSpace but spanned by only the elementary error generators corresponding to
-#    low-weight (up to some maximum weight) basis elements
-#    (so far, only Pauli-product bases work for this, since `Basis` objects don't keep track of each
-#    element's weight (?)).
-#    """
-#    pass
+            self.vectors[:, j] /= sign * _np.linalg.norm(self.vectors[:, j], ord=norm_order)
\ No newline at end of file
diff --git a/pygsti/models/explicitmodel.py b/pygsti/models/explicitmodel.py
index 179145c1f..95b2314e8 100644
--- a/pygsti/models/explicitmodel.py
+++ b/pygsti/models/explicitmodel.py
@@ -186,17 +186,7 @@ def _excalc(self):
 
         return _explicitcalc.ExplicitOpModelCalc(self.state_space.dim, simplified_preps, simplified_ops,
                                                  simplified_effects, self.num_params, self._param_interposer)
-
-    #Unneeded - just use string processing & rely on effect labels *not* having underscores in them
-    #def simplify_spamtuple_to_outcome_label(self, simplified_spamTuple):
-    #    #TODO: make this more efficient (prep lbl isn't even used!)
-    #    for prep_lbl in self.preps:
-    #        for povm_lbl in self.povms:
-    #            for elbl in self.povms[povm_lbl]:
-    #                if simplified_spamTuple == (prep_lbl, povm_lbl + "_" + elbl):
-    #                    return (elbl,) # outcome "label" (a tuple)
-    #    raise ValueError("No outcome label found for simplified spam_tuple: ", simplified_spamTuple)
-
+    
     def _embed_operation(self, op_target_labels, op_val, force=False):
         """
         Called by OrderedMemberDict._auto_embed to create an embedded-gate
@@ -917,7 +907,7 @@ def __str__(self):
         for lbl, inst in self.instruments.items():
             s += "%s = " % str(lbl) + str(inst) + "\n"
         for lbl, factory in self.factories.items():
-            s += "%s = (factory)" % lbl + '\n'
+            s += "%s = (factory)" % str(lbl) + '\n'
         s += "\n"
 
         return s
diff --git a/pygsti/models/model.py b/pygsti/models/model.py
index 53ae2aa0c..49223cf3f 100644
--- a/pygsti/models/model.py
+++ b/pygsti/models/model.py
@@ -2521,11 +2521,9 @@ def _add_reparameterization(self, primitive_op_labels, fogi_dirs, errgenset_spac
 
     def setup_fogi(self, initial_gauge_basis, create_complete_basis_fn=None,
                    op_label_abbrevs=None, reparameterize=False, reduce_to_model_space=True,
-                   dependent_fogi_action='drop', include_spam=True, primitive_op_labels=None, include_fogv = False):
+                   dependent_fogi_action='drop', include_spam=True, primitive_op_labels=None):
 
         from pygsti.baseobjs.errorgenbasis import CompleteElementaryErrorgenBasis as _CompleteElementaryErrorgenBasis
-        from pygsti.baseobjs.errorgenbasis import ExplicitElementaryErrorgenBasis as _ExplicitElementaryErrorgenBasis
-        from pygsti.baseobjs.errorgenspace import ErrorgenSpace as _ErrorgenSpace
 
         from pygsti.tools import basistools as _bt
         from pygsti.tools import fogitools as _fogit
@@ -2680,29 +2678,12 @@ def create_complete_basis_fn(target_sslbls):
                                      errorgen_coefficient_labels,  # gauge_errgen_space_labels,
                                      op_label_abbrevs, reduce_to_model_space, dependent_fogi_action,
                                      norm_order=norm_order)
-        #print(f'{primitive_op_labels=}')
-        #print(f'{primitive_prep_labels=}')
-        #print(f'{primitive_povm_labels=}')
-        #print(f'{self.fogi_store.fogi_directions.toarray()=}')
-        #print(f'{self.fogi_store.errorgen_space_op_elem_labels=}')
         if reparameterize:
-            fogi_interposer = self._add_reparameterization(
+            self.param_interposer = self._add_reparameterization(
                 primitive_op_labels + primitive_prep_labels + primitive_povm_labels,
                 self.fogi_store.fogi_directions.toarray(),  # DENSE now (leave sparse in FUTURE?)
                 self.fogi_store.errorgen_space_op_elem_labels)
             
-            if include_fogv:
-                fogv_interposer = self._add_reparameterization(
-                primitive_op_labels + primitive_prep_labels + primitive_povm_labels,
-                self.fogi_store.fogv_directions.toarray(),  # DENSE now (leave sparse in FUTURE?)
-                self.fogi_store.errorgen_space_op_elem_labels)
-            
-            self.param_interposer = fogi_interposer
-
-            if include_fogv:
-                self.param_interposer.gaugefull_inv_transform_matrix = _np.vstack([fogi_interposer.inv_transform_matrix, fogv_interposer.inv_transform_matrix])
-                self.param_interposer.gaugefull_transform_matrix = _np.linalg.inv(self.param_interposer.gaugefull_inv_transform_matrix)
-
     def fogi_errorgen_component_labels(self, include_fogv=False, typ='normal'):
         labels = self.fogi_store.fogi_errorgen_direction_labels(typ)
         if include_fogv:

From 67f31f7e6d646af748e5134b3afa6b896f77bbb4 Mon Sep 17 00:00:00 2001
From: Juan G Mendoza <jgonz19@sandia.gov>
Date: Tue, 24 Jun 2025 15:02:58 -0600
Subject: [PATCH 03/34] small white space corrections

---
 pygsti/baseobjs/errorgenbasis.py | 8 +++-----
 pygsti/baseobjs/errorgenspace.py | 5 +----
 pygsti/models/explicitmodel.py   | 2 +-
 pygsti/models/fogistore.py       | 2 +-
 pygsti/models/model.py           | 5 +----
 5 files changed, 7 insertions(+), 15 deletions(-)

diff --git a/pygsti/baseobjs/errorgenbasis.py b/pygsti/baseobjs/errorgenbasis.py
index 24be81bf5..b8b6c90ef 100644
--- a/pygsti/baseobjs/errorgenbasis.py
+++ b/pygsti/baseobjs/errorgenbasis.py
@@ -123,9 +123,7 @@ def _to_nice_serialization(self):
         state.update({'state_space' : self.state_space._to_nice_serialization(),
                       'labels' : [label.__str__() for label in self.labels],
                       '_basis_1q' : self._basis_1q if isinstance(self._basis_1q, str) else self._basis_1q._to_nice_serialization()
-
-        }
-        )
+        })
         return state
     @classmethod
     def from_nice_serialization(cls, state):
@@ -273,6 +271,7 @@ def union(self, other_basis):
         #assert that these two bases have compatible label types.
         msg = 'Incompatible `ElementaryErrrogenLabel` types, the two `ElementaryErrorgenBasis` should have the same label type.'
         assert type(self._labels[0]) == type(other_basis.labels[0]), msg
+
         #Get the union of the two bases labels.
         union_labels = set(self._labels) | set(other_basis.labels)
         union_state_space = self.state_space.union(other_basis.state_space)
@@ -754,7 +753,6 @@ def label_index(self, label, ok_if_missing=False, identity_label='I'):
         identity_label : str, optional (default 'I')
             An optional string specifying the label used to denote the identity in basis element labels.
         """
-
         if isinstance(label, _LocalElementaryErrorgenLabel):
             label = _GlobalElementaryErrorgenLabel.cast(label, self.sslbls, identity_label=identity_label)
 
@@ -788,7 +786,7 @@ def label_index(self, label, ok_if_missing=False, identity_label='I'):
             raise ValueError("Invalid elementary errorgen type: %s" % str(eetype))
 
         return base + indices[label]
-        
+
     def create_subbasis(self, sslbl_overlap, retain_max_weights=True):
         """
         Create a sub-basis of this basis by including only the elements
diff --git a/pygsti/baseobjs/errorgenspace.py b/pygsti/baseobjs/errorgenspace.py
index 9ede840e4..367515bcc 100644
--- a/pygsti/baseobjs/errorgenspace.py
+++ b/pygsti/baseobjs/errorgenspace.py
@@ -35,10 +35,7 @@ def _to_nice_serialization(self):
         state = super()._to_nice_serialization()
         state.update({'vectors' : self._encodemx(self.vectors),
                       'basis': self.elemgen_basis._to_nice_serialization()
-
-        }
-
-        )
+        })
         return state
     @classmethod
     def from_nice_serialization(cls, state):
diff --git a/pygsti/models/explicitmodel.py b/pygsti/models/explicitmodel.py
index 95b2314e8..fef21c987 100644
--- a/pygsti/models/explicitmodel.py
+++ b/pygsti/models/explicitmodel.py
@@ -186,7 +186,7 @@ def _excalc(self):
 
         return _explicitcalc.ExplicitOpModelCalc(self.state_space.dim, simplified_preps, simplified_ops,
                                                  simplified_effects, self.num_params, self._param_interposer)
-    
+
     def _embed_operation(self, op_target_labels, op_val, force=False):
         """
         Called by OrderedMemberDict._auto_embed to create an embedded-gate
diff --git a/pygsti/models/fogistore.py b/pygsti/models/fogistore.py
index de826fadf..2dfd41b61 100644
--- a/pygsti/models/fogistore.py
+++ b/pygsti/models/fogistore.py
@@ -26,6 +26,7 @@
 from pygsti.baseobjs.errorgenlabel import GlobalElementaryErrorgenLabel as _GlobalElementaryErrorgenLabel
 from pygsti.tools.slicetools import slice_hash as _slice_hash
 from pygsti.baseobjs.errorgenspace import ErrorgenSpace as _ErrorgenSpace
+
 class FirstOrderGaugeInvariantStore(_NicelySerializable):
     """
     An object that computes and stores the first-order-gauge-invariant quantities of a model.
@@ -70,7 +71,6 @@ def compute_fogis(cls, gauge_action_matrices_by_op, gauge_action_gauge_spaces_by
         common_gauge_space = None
         for op_label, gauge_space in gauge_action_gauge_spaces_by_op.items():
             
-            #FOGI DEBUG print("DEBUG gauge space of ", op_label, "has dim", gauge_space.vectors.shape[1])
             if common_gauge_space is None:
                 common_gauge_space = gauge_space
             else:
diff --git a/pygsti/models/model.py b/pygsti/models/model.py
index 49223cf3f..1e124f94b 100644
--- a/pygsti/models/model.py
+++ b/pygsti/models/model.py
@@ -1319,9 +1319,6 @@ def set_parameter_values(self, indices, values, close=False):
 
         if OpModel._pcheck: self._check_paramvec()
 
-        #for (_, obj), (_, obj2) in zip(test_model._iter_parameterized_objs(), self._iter_parameterized_objs()):
-        #    assert _np.allclose(obj.to_vector(), obj2.to_vector())
-        #    print('checked')
     @property
     def param_interposer(self):
         return self._param_interposer
@@ -2683,7 +2680,7 @@ def create_complete_basis_fn(target_sslbls):
                 primitive_op_labels + primitive_prep_labels + primitive_povm_labels,
                 self.fogi_store.fogi_directions.toarray(),  # DENSE now (leave sparse in FUTURE?)
                 self.fogi_store.errorgen_space_op_elem_labels)
-            
+
     def fogi_errorgen_component_labels(self, include_fogv=False, typ='normal'):
         labels = self.fogi_store.fogi_errorgen_direction_labels(typ)
         if include_fogv:

From f358b32c63f8bd7d067252aedffa8f703cc8ea35 Mon Sep 17 00:00:00 2001
From: Juan G Mendoza <jgonz19@sandia.gov>
Date: Tue, 24 Jun 2025 16:05:40 -0600
Subject: [PATCH 04/34] fixed fogistore attribute name

---
 pygsti/models/fogistore.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pygsti/models/fogistore.py b/pygsti/models/fogistore.py
index 2dfd41b61..11dbc8d55 100644
--- a/pygsti/models/fogistore.py
+++ b/pygsti/models/fogistore.py
@@ -46,7 +46,7 @@ def __init__(self, primitive_op_labels, gauge_space, elem_errorgen_labels_by_op,
         self.allop_gauge_action = allop_gauge_action
         self.gauge_space_directions = gauge_space_directions
         self.norm_order = norm_order
-        self.dependent_fogi_action = dependent_fogi_action
+        self._dependent_fogi_action = dependent_fogi_action
 
         self.errorgen_space_op_elem_labels = tuple([(op_label, elem_lbl) for op_label in self.primitive_op_labels
                                                     for elem_lbl in self.elem_errorgen_labels_by_op[op_label]])

From 8fc9fa285ea1216022fe0b684234fc29d838310d Mon Sep 17 00:00:00 2001
From: Juan G Mendoza <jgonz19@sandia.gov>
Date: Wed, 25 Jun 2025 10:40:30 -0600
Subject: [PATCH 05/34] fixed other instances of dependent_fogi_action being
 incorrectly named

---
 pygsti/models/fogistore.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/pygsti/models/fogistore.py b/pygsti/models/fogistore.py
index 11dbc8d55..347634644 100644
--- a/pygsti/models/fogistore.py
+++ b/pygsti/models/fogistore.py
@@ -189,7 +189,7 @@ def _check_fogi_store(self):
         # elem_errogen_vec = dot(pinv_fogi_dirs_T, fogi_coeffs), where dot(fogi_directions.T, pinv_fogi_dirs_T) == I
         # (This will only be the case when fogi_vecs are linearly independent, so when dependent_indices == 'drop')
 
-        if self.dependent_fogi_action == 'drop':
+        if self._dependent_fogi_action == 'drop':
             #assert(_mt.columns_are_orthogonal(self.fogi_directions))  # not true unless we construct them so...
             assert(_np.linalg.norm(_np.dot(fogi_dirs.T, _np.linalg.pinv(fogi_dirs.T))
                                     - _np.identity(fogi_dirs.shape[1], 'd')) < 1e-6)
@@ -223,7 +223,7 @@ def _to_nice_serialization(self):
                       'allop_gauge_action' : self._encodemx(self.allop_gauge_action),
                       'gauge_space_directions' : self._encodemx(self.gauge_space_directions),
                       'norm_order' : self.norm_order,
-                      'dependent_fogi_action' : self.dependent_fogi_action,
+                      'dependent_fogi_action' : self._dependent_fogi_action,
                       'gauge_space' : self.gauge_space._to_nice_serialization()
                       })
         return state

From b15ff9b7da876588e47f4e4d9b9826a64efc82f0 Mon Sep 17 00:00:00 2001
From: Juan G Mendoza <jgonz19@sandia.gov>
Date: Wed, 25 Jun 2025 17:13:28 -0600
Subject: [PATCH 06/34] added serialization methods for completerrgenbasis,
 added support for local errgenbasis labels

---
 pygsti/baseobjs/errorgenbasis.py | 28 ++++++++++++++++++++++++++--
 pygsti/models/fogistore.py       |  2 +-
 pygsti/models/model.py           |  2 +-
 3 files changed, 28 insertions(+), 4 deletions(-)

diff --git a/pygsti/baseobjs/errorgenbasis.py b/pygsti/baseobjs/errorgenbasis.py
index b8b6c90ef..3a39af444 100644
--- a/pygsti/baseobjs/errorgenbasis.py
+++ b/pygsti/baseobjs/errorgenbasis.py
@@ -50,6 +50,7 @@ def __len__(self):
         Number of elementary errorgen elements in this basis.
         """
         return len(self.labels)
+    
 
 #helper function for checking label types.
 def _all_elements_same_type(lst):
@@ -122,12 +123,17 @@ def _to_nice_serialization(self):
         state = super()._to_nice_serialization()
         state.update({'state_space' : self.state_space._to_nice_serialization(),
                       'labels' : [label.__str__() for label in self.labels],
+                      'label_type' : 'global' if isinstance(self.labels[0], _GlobalElementaryErrorgenLabel) else 'local',
                       '_basis_1q' : self._basis_1q if isinstance(self._basis_1q, str) else self._basis_1q._to_nice_serialization()
         })
         return state
     @classmethod
     def from_nice_serialization(cls, state):
-        return cls(_StateSpace.from_nice_serialization(state['state_space']), [_GlobalElementaryErrorgenLabel.cast(label) for label in state['labels']], state['_basis_1q'] if isinstance(state['_basis_1q'], str) else _Basis.from_nice_serialization(state['_basis_1q']))
+        if state['label_type'] == 'global':
+            cast = _GlobalElementaryErrorgenLabel.cast
+        else:
+            cast = _LocalElementaryErrorgenLabel.cast
+        return cls(_StateSpace.from_nice_serialization(state['state_space']), [cast(label) for label in state['labels']], state['_basis_1q'] if isinstance(state['_basis_1q'], str) else _Basis.from_nice_serialization(state['_basis_1q']))
     @property
     def labels(self):
         return self._labels
@@ -528,7 +534,7 @@ def __init__(self, basis_1q, state_space, elementary_errorgen_types=('H', 'S', '
             `GlobalElementaryErrorgenLabel` and `LocalElementaryErrorgenLabel`,
             respectively.
         """
-
+        super().__init__()
         if isinstance(basis_1q, _Basis):
             self._basis_1q = basis_1q
         elif isinstance(basis_1q, str):
@@ -601,6 +607,24 @@ def __init__(self, basis_1q, state_space, elementary_errorgen_types=('H', 'S', '
         #                      (IXX,XXI) #   on right, loop over all possible choices of at least one, an at most m,
         #                      (IXX,XXX) #    nontrivial indices to place within the m nontriv left indices (1 & 2 here)
 
+    def _to_nice_serialization(self):
+         
+        state = super()._to_nice_serialization()
+        state.update({'basis_1q' : self._basis_1q.to_nice_serialization(),
+                      'state_space' : self.state_space.to_nice_serialization(),
+                      'elementary_errorgen_types': self._elementary_errorgen_types,
+                      'max_weights' : self.max_weights,
+                      'sslbl_overlap' : self._sslbl_overlap,
+                      'default_label_type': self._default_lbl_typ 
+        })
+        return state
+    
+    @classmethod
+    def from_nice_serialization(cls, state):
+        
+        return cls(_Basis.from_nice_serialization(state['basis_1q']), _StateSpace.from_nice_serialization(state['state_space']), \
+                   state['elementary_errorgen_types'], max_weights=state['max_weights'], sslbl_overlap=state['sslbl_overlap'], \
+                    default_label_type=state['default_label_type'])
     def __len__(self):
         """ Number of elementary errorgen elements in this basis """
         return self._offsets[self._elementary_errorgen_types[-1]]['END']
diff --git a/pygsti/models/fogistore.py b/pygsti/models/fogistore.py
index 347634644..ecdac056e 100644
--- a/pygsti/models/fogistore.py
+++ b/pygsti/models/fogistore.py
@@ -53,7 +53,7 @@ def __init__(self, primitive_op_labels, gauge_space, elem_errorgen_labels_by_op,
         self.fogv_labels = ["%d gauge action" % i for i in range(self.fogv_directions.shape[1])]
     
     @classmethod
-    def compute_fogis(cls, gauge_action_matrices_by_op, gauge_action_gauge_spaces_by_op, errorgen_coefficient_labels_by_op,
+    def from_gauge_action_matrices(cls, gauge_action_matrices_by_op, gauge_action_gauge_spaces_by_op, errorgen_coefficient_labels_by_op,
                  op_label_abbrevs=None, reduce_to_model_space=True,
                  dependent_fogi_action='drop', norm_order=None):
         """
diff --git a/pygsti/models/model.py b/pygsti/models/model.py
index 1e124f94b..88a5b1ba5 100644
--- a/pygsti/models/model.py
+++ b/pygsti/models/model.py
@@ -2671,7 +2671,7 @@ def create_complete_basis_fn(target_sslbls):
         # so 'auto' utilizes intelligence within FOGIStore
 
 
-        self.fogi_store = _FOGIStore.compute_fogis(gauge_action_matrices, gauge_action_gauge_spaces,
+        self.fogi_store = _FOGIStore.from_gauge_action_matrices(gauge_action_matrices, gauge_action_gauge_spaces,
                                      errorgen_coefficient_labels,  # gauge_errgen_space_labels,
                                      op_label_abbrevs, reduce_to_model_space, dependent_fogi_action,
                                      norm_order=norm_order)

From 5b4679e946fe43994f677627ee670336d6f1dc2e Mon Sep 17 00:00:00 2001
From: Juan G Mendoza <jgonz19@sandia.gov>
Date: Wed, 25 Jun 2025 17:14:56 -0600
Subject: [PATCH 07/34] removed print statements

---
 pygsti/models/model.py | 2 --
 1 file changed, 2 deletions(-)

diff --git a/pygsti/models/model.py b/pygsti/models/model.py
index 88a5b1ba5..7ff38dd7c 100644
--- a/pygsti/models/model.py
+++ b/pygsti/models/model.py
@@ -1266,7 +1266,6 @@ def set_parameter_values(self, indices, values, close=False):
 
         if self._index_mm_map is None:
             self.from_vector(self._paramvec)
-            #print("optimized code was skipped")
             return
 
         if self._param_interposer is not None:
@@ -1307,7 +1306,6 @@ def set_parameter_values(self, indices, values, close=False):
 
         # Call from_vector on elements of the cache
         if self._call_fromvector_on_cache:
-            #print(f'{self._opcaches=}')
             for opcache in self._opcaches.values():
                 for obj in opcache.values():
                     opcache_elem_gpindices = _slct.indices(obj.gpindices) if isinstance(obj.gpindices, slice) else obj.gpindices

From dd679e8c00c2759819ce90b967cca8082544069d Mon Sep 17 00:00:00 2001
From: Juan G Mendoza <jgonz19@sandia.gov>
Date: Wed, 25 Jun 2025 17:37:27 -0600
Subject: [PATCH 08/34] started docstring

---
 pygsti/models/fogistore.py | 22 ++++++++++++++++++++
 pygsti/models/model.py     | 42 +++++++++++++++++++-------------------
 2 files changed, 43 insertions(+), 21 deletions(-)

diff --git a/pygsti/models/fogistore.py b/pygsti/models/fogistore.py
index ecdac056e..a9fe0ebc6 100644
--- a/pygsti/models/fogistore.py
+++ b/pygsti/models/fogistore.py
@@ -34,6 +34,27 @@ class FirstOrderGaugeInvariantStore(_NicelySerializable):
     Currently, it is only compatible with :class:`ExplicitOpModel` objects.
     """
     def __init__(self, primitive_op_labels, gauge_space, elem_errorgen_labels_by_op, op_errorgen_indices, fogi_directions, fogi_metadata, dependent_dir_indices, fogv_directions, allop_gauge_action, gauge_space_directions, norm_order=None, dependent_fogi_action='drop'):
+        """
+
+        Parameters
+        ----------
+        primitive_op_labels : tuple of Label
+            Labels describing the gate set operations
+
+        gauge_space : ErrorgenSpace
+            Special way of intersecting the gauge spaces of all ops, see the beginning of from_gauge_action_matrices 
+            for a more detailed description
+        elem_errorgen_labels_by_op (_type_): _description_
+        op_errorgen_indices (_type_): _description_
+        fogi_directions (_type_): _description_
+        fogi_metadata (_type_): _description_
+        dependent_dir_indices (_type_): _description_
+        fogv_directions (_type_): _description_
+        allop_gauge_action (_type_): _description_
+        gauge_space_directions (_type_): _description_
+        norm_order (_type_, optional): _description_. Defaults to None.
+        dependent_fogi_action (str, optional): _description_. Defaults to 'drop'.
+        """
         super().__init__()
         self.primitive_op_labels = primitive_op_labels
         self.gauge_space = gauge_space
@@ -62,6 +83,7 @@ def from_gauge_action_matrices(cls, gauge_action_matrices_by_op, gauge_action_ga
 
         primitive_op_labels = tuple(gauge_action_matrices_by_op.keys())
 
+        
         # Construct common gauge space by special way of intersecting the gauge spaces for all the ops
         # Note: the gauge_space of each op is constructed (see `setup_fogi`) so that the gauge action is
         #  zero on any elementary error generator not in the elementary-errorgen basis associated with
diff --git a/pygsti/models/model.py b/pygsti/models/model.py
index 7ff38dd7c..7bbe62d19 100644
--- a/pygsti/models/model.py
+++ b/pygsti/models/model.py
@@ -434,27 +434,6 @@ class OpModel(Model):
     possible to build up layer operations in an on-demand fashion from pieces
     within the model.
 
-    Parameters
-    ----------
-    state_space : StateSpace
-        The state space for this model.
-
-    basis : Basis
-        The basis used for the state space by dense operator representations.
-
-    evotype : Evotype or str, optional
-        The evolution type of this model, describing how states are
-        represented.  The special value `"default"` is equivalent
-        to specifying the value of `pygsti.evotypes.Evotype.default_evotype`.
-
-    layer_rules : LayerRules
-        The "layer rules" used for constructing operators for circuit
-        layers.  This functionality is essential to using this model to
-        simulate ciruits, and is typically supplied by derived classes.
-
-    simulator : ForwardSimulator or {"auto", "matrix", "map"}
-        The forward simulator (or typ) that this model should use.  `"auto"`
-        tries to determine the best type automatically.
     """
 
     #Whether to perform extra parameter-vector integrity checks
@@ -466,6 +445,27 @@ class OpModel(Model):
     def __init__(self, state_space, basis, evotype, layer_rules, simulator="auto"):
         """
         Creates a new OpModel.  Rarely used except from derived classes `__init__` functions.
+        Parameters
+        ----------
+        state_space : StateSpace
+            The state space for this model.
+
+        basis : Basis
+            The basis used for the state space by dense operator representations.
+
+        evotype : Evotype or str, optional
+            The evolution type of this model, describing how states are
+            represented.  The special value `"default"` is equivalent
+            to specifying the value of `pygsti.evotypes.Evotype.default_evotype`.
+
+        layer_rules : LayerRules
+            The "layer rules" used for constructing operators for circuit
+            layers.  This functionality is essential to using this model to
+            simulate ciruits, and is typically supplied by derived classes.
+
+        simulator : ForwardSimulator or {"auto", "matrix", "map"}
+            The forward simulator (or typ) that this model should use.  `"auto"`
+            tries to determine the best type automatically.
         """
         self._set_state_space(state_space, basis)
         #sets self._state_space, self._basis

From 8b2f67685b087886ebb1f2c0719acbf47b0219f8 Mon Sep 17 00:00:00 2001
From: Juan G Mendoza <jgonz19@sandia.gov>
Date: Thu, 26 Jun 2025 15:25:15 -0600
Subject: [PATCH 09/34] finished fogistore docstring

---
 pygsti/models/fogistore.py | 91 +++++++++++++++++++++++++++++---------
 pygsti/tools/fogitools.py  |  9 +++-
 2 files changed, 79 insertions(+), 21 deletions(-)

diff --git a/pygsti/models/fogistore.py b/pygsti/models/fogistore.py
index a9fe0ebc6..3da0415e8 100644
--- a/pygsti/models/fogistore.py
+++ b/pygsti/models/fogistore.py
@@ -33,7 +33,9 @@ class FirstOrderGaugeInvariantStore(_NicelySerializable):
 
     Currently, it is only compatible with :class:`ExplicitOpModel` objects.
     """
-    def __init__(self, primitive_op_labels, gauge_space, elem_errorgen_labels_by_op, op_errorgen_indices, fogi_directions, fogi_metadata, dependent_dir_indices, fogv_directions, allop_gauge_action, gauge_space_directions, norm_order=None, dependent_fogi_action='drop'):
+    def __init__(self, primitive_op_labels, gauge_space, elem_errorgen_labels_by_op, op_errorgen_indices, fogi_directions, 
+                 fogi_metadata, dependent_dir_indices, fogv_directions, allop_gauge_action, gauge_space_directions, 
+                 norm_order='auto', dependent_fogi_action='drop'):
         """
 
         Parameters
@@ -42,18 +44,70 @@ def __init__(self, primitive_op_labels, gauge_space, elem_errorgen_labels_by_op,
             Labels describing the gate set operations
 
         gauge_space : ErrorgenSpace
-            Special way of intersecting the gauge spaces of all ops, see the beginning of from_gauge_action_matrices 
-            for a more detailed description
-        elem_errorgen_labels_by_op (_type_): _description_
-        op_errorgen_indices (_type_): _description_
-        fogi_directions (_type_): _description_
-        fogi_metadata (_type_): _description_
-        dependent_dir_indices (_type_): _description_
-        fogv_directions (_type_): _description_
-        allop_gauge_action (_type_): _description_
-        gauge_space_directions (_type_): _description_
-        norm_order (_type_, optional): _description_. Defaults to None.
-        dependent_fogi_action (str, optional): _description_. Defaults to 'drop'.
+            Special way of intersecting the gauge spaces of all ops, see the beginning of from_gauge_action_matrices()
+            for a more detailed description.
+
+        elem_errorgen_labels_by_op : dict of Label keys with GlobalElementaryErrorgenLabel or GlobalElementaryErrorgenLabel values
+            Elementary error generator labels for every operation in the gate set
+
+        op_errorgen_indices : dict of Label keys with Slices values
+            Specifies which indices from the error generator space labels correspond to which operation. For example,
+            these allows us to index accordingly into the rows of allop_gauge_action which contain all error generators
+            from every operation stacked on top of each other.
+
+        fogi_directions : numpy array
+            Array where each column is a FOGI direction, the rows represent the modelmember parameters from which they were
+            derived (typically elementary error generators)
+
+        fogi_metadata : list of dictionaries 
+            Entry number k of this list contains metadata information about fogi quantity k (column k of fogi_directions). 
+            It is organized as a dictionary with the following entries:
+
+                - 'name' (str) : Label describing the FOGI quantity
+
+                - 'abbrev' (str) : An abbreviated version of the label above TODO: more details 
+
+                - 'r' (float) : constant that allows us to convert between the gauge-space-normalized
+                       FOGI direction to the errgen-space-normalized version of it.
+                       fogi_dir = r * dot(M, epsilon) = r * dot(inv_diff_gauge_action, int_vec)
+
+                - 'gaugespace_dir' (numpy array) : A corresponding gauge direction TODO: more details on their relationship
+                                    to FOGI quantities
+
+                - 'opset' (list of Label) : The set of model member operations that are affected by the corresponding FOGI
+                           quantity
+
+                - 'raw' (str) : Another version of 'name' without abbreviations
+
+        dependent_dir_indices : list of int
+            Specifies the indices of a set of FOGI quantities which, if removed, the remaining FOGI quantities are linearly independent.
+
+        fogv_directions : numpy array
+            A list of gauge directions that span the gauge space of the target gate set
+
+        allop_gauge_action : lil_matrix
+            A restricted version of each operation's gauge_action stacked together. Ignores elemgens that are not included in the 
+            common gauge space. The common gauge space is the intersection of all the individual gauge spaces. TODO: I believe 
+            this is done to remove the "trivial"/"meta" gauge from SPAM operations. We should verify this and if so, include
+            it in this docstring
+
+        gauge_space_directions : numpy array 
+            Contains, when applicable, a gauge-space direction that
+            correspondings to the FOGI quantity in fogi_directions.  Such a gauge-space direction
+            exists for relational FOGI quantities, where the FOGI quantity is constructed by taking the
+            *difference* of a gauge-space action (given by the gauge-space direction) on two operations
+            (or sets of operations). TODO: How does one know the corresponding fogi direction? each
+            one of these is computed from fogv_directions.
+
+        norm_order : str or int
+
+            Defines the order of the norm to normalize FOGI directions. It should be 1 for normalizing 'S' quantities and 2 for 'H',
+            so 'auto' utilizes intelligence.
+
+        dependent_fogi_action : 'drop' or 'mark'
+
+            If 'drop', all linearly dependent FOGI directions are not returned, resulting in a linearly independent set of quantities. 
+            If 'mark' linearly dependent FOGI directions are kept and marked by dependent_dir_indices.
         """
         super().__init__()
         self.primitive_op_labels = primitive_op_labels
@@ -76,7 +130,7 @@ def __init__(self, primitive_op_labels, gauge_space, elem_errorgen_labels_by_op,
     @classmethod
     def from_gauge_action_matrices(cls, gauge_action_matrices_by_op, gauge_action_gauge_spaces_by_op, errorgen_coefficient_labels_by_op,
                  op_label_abbrevs=None, reduce_to_model_space=True,
-                 dependent_fogi_action='drop', norm_order=None):
+                 dependent_fogi_action='drop', norm_order='auto'):
         """
         TODO: docstring
         """
@@ -176,23 +230,20 @@ def from_gauge_action_matrices(cls, gauge_action_matrices_by_op, gauge_action_ga
 
         #Get gauge-space directions corresponding to the fogv directions
         # (pinv_allop_gauge_action takes errorgen-set -> gauge-gen space)
-        allop_gauge_action = allop_gauge_action
         pinv_allop_gauge_action = _np.linalg.pinv(allop_gauge_action.toarray(), rcond=1e-7)
-        gauge_space_directions = _np.dot(pinv_allop_gauge_action,
-                                         fogv_directions.toarray())  # in gauge-generator space
-        gauge_space_directions = gauge_space_directions
+        gauge_space_directions = pinv_allop_gauge_action @ fogv_directions.toarray()  # in gauge-generator space
 
         #Notes on error-gen vs gauge-gen space:
         # self.fogi_directions and self.fogv_directions are dual vectors in error-generator space,
         # i.e. with elements corresponding to the elementary error generators given by
         # self.errorgen_space_op_elem_labels.
-        # self.fogi_gaugespace_directions contains, when applicable, a gauge-space direction that
+        # self.fogi_gauge_space_directions contains, when applicable, a gauge-space direction that
         # correspondings to the FOGI quantity in self.fogi_directions.  Such a gauge-space direction
         # exists for relational FOGI quantities, where the FOGI quantity is constructed by taking the
         # *difference* of a gauge-space action (given by the gauge-space direction) on two operations
         # (or sets of operations).
 
-        #Store auxiliary info for later use
+        
         fogi_store = cls(primitive_op_labels, gauge_space, elem_errorgen_labels_by_op, op_errorgen_indices, fogi_directions, fogi_metadata, dependent_dir_indices, fogv_directions, allop_gauge_action, gauge_space_directions, norm_order, dependent_fogi_action)
         fogi_store._check_fogi_store()
         return fogi_store
diff --git a/pygsti/tools/fogitools.py b/pygsti/tools/fogitools.py
index f241867ef..e51c3add6 100644
--- a/pygsti/tools/fogitools.py
+++ b/pygsti/tools/fogitools.py
@@ -339,7 +339,14 @@ def _create_op_errgen_indices_dict(primitive_op_labels, errorgen_coefficient_lab
 def construct_fogi_quantities(primitive_op_labels, gauge_action_matrices,
                               errorgen_coefficient_labels, op_errgen_indices, gauge_space,
                               op_label_abbrevs=None, dependent_fogi_action='drop', norm_order=None):
-    """ TODO: docstring """
+    """ TODO: docstring 
+    Parameters
+    ----------
+    dependent_fogi_action : 'drop' or 'mark'
+        'drop' ensures the set of FOGI quantities returned is linearly independent. 'mark' ensures 
+
+
+    """
     assert(dependent_fogi_action in ('drop', 'mark'))
     orthogonalize_relationals = True
 

From 2a46692eeab2979805258806f20d431fcc3f14ed Mon Sep 17 00:00:00 2001
From: Juan G Mendoza <jgonz19@sandia.gov>
Date: Fri, 27 Jun 2025 09:41:04 -0600
Subject: [PATCH 10/34] implemented and tested __eq__ method for fogistore

---
 pygsti/baseobjs/errorgenbasis.py      |  3 ++
 pygsti/baseobjs/errorgenspace.py      |  3 ++
 pygsti/models/fogistore.py            | 76 ++++++++++++++++++++++++++-
 pygsti/models/modelparaminterposer.py |  6 ++-
 4 files changed, 86 insertions(+), 2 deletions(-)

diff --git a/pygsti/baseobjs/errorgenbasis.py b/pygsti/baseobjs/errorgenbasis.py
index 3a39af444..6e8b1d974 100644
--- a/pygsti/baseobjs/errorgenbasis.py
+++ b/pygsti/baseobjs/errorgenbasis.py
@@ -119,6 +119,9 @@ def __init__(self, state_space, labels, basis_1q=None):
         self._cached_dual_matrices = None
         self._cached_supports = None
 
+    def __eq__(self, other):
+        return self.state_space.__eq__(other.state_space) and [label.__str__() for label in self.labels] == [label.__str__() for label in other.labels] and self._basis_1q.__eq__(other._basis_1q)
+    
     def _to_nice_serialization(self):
         state = super()._to_nice_serialization()
         state.update({'state_space' : self.state_space._to_nice_serialization(),
diff --git a/pygsti/baseobjs/errorgenspace.py b/pygsti/baseobjs/errorgenspace.py
index 367515bcc..a75070368 100644
--- a/pygsti/baseobjs/errorgenspace.py
+++ b/pygsti/baseobjs/errorgenspace.py
@@ -84,6 +84,9 @@ def intersection(self, other_space, free_on_unspecified_space=False, use_nice_nu
 
         return ErrorgenSpace(intersection_vecs, common_basis)
 
+    def __eq__(self, other):
+
+        return _np.allclose(self.vectors, other.vectors) and self.elemgen_basis.__eq__(other.elemgen_basis)
     def union(self, other_space):
         """
         TODO: docstring
diff --git a/pygsti/models/fogistore.py b/pygsti/models/fogistore.py
index 3da0415e8..6a5eeb295 100644
--- a/pygsti/models/fogistore.py
+++ b/pygsti/models/fogistore.py
@@ -86,7 +86,7 @@ def __init__(self, primitive_op_labels, gauge_space, elem_errorgen_labels_by_op,
             A list of gauge directions that span the gauge space of the target gate set
 
         allop_gauge_action : lil_matrix
-            A restricted version of each operation's gauge_action stacked together. Ignores elemgens that are not included in the 
+            Stands for "All operations gauge action". It is a restricted version of each operation's gauge_action stacked together. Ignores elemgens that are not included in the 
             common gauge space. The common gauge space is the intersection of all the individual gauge spaces. TODO: I believe 
             this is done to remove the "trivial"/"meta" gauge from SPAM operations. We should verify this and if so, include
             it in this docstring
@@ -359,7 +359,81 @@ def find_nice_fogiv_directions(self):
         #gauge_space_directions = _np.dot(pinv_allop_gauge_action, self.fogv_directions)  # in gauge-generator space
         #assert(_np.linalg.matrix_rank(gauge_space_directions) <= self._gauge_space_dim)  # should be nearly full rank
         pass
+    
+    def __eq__(self, other):
 
+        
+        if [label.__str__() for label in self.primitive_op_labels] != [label.__str__() for label in other.primitive_op_labels]:
+            return False
+        
+        if not self.gauge_space.__eq__(other.gauge_space):
+            return False
+        
+        if [op.__str__() for op in self.elem_errorgen_labels_by_op.keys()] != [op.__str__() for op in other.elem_errorgen_labels_by_op.keys()]:
+            return False
+        
+        for errgens_self, errgens_other in zip(list(self.elem_errorgen_labels_by_op.values()), list(other.elem_errorgen_labels_by_op.values())):
+            for errgen_self, errgen_other in zip(errgens_self, errgens_other):
+                if errgen_self.__str__() != errgen_other.__str__():
+                    return False
+        
+        if not (self.fogi_directions != other.fogi_directions).nnz == 0:
+            return False
+        
+        for object_self, object_other in zip(self.fogi_metadata, other.fogi_metadata):
+
+            if (object_self['name'] != object_other['name']) or (object_self['abbrev'] != object_other['abbrev']) or (object_self['r'] != object_other['r']) \
+                or (object_self['raw'] != object_other['raw']):
+                return False
+            #Case for when one of the spaces is None but the other isn't
+            if (object_self['gaugespace_dir'] is not None and object_other['gaugespace_dir'] is None) or \
+            (object_self['gaugespace_dir'] is None and object_other['gaugespace_dir'] is not None):
+                return False
+            #If they are both not none, check that their matrices are equal
+            if object_self['gaugespace_dir'] is not None and object_other['gaugespace_dir'] is not None:
+                if not _np.allclose(object_self['gaugespace_dir'], object_other['gaugespace_dir']):
+                    return False
+            if [label.__str__() for label in object_self['opset']] != [label.__str__() for label in object_other['opset']]:
+                return False
+        
+        if not _np.allclose(self.dependent_dir_indices, other.dependent_dir_indices):
+            return False
+        
+        if not (self.fogv_directions != other.fogv_directions).nnz == 0:
+            return False 
+        
+        if not (self.allop_gauge_action != other.allop_gauge_action).nnz == 0:
+            return False
+        #Case for when one of the spaces is none but the other isn't
+        if (self.gauge_space_directions is not None and other.gauge_space_directions is None) or \
+            (self.gauge_space_directions is None and other.gauge_space_directions is not None):
+            return False
+        #If they are both not none, check that their matrices are equal
+        if self.gauge_space_directions is not None and other.gauge_space_directions is not None:
+            if not _np.allclose(self.gauge_space_directions, other.gauge_space_directions):
+                return False
+        
+        if self.norm_order != other.norm_order or self._dependent_fogi_action != other._dependent_fogi_action:
+            return False
+        return True
+        
+        '''
+        self.primitive_op_labels = primitive_op_labels
+        self.gauge_space = gauge_space
+        self.elem_errorgen_labels_by_op = elem_errorgen_labels_by_op
+        self.op_errorgen_indices = op_errorgen_indices
+        self.fogi_directions = fogi_directions
+        self.fogi_metadata = fogi_metadata
+        self.dependent_dir_indices = dependent_dir_indices
+        self.fogv_directions = fogv_directions
+        self.allop_gauge_action = allop_gauge_action
+        self.gauge_space_directions = gauge_space_directions
+        self.norm_order = norm_order
+        self._dependent_fogi_action = dependent_fogi_action
+
+        self.errorgen_space_op_elem_labels = tuple([(op_label, elem_lbl) for op_label in self.primitive_op_labels
+                                                    for elem_lbl in self.elem_errorgen_labels_by_op[op_label]])
+        self.fogv_labels = ["%d gauge action" % i for i in range(self.fogv_directions.shape[1])]'''
     @property
     def errorgen_space_dim(self):
         return self.fogi_directions.shape[0]
diff --git a/pygsti/models/modelparaminterposer.py b/pygsti/models/modelparaminterposer.py
index 4b92629c9..1d71e5f5c 100644
--- a/pygsti/models/modelparaminterposer.py
+++ b/pygsti/models/modelparaminterposer.py
@@ -99,4 +99,8 @@ def _to_nice_serialization(self):
     @classmethod
     def _from_nice_serialization(cls, state):  # memo holds already de-serialized objects
         return cls(cls._decodemx(state['transform_matrix']))
-
+    
+    def __eq__(self, other):
+        assert isinstance(other, LinearInterposer)
+        return _np.allclose(self.transform_matrix, other.transform_matrix)
+    
\ No newline at end of file

From 6703ab8b9d12fa79db5d93e5b95cc2e85199265b Mon Sep 17 00:00:00 2001
From: Juan G Mendoza <jgonz19@sandia.gov>
Date: Fri, 27 Jun 2025 10:52:51 -0600
Subject: [PATCH 11/34] added test for FOGI equal method

---
 pygsti/models/fogistore.py            |  1 +
 pygsti/models/modelparaminterposer.py |  5 +++-
 test/unit/objects/test_fogi.py        | 35 +++++++++++++++++++++++++++
 3 files changed, 40 insertions(+), 1 deletion(-)

diff --git a/pygsti/models/fogistore.py b/pygsti/models/fogistore.py
index 6a5eeb295..ab08b0c24 100644
--- a/pygsti/models/fogistore.py
+++ b/pygsti/models/fogistore.py
@@ -362,6 +362,7 @@ def find_nice_fogiv_directions(self):
     
     def __eq__(self, other):
 
+        assert isinstance(other, FirstOrderGaugeInvariantStore), 'Object provided is not of type FirstOrderGaugeInvariantStore'
         
         if [label.__str__() for label in self.primitive_op_labels] != [label.__str__() for label in other.primitive_op_labels]:
             return False
diff --git a/pygsti/models/modelparaminterposer.py b/pygsti/models/modelparaminterposer.py
index 1d71e5f5c..d6e53ba3c 100644
--- a/pygsti/models/modelparaminterposer.py
+++ b/pygsti/models/modelparaminterposer.py
@@ -101,6 +101,9 @@ def _from_nice_serialization(cls, state):  # memo holds already de-serialized ob
         return cls(cls._decodemx(state['transform_matrix']))
     
     def __eq__(self, other):
-        assert isinstance(other, LinearInterposer)
+        assert isinstance(other, LinearInterposer), 'Object provided is not of LinearInterposer type'
+
+        if self.transform_matrix.shape != other.transform_matrix.shape:
+            return False
         return _np.allclose(self.transform_matrix, other.transform_matrix)
     
\ No newline at end of file
diff --git a/test/unit/objects/test_fogi.py b/test/unit/objects/test_fogi.py
index d55314aa2..67cd14d44 100644
--- a/test/unit/objects/test_fogi.py
+++ b/test/unit/objects/test_fogi.py
@@ -5,6 +5,8 @@
 
 from ..util import BaseCase
 from pygsti.modelpacks import smq1Q_XYI as std
+from pygsti.modelpacks import smq1Q_XY as std2
+from pygsti.modelpacks import smq1Q_XZ as std3
 from pygsti.baseobjs import Basis, CompleteElementaryErrorgenBasis
 from pygsti.processors import QubitProcessorSpec
 from pygsti.models import create_crosstalk_free_model
@@ -228,3 +230,36 @@ def test_cloud_crosstalk_fogi(self):
         w = np.random.rand(mdl.num_params)
         w[0:nprefix] = 0 # zero out all unused params (these can be SPAM and can't be any value?)
         mdl.from_vector(w)
+
+    def test_equal_method(self):
+
+        def equal_fogi_models(fogi_model, fogi_model2):
+
+            return fogi_model.fogi_store.__eq__(fogi_model2.fogi_store) and fogi_model.param_interposer.__eq__(fogi_model2.param_interposer)
+        model = std.target_model('GLND')
+        model2 = std2.target_model('GLND')
+        model3 = std3.target_model('GLND')
+
+        basis1q = Basis.cast('pp', 4)
+        gauge_basis = CompleteElementaryErrorgenBasis(
+            basis1q, model.state_space, elementary_errorgen_types='HSCA')
+        gauge_basis2 = CompleteElementaryErrorgenBasis(
+            basis1q, model2.state_space, elementary_errorgen_types='HSCA')
+        gauge_basis3 = CompleteElementaryErrorgenBasis(
+            basis1q, model3.state_space, elementary_errorgen_types='HSCA')
+        
+        model.setup_fogi(gauge_basis, None, None, reparameterize=True, dependent_fogi_action='drop', include_spam=True)
+        model2.setup_fogi(gauge_basis2, None, None, reparameterize=True, dependent_fogi_action='drop', include_spam=True)
+        model3.setup_fogi(gauge_basis3, None, None, reparameterize=True, dependent_fogi_action='drop', include_spam=True)
+        
+        msg = 'FOGI models that are the same are identified as different by __eq__ methods'
+        self.assertTrue(equal_fogi_models(model, model), msg=msg)
+        self.assertTrue(equal_fogi_models(model2, model2), msg=msg)
+        self.assertTrue(equal_fogi_models(model3, model3), msg=msg)
+
+        msg = 'FOGI models that are different are not recognized as different by __eq__ methods'
+        self.assertFalse(equal_fogi_models(model, model2), msg=msg)
+        self.assertFalse(equal_fogi_models(model, model3), msg=msg)
+        self.assertFalse(equal_fogi_models(model2, model3), msg=msg)
+
+

From bc39b125e08e406e1705ecaf78a012a112c1af88 Mon Sep 17 00:00:00 2001
From: Juan G Mendoza <jgonz19@sandia.gov>
Date: Fri, 27 Jun 2025 11:24:27 -0600
Subject: [PATCH 12/34] added tests for fogi serialization

---
 test/unit/io/test_nice_serialization.py |  2 ++
 test/unit/objects/test_fogi.py          | 32 +++++++++++++++++++++----
 2 files changed, 30 insertions(+), 4 deletions(-)

diff --git a/test/unit/io/test_nice_serialization.py b/test/unit/io/test_nice_serialization.py
index 3d656a2ce..53a2c3437 100644
--- a/test/unit/io/test_nice_serialization.py
+++ b/test/unit/io/test_nice_serialization.py
@@ -69,6 +69,8 @@ def test_cloudnoise_model(self, pth):
         self.assertTrue(mdl_cloud.is_equivalent(mdl_cloud2))
 
 
+#Somewhat random question for other pyGSTi developers: 
+#Am I missing something or does this class not test anything related to serialization?
 class ModelEquivalenceTester(BaseCase):
 
     def setUp(self):
diff --git a/test/unit/objects/test_fogi.py b/test/unit/objects/test_fogi.py
index 67cd14d44..e05326bc7 100644
--- a/test/unit/objects/test_fogi.py
+++ b/test/unit/objects/test_fogi.py
@@ -2,15 +2,14 @@
 
 import sys
 import numpy as np
-
-from ..util import BaseCase
+from ..util import BaseCase, with_temp_path
 from pygsti.modelpacks import smq1Q_XYI as std
 from pygsti.modelpacks import smq1Q_XY as std2
 from pygsti.modelpacks import smq1Q_XZ as std3
 from pygsti.baseobjs import Basis, CompleteElementaryErrorgenBasis
 from pygsti.processors import QubitProcessorSpec
 from pygsti.models import create_crosstalk_free_model
-from pygsti.models import create_cloud_crosstalk_model_from_hops_and_weights
+from pygsti.models import create_cloud_crosstalk_model_from_hops_and_weights, Model
 
 
 class FogiTester(BaseCase):
@@ -234,8 +233,8 @@ def test_cloud_crosstalk_fogi(self):
     def test_equal_method(self):
 
         def equal_fogi_models(fogi_model, fogi_model2):
-
             return fogi_model.fogi_store.__eq__(fogi_model2.fogi_store) and fogi_model.param_interposer.__eq__(fogi_model2.param_interposer)
+        
         model = std.target_model('GLND')
         model2 = std2.target_model('GLND')
         model3 = std3.target_model('GLND')
@@ -262,4 +261,29 @@ def equal_fogi_models(fogi_model, fogi_model2):
         self.assertFalse(equal_fogi_models(model, model3), msg=msg)
         self.assertFalse(equal_fogi_models(model2, model3), msg=msg)
 
+    #TODO: should this be in test_nice_serialization instead?
+    @with_temp_path
+    def test_fogi_serialization(self, temp_pth):
+        def equal_fogi_models(fogi_model, fogi_model2):
+            return fogi_model.fogi_store.__eq__(fogi_model2.fogi_store) and fogi_model.param_interposer.__eq__(fogi_model2.param_interposer)
+        
+        model = std.target_model('GLND')
+        model2 = std2.target_model('GLND')
+        model3 = std3.target_model('GLND')
+
+        model.write(temp_pth + '.json')
+        loaded_model = Model.read(temp_pth + '.json')
+        model2.write(temp_pth + '.json')
+        loaded_model2 = Model.read(temp_pth + '.json')
+        model3.write(temp_pth + '.json')
+        loaded_model3 = Model.read(temp_pth + '.json')
+
+        self.assertTrue(equal_fogi_models(model, loaded_model))
+        self.assertTrue(equal_fogi_models(model2, loaded_model2))
+        self.assertTrue(equal_fogi_models(model3, loaded_model3))
+        
+
+        
+
+
 

From f40ee352fa8928e95d17025a79b0882520d2e2d8 Mon Sep 17 00:00:00 2001
From: Juan G Mendoza <jgonz19@sandia.gov>
Date: Fri, 27 Jun 2025 13:00:39 -0600
Subject: [PATCH 13/34] added set_param_values unit tests for
 param_interposer/fogi models

---
 test/unit/io/test_nice_serialization.py |  4 +---
 test/unit/objects/test_fogi.py          | 31 ++++++++++++++++++++++++-
 test/unit/objects/test_model.py         | 17 ++++++++++++--
 3 files changed, 46 insertions(+), 6 deletions(-)

diff --git a/test/unit/io/test_nice_serialization.py b/test/unit/io/test_nice_serialization.py
index 53a2c3437..c71c25b41 100644
--- a/test/unit/io/test_nice_serialization.py
+++ b/test/unit/io/test_nice_serialization.py
@@ -121,6 +121,4 @@ def test_cloud_model_equal(self):
                                                  stochastic_error_probs={'Gy': (0.01, 0.02, 0.03)},
                                                  lindblad_error_coeffs={'Gcnot': {('H','ZZ'): 0.07, ('S','XX'): 0.10}},
                                                  independent_gates=True, independent_spam=True, verbosity=2)
-        self.check_model(mdl_cloud)
-
-
+        self.check_model(mdl_cloud)
\ No newline at end of file
diff --git a/test/unit/objects/test_fogi.py b/test/unit/objects/test_fogi.py
index e05326bc7..997cac98e 100644
--- a/test/unit/objects/test_fogi.py
+++ b/test/unit/objects/test_fogi.py
@@ -271,6 +271,18 @@ def equal_fogi_models(fogi_model, fogi_model2):
         model2 = std2.target_model('GLND')
         model3 = std3.target_model('GLND')
 
+        basis1q = Basis.cast('pp', 4)
+        gauge_basis = CompleteElementaryErrorgenBasis(
+            basis1q, model.state_space, elementary_errorgen_types='HSCA')
+        gauge_basis2 = CompleteElementaryErrorgenBasis(
+            basis1q, model2.state_space, elementary_errorgen_types='HSCA')
+        gauge_basis3 = CompleteElementaryErrorgenBasis(
+            basis1q, model3.state_space, elementary_errorgen_types='HSCA')
+        
+        model.setup_fogi(gauge_basis, None, None, reparameterize=True, dependent_fogi_action='drop', include_spam=True)
+        model2.setup_fogi(gauge_basis2, None, None, reparameterize=True, dependent_fogi_action='drop', include_spam=True)
+        model3.setup_fogi(gauge_basis3, None, None, reparameterize=True, dependent_fogi_action='drop', include_spam=True)
+
         model.write(temp_pth + '.json')
         loaded_model = Model.read(temp_pth + '.json')
         model2.write(temp_pth + '.json')
@@ -282,8 +294,25 @@ def equal_fogi_models(fogi_model, fogi_model2):
         self.assertTrue(equal_fogi_models(model2, loaded_model2))
         self.assertTrue(equal_fogi_models(model3, loaded_model3))
         
-
+    def test_set_param_values(self):
         
+        model = std.target_model('GLND')
+        cp = model.copy()
+
+        basis1q = Basis.cast('pp', 4)
+        gauge_basis = CompleteElementaryErrorgenBasis(
+            basis1q, model.state_space, elementary_errorgen_types='HSCA')
 
+        model.setup_fogi(gauge_basis, None, None, reparameterize=True, dependent_fogi_action='drop', include_spam=True)
+        cp.setup_fogi(gauge_basis, None, None, reparameterize=True, dependent_fogi_action='drop', include_spam=True)
+
+        test_vec = np.arange(model.num_params) * 1e-3
+        cp.set_parameter_values(np.arange(model.num_params), test_vec)
+        model.from_vector(test_vec)
+        self.assertAlmostEqual(np.linalg.norm(model.to_vector() - cp.to_vector()), 0)
+        #TODO: Uncomment when issue #600 is resolved, remove line above
+        #self.assertTrue(cp.is_equivalent(cp2))
+        
+        
 
 
diff --git a/test/unit/objects/test_model.py b/test/unit/objects/test_model.py
index 8f5f0e0c9..66835a1d3 100644
--- a/test/unit/objects/test_model.py
+++ b/test/unit/objects/test_model.py
@@ -16,7 +16,6 @@
 from pygsti.circuits import Circuit
 from pygsti.models.gaugegroup import FullGaugeGroupElement
 from ..util import BaseCase, needs_cvxpy
-
 SKIP_DIAMONDIST_ON_WIN = True
 
 
@@ -191,7 +190,9 @@ def test_strdiff(self):
 
     def test_copy(self):
         gs2 = self.model.copy()
-        # TODO assert correctness
+        #TODO: Uncomment when issue #600 is resolved
+        #self.assertTrue(gs2.is_equivalent(self.model))
+        
 
     def test_deriv_wrt_params(self):
         deriv = self.model.deriv_wrt_params()
@@ -221,6 +222,18 @@ def test_vectorize(self):
         cp.from_vector(v)
         self.assertAlmostEqual(self.model.frobeniusdist(cp), 0)
 
+    def test_set_parameter_values(self):
+        if self.model.num_params > 0:
+            cp = self.model.copy()
+            cp2 = self.model.copy()
+            test_vec = np.arange(self.model.num_params) * 1e-3
+            cp.set_parameter_values(np.arange(self.model.num_params), test_vec)
+            cp2.from_vector(test_vec)
+            self.assertAlmostEqual(np.linalg.norm(cp2.to_vector() - cp.to_vector()), 0)
+            #TODO: Uncomment when issue #600 is resolved, remove line above
+            #self.assertTrue(cp.is_equivalent(cp2))
+
+
     def test_pickle(self):
         # XXX what exactly does this cover and is it needed?  EGN: this tests that the individual pieces (~dicts) within a model can be pickled; it's useful for debuggin b/c often just one of these will break.
         p = pickle.dumps(self.model.preps)

From acf96f1c38c03b5888c594a66a1f342c3fe38dfd Mon Sep 17 00:00:00 2001
From: Juan G Mendoza <jgonz19@sandia.gov>
Date: Mon, 30 Jun 2025 11:57:29 -0600
Subject: [PATCH 14/34] added a docstring and fixed vector decoding for PR

---
 pygsti/baseobjs/errorgenspace.py | 10 +++++++++-
 1 file changed, 9 insertions(+), 1 deletion(-)

diff --git a/pygsti/baseobjs/errorgenspace.py b/pygsti/baseobjs/errorgenspace.py
index a75070368..f378278a2 100644
--- a/pygsti/baseobjs/errorgenspace.py
+++ b/pygsti/baseobjs/errorgenspace.py
@@ -21,6 +21,14 @@ class ErrorgenSpace(_NicelySerializable):
 
     This object collects the information needed to specify a space
     within the space of all error generators.
+
+    Parameters
+    ----------
+    vectors : numpy array
+        List of vectors that span the space
+    
+        elemgen_basis : ElementaryErrorgenBasis
+            The elementary error generator basis that define the entries of self.vectors 
     """
 
     def __init__(self, vectors, basis):
@@ -39,7 +47,7 @@ def _to_nice_serialization(self):
         return state
     @classmethod
     def from_nice_serialization(cls, state):
-        return cls(state['vectors'], ExplicitElementaryErrorgenBasis.from_nice_serialization(state['basis']))
+        return cls(cls._decodemx(state['vectors']), ExplicitElementaryErrorgenBasis.from_nice_serialization(state['basis']))
     def intersection(self, other_space, free_on_unspecified_space=False, use_nice_nullspace=False):
         """
         TODO: docstring

From 17dab7528ef686cf275bcb75d42b5b350defdec8 Mon Sep 17 00:00:00 2001
From: Juan G Mendoza <jgonz19@sandia.gov>
Date: Mon, 30 Jun 2025 11:58:23 -0600
Subject: [PATCH 15/34] fixed indentation

---
 pygsti/baseobjs/errorgenspace.py | 27 ++++++++++++++-------------
 1 file changed, 14 insertions(+), 13 deletions(-)

diff --git a/pygsti/baseobjs/errorgenspace.py b/pygsti/baseobjs/errorgenspace.py
index f378278a2..de949a221 100644
--- a/pygsti/baseobjs/errorgenspace.py
+++ b/pygsti/baseobjs/errorgenspace.py
@@ -16,22 +16,23 @@
 from pygsti.baseobjs.errorgenbasis import ExplicitElementaryErrorgenBasis
 
 class ErrorgenSpace(_NicelySerializable):
-    """
-    A vector space of error generators, spanned by some basis.
 
-    This object collects the information needed to specify a space
-    within the space of all error generators.
-
-    Parameters
-    ----------
-    vectors : numpy array
-        List of vectors that span the space
-    
-        elemgen_basis : ElementaryErrorgenBasis
-            The elementary error generator basis that define the entries of self.vectors 
-    """
 
     def __init__(self, vectors, basis):
+        """
+        A vector space of error generators, spanned by some basis.
+
+        This object collects the information needed to specify a space
+        within the space of all error generators.
+
+        Parameters
+        ----------
+        vectors : numpy array
+            List of vectors that span the space
+        
+            elemgen_basis : ElementaryErrorgenBasis
+                The elementary error generator basis that define the entries of self.vectors 
+        """
         super().__init__()
         self.vectors = vectors
         self.elemgen_basis = basis

From 3d5b33da71482e37e144069165749057905ea24d Mon Sep 17 00:00:00 2001
From: Juan G Mendoza <jgonz19@sandia.gov>
Date: Mon, 30 Jun 2025 15:25:37 -0600
Subject: [PATCH 16/34] finished docstring for FOGIStore

---
 pygsti/models/fogistore.py | 42 ++++++++++++++++++++++++++++++++++----
 1 file changed, 38 insertions(+), 4 deletions(-)

diff --git a/pygsti/models/fogistore.py b/pygsti/models/fogistore.py
index ab08b0c24..5483f95f8 100644
--- a/pygsti/models/fogistore.py
+++ b/pygsti/models/fogistore.py
@@ -129,10 +129,44 @@ def __init__(self, primitive_op_labels, gauge_space, elem_errorgen_labels_by_op,
     
     @classmethod
     def from_gauge_action_matrices(cls, gauge_action_matrices_by_op, gauge_action_gauge_spaces_by_op, errorgen_coefficient_labels_by_op,
-                 op_label_abbrevs=None, reduce_to_model_space=True,
-                 dependent_fogi_action='drop', norm_order='auto'):
+                 op_label_abbrevs=None, dependent_fogi_action='drop', norm_order='auto'):
         """
-        TODO: docstring
+        Receives the matrices describing the individual gauge action of each operation within the gate set, and massages
+        them into appropriate from to be fed into construct_fogi_quantities() which finds FOGI quantities for the corresponding
+        gate set. Aditionally, gauge space vectors are computed, and gauge directions that "correspond" to FOGI quantities are 
+        identified. With all of this, a FOGIStore object is created and returned.
+
+        Parameters
+        ----------
+        gauge_action_matrices_by_op : dict of Label keys and numpy array values
+            Gauge action of the primitive operator described by Label. For SPAM, 
+            "gauge action" matrix is just the identity on the *relevant* space of gauge transformations.
+            For gates, it is K - UKU^dag where K is a gauge generator, and U is the ideal gate.
+
+        gauge_action_gauge_spaces_by_op : dict of Label keys and numpy array values
+            "Formatted" gauge spaces corresponding to each gauge action from gauge_action_matrices_by_op. 
+            This formatting occurs in model.py:_format_gauge_action_matrix() TODO: More details
+
+        errorgen_coefficient_labels_by_op : dict of Label keys and list of Label values
+            If indexed by an operation Label, this variable returns a list of error generator labels
+            that identify every row of their corresponding gauge action matrix in gauge_action_matrices_by_op
+
+        op_label_abbrevs : dict of Label keys with str values
+            The user is able to fill the values of this dictionary with any strings to represent
+            the operations in the gate set. typically used to choose shorter names to avoid
+            FOGI names that are too long/complicated to read. 
+
+
+        norm_order : int or 'auto' (Defaults to 'auto')
+
+            Defines the order of the norm to normalize FOGI directions. It should be 1 for 
+            normalizing 'S' quantities and 2 for 'H', so 'auto' utilizes intelligence.
+
+        dependent_fogi_action : 'drop' or 'mark' (Defaults to 'drop')
+
+            If 'drop', all linearly dependent FOGI directions are not returned, resulting in a linearly 
+            independent set of quantities. If 'mark' linearly dependent FOGI directions are kept and
+            marked by dependent_dir_indices.
         """
 
         primitive_op_labels = tuple(gauge_action_matrices_by_op.keys())
@@ -363,7 +397,7 @@ def find_nice_fogiv_directions(self):
     def __eq__(self, other):
 
         assert isinstance(other, FirstOrderGaugeInvariantStore), 'Object provided is not of type FirstOrderGaugeInvariantStore'
-        
+
         if [label.__str__() for label in self.primitive_op_labels] != [label.__str__() for label in other.primitive_op_labels]:
             return False
         

From 6ffdf9e7eaaa8b9b58afa12910ba9b544d9f1677 Mon Sep 17 00:00:00 2001
From: Juan G Mendoza <jgonz19@sandia.gov>
Date: Mon, 30 Jun 2025 16:17:31 -0600
Subject: [PATCH 17/34] fixed issue with FOGIStore method calling

---
 pygsti/models/model.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pygsti/models/model.py b/pygsti/models/model.py
index 7bbe62d19..f8e83d4e9 100644
--- a/pygsti/models/model.py
+++ b/pygsti/models/model.py
@@ -2671,7 +2671,7 @@ def create_complete_basis_fn(target_sslbls):
 
         self.fogi_store = _FOGIStore.from_gauge_action_matrices(gauge_action_matrices, gauge_action_gauge_spaces,
                                      errorgen_coefficient_labels,  # gauge_errgen_space_labels,
-                                     op_label_abbrevs, reduce_to_model_space, dependent_fogi_action,
+                                     op_label_abbrevs, dependent_fogi_action,
                                      norm_order=norm_order)
         if reparameterize:
             self.param_interposer = self._add_reparameterization(

From b6fde07f0212d360664778da6068e8910f28cb16 Mon Sep 17 00:00:00 2001
From: Juan G Mendoza <jgonz19@sandia.gov>
Date: Mon, 30 Jun 2025 16:18:15 -0600
Subject: [PATCH 18/34] fixed fogistore

---
 pygsti/models/model.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pygsti/models/model.py b/pygsti/models/model.py
index f8e83d4e9..7a92a7f45 100644
--- a/pygsti/models/model.py
+++ b/pygsti/models/model.py
@@ -2670,7 +2670,7 @@ def create_complete_basis_fn(target_sslbls):
 
 
         self.fogi_store = _FOGIStore.from_gauge_action_matrices(gauge_action_matrices, gauge_action_gauge_spaces,
-                                     errorgen_coefficient_labels,  # gauge_errgen_space_labels,
+                                     errorgen_coefficient_labels, 
                                      op_label_abbrevs, dependent_fogi_action,
                                      norm_order=norm_order)
         if reparameterize:

From 4f2ce405fa9d9518fe43df79b1ad8197dc4ec3ce Mon Sep 17 00:00:00 2001
From: Juan G Mendoza <jgonz19@sandia.gov>
Date: Wed, 2 Jul 2025 20:50:38 -0600
Subject: [PATCH 19/34] small patch on label_index problem that has been
 identified before

---
 pygsti/baseobjs/errorgenbasis.py | 12 ++++++++++--
 1 file changed, 10 insertions(+), 2 deletions(-)

diff --git a/pygsti/baseobjs/errorgenbasis.py b/pygsti/baseobjs/errorgenbasis.py
index 6e8b1d974..62282c028 100644
--- a/pygsti/baseobjs/errorgenbasis.py
+++ b/pygsti/baseobjs/errorgenbasis.py
@@ -780,7 +780,15 @@ def label_index(self, label, ok_if_missing=False, identity_label='I'):
         identity_label : str, optional (default 'I')
             An optional string specifying the label used to denote the identity in basis element labels.
         """
-        if isinstance(label, _LocalElementaryErrorgenLabel):
+        try:
+            return self.labels.index(label)
+        except ValueError as error:
+            
+            if ok_if_missing:
+                return None
+            else:
+                raise error
+        '''if isinstance(label, _LocalElementaryErrorgenLabel):
             label = _GlobalElementaryErrorgenLabel.cast(label, self.sslbls, identity_label=identity_label)
 
         support = label.sslbls
@@ -812,7 +820,7 @@ def label_index(self, label, ok_if_missing=False, identity_label='I'):
         else:
             raise ValueError("Invalid elementary errorgen type: %s" % str(eetype))
 
-        return base + indices[label]
+        return base + indices[label]'''
 
     def create_subbasis(self, sslbl_overlap, retain_max_weights=True):
         """

From 2f00ce8d9f7de2ff0bddefb64a22a564ea598682 Mon Sep 17 00:00:00 2001
From: Juan G Mendoza <jgonz19@sandia.gov>
Date: Thu, 3 Jul 2025 10:37:45 -0600
Subject: [PATCH 20/34] created amalgamation of Erik's label_index and mine, it
 uses his checks but does not return his result

---
 pygsti/baseobjs/errorgenbasis.py        | 22 ++++++++++++----------
 test/unit/objects/test_errorgenbasis.py |  1 +
 2 files changed, 13 insertions(+), 10 deletions(-)

diff --git a/pygsti/baseobjs/errorgenbasis.py b/pygsti/baseobjs/errorgenbasis.py
index 62282c028..923d60a61 100644
--- a/pygsti/baseobjs/errorgenbasis.py
+++ b/pygsti/baseobjs/errorgenbasis.py
@@ -780,17 +780,12 @@ def label_index(self, label, ok_if_missing=False, identity_label='I'):
         identity_label : str, optional (default 'I')
             An optional string specifying the label used to denote the identity in basis element labels.
         """
-        try:
-            return self.labels.index(label)
-        except ValueError as error:
-            
-            if ok_if_missing:
-                return None
-            else:
-                raise error
-        '''if isinstance(label, _LocalElementaryErrorgenLabel):
+        
+        if isinstance(label, _LocalElementaryErrorgenLabel):
             label = _GlobalElementaryErrorgenLabel.cast(label, self.sslbls, identity_label=identity_label)
 
+        
+
         support = label.sslbls
         eetype = label.errorgen_type
         bels = label.basis_element_labels
@@ -820,7 +815,14 @@ def label_index(self, label, ok_if_missing=False, identity_label='I'):
         else:
             raise ValueError("Invalid elementary errorgen type: %s" % str(eetype))
 
-        return base + indices[label]'''
+        try:
+            return self.labels.index(label)
+        except ValueError as error:
+            
+            if ok_if_missing:
+                return None
+            else:
+                raise error
 
     def create_subbasis(self, sslbl_overlap, retain_max_weights=True):
         """
diff --git a/test/unit/objects/test_errorgenbasis.py b/test/unit/objects/test_errorgenbasis.py
index 4d2f1da7f..5602641cf 100644
--- a/test/unit/objects/test_errorgenbasis.py
+++ b/test/unit/objects/test_errorgenbasis.py
@@ -112,6 +112,7 @@ def test_label_index(self):
 
         lbl_idx = self.complete_errorgen_basis_default_1Q.label_index(test_eg)
         lbl_idx_1 = self.complete_errorgen_basis_default_1Q.label_index(test_eg_local)
+        print(lbl_idx, lbl_idx_1, labels.index(test_eg))
         assert lbl_idx == lbl_idx_1
         assert lbl_idx ==  labels.index(test_eg)
 

From ad2526c715362514f3d1bfc3d00eee8c4d99c6b1 Mon Sep 17 00:00:00 2001
From: Juan G Mendoza <jgonz19@sandia.gov>
Date: Thu, 3 Jul 2025 10:56:38 -0600
Subject: [PATCH 21/34] removed debug print

---
 test/unit/objects/test_errorgenbasis.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/test/unit/objects/test_errorgenbasis.py b/test/unit/objects/test_errorgenbasis.py
index 5602641cf..4d2f1da7f 100644
--- a/test/unit/objects/test_errorgenbasis.py
+++ b/test/unit/objects/test_errorgenbasis.py
@@ -112,7 +112,6 @@ def test_label_index(self):
 
         lbl_idx = self.complete_errorgen_basis_default_1Q.label_index(test_eg)
         lbl_idx_1 = self.complete_errorgen_basis_default_1Q.label_index(test_eg_local)
-        print(lbl_idx, lbl_idx_1, labels.index(test_eg))
         assert lbl_idx == lbl_idx_1
         assert lbl_idx ==  labels.index(test_eg)
 

From ca32580e6983d991e321417597b3a4bbd2555666 Mon Sep 17 00:00:00 2001
From: Juan G Mendoza <jgonz19@sandia.gov>
Date: Thu, 3 Jul 2025 16:51:28 -0600
Subject: [PATCH 22/34] finally properly fixed FOGI label_index bug

---
 pygsti/baseobjs/errorgenbasis.py | 18 +++++-------------
 1 file changed, 5 insertions(+), 13 deletions(-)

diff --git a/pygsti/baseobjs/errorgenbasis.py b/pygsti/baseobjs/errorgenbasis.py
index 923d60a61..276316bbf 100644
--- a/pygsti/baseobjs/errorgenbasis.py
+++ b/pygsti/baseobjs/errorgenbasis.py
@@ -784,8 +784,6 @@ def label_index(self, label, ok_if_missing=False, identity_label='I'):
         if isinstance(label, _LocalElementaryErrorgenLabel):
             label = _GlobalElementaryErrorgenLabel.cast(label, self.sslbls, identity_label=identity_label)
 
-        
-
         support = label.sslbls
         eetype = label.errorgen_type
         bels = label.basis_element_labels
@@ -804,8 +802,10 @@ def label_index(self, label, ok_if_missing=False, identity_label='I'):
         elif eetype in ('C', 'A'):
             assert(len(trivial_bel) == 1)  # assumes this is a single character
             nontrivial_inds = [i for i, letter in enumerate(bels[0]) if letter != trivial_bel]
-            left_support = tuple([self.sslbls[i] for i in nontrivial_inds])
-
+            left_support = tuple([label.sslbls[i] for i in nontrivial_inds])
+            #left_support = tuple([self.sslbls[i] for i in nontrivial_inds])
+            #This line above is supposed to be a bugfix, I want to verify this with
+            #Corey before removing this comment
             if ok_if_missing and (support, left_support) not in self._offsets[eetype]:
                 return None
             base = self._offsets[eetype][(support, left_support)]
@@ -814,15 +814,7 @@ def label_index(self, label, ok_if_missing=False, identity_label='I'):
                 support, left_support, eetype, [trivial_bel], nontrivial_bels))}
         else:
             raise ValueError("Invalid elementary errorgen type: %s" % str(eetype))
-
-        try:
-            return self.labels.index(label)
-        except ValueError as error:
-            
-            if ok_if_missing:
-                return None
-            else:
-                raise error
+        return base + indices[label]
 
     def create_subbasis(self, sslbl_overlap, retain_max_weights=True):
         """

From f5aa81e42c5fa7ccc84679f59b5c26da1c65fdc1 Mon Sep 17 00:00:00 2001
From: Juan G Mendoza <jgonz19@sandia.gov>
Date: Tue, 8 Jul 2025 16:32:21 -0600
Subject: [PATCH 23/34] moved unit tests to correct file

---
 test/unit/io/test_nice_serialization.py | 53 -----------------------
 test/unit/objects/test_model.py         | 57 ++++++++++++++++++++++++-
 2 files changed, 56 insertions(+), 54 deletions(-)

diff --git a/test/unit/io/test_nice_serialization.py b/test/unit/io/test_nice_serialization.py
index c71c25b41..0500dcf88 100644
--- a/test/unit/io/test_nice_serialization.py
+++ b/test/unit/io/test_nice_serialization.py
@@ -69,56 +69,3 @@ def test_cloudnoise_model(self, pth):
         self.assertTrue(mdl_cloud.is_equivalent(mdl_cloud2))
 
 
-#Somewhat random question for other pyGSTi developers: 
-#Am I missing something or does this class not test anything related to serialization?
-class ModelEquivalenceTester(BaseCase):
-
-    def setUp(self):
-        nQubits = 2
-        self.pspec_2Q = QubitProcessorSpec(nQubits, ('Gx', 'Gy', 'Gcnot'), geometry="line",
-                                           qubit_labels=['qb{}'.format(i) for i in range(nQubits)])
-
-    def check_model(self, mdl):
-        mcopy = mdl.copy()
-        self.assertFalse(mdl is mcopy)
-        self.assertTrue(mcopy.is_similar(mdl))
-        self.assertTrue(mcopy.is_equivalent(mdl))
-
-        if mdl.num_params > 0:
-            r = np.random.random(mdl.num_params)
-            if np.linalg.norm(r) < 1e6:  # just in case we randomly get all zeros!
-                r = 0.1 * np.ones(mcopy.num_params)
-            v_prime = mdl.to_vector() + r
-            mcopy.from_vector(v_prime)
-            self.assertTrue(mcopy.is_similar(mdl))
-            self.assertFalse(mcopy.is_equivalent(mdl))
-
-    def test_explicit_model_equal(self):
-        mdl_explicit = smq1Q_XYI.target_model()
-        self.check_model(mdl_explicit)
-
-    def test_local_model_equal(self):
-        mdl_local = create_crosstalk_free_model(self.pspec_2Q,
-                                                ideal_gate_type='H+S', ideal_spam_type='tensor product H+S',
-                                                independent_gates=False,
-                                                ensure_composed_gates=False)
-        self.check_model(mdl_local)
-
-        mdl_local = create_crosstalk_free_model(self.pspec_2Q,
-                                                ideal_gate_type='H+S', ideal_spam_type='computational',
-                                                independent_gates=True,
-                                                ensure_composed_gates=True)
-        self.check_model(mdl_local)
-
-    def test_cloud_model_equal(self):
-        mdl_cloud = create_cloud_crosstalk_model(self.pspec_2Q, depolarization_strengths={'Gx': 0.05},
-                                                 stochastic_error_probs={'Gy': (0.01, 0.02, 0.03)},
-                                                 lindblad_error_coeffs={'Gcnot': {('H','ZZ'): 0.07, ('S','XX'): 0.10}},
-                                                 independent_gates=False, independent_spam=True, verbosity=2)
-        self.check_model(mdl_cloud)
-
-        mdl_cloud = create_cloud_crosstalk_model(self.pspec_2Q, depolarization_strengths={'Gx': 0.05},
-                                                 stochastic_error_probs={'Gy': (0.01, 0.02, 0.03)},
-                                                 lindblad_error_coeffs={'Gcnot': {('H','ZZ'): 0.07, ('S','XX'): 0.10}},
-                                                 independent_gates=True, independent_spam=True, verbosity=2)
-        self.check_model(mdl_cloud)
\ No newline at end of file
diff --git a/test/unit/objects/test_model.py b/test/unit/objects/test_model.py
index 66835a1d3..ee9369f80 100644
--- a/test/unit/objects/test_model.py
+++ b/test/unit/objects/test_model.py
@@ -15,6 +15,9 @@
 from pygsti.models import ExplicitOpModel
 from pygsti.circuits import Circuit
 from pygsti.models.gaugegroup import FullGaugeGroupElement
+from pygsti.modelpacks import smq1Q_XYI
+from pygsti.processors import QubitProcessorSpec
+from pygsti.models import create_crosstalk_free_model, create_cloud_crosstalk_model
 from ..util import BaseCase, needs_cvxpy
 SKIP_DIAMONDIST_ON_WIN = True
 
@@ -210,7 +213,7 @@ def test_jtracedist(self):
 
     @needs_cvxpy
     def test_diamonddist(self):
-        if SKIP_DIAMONDIST_ON_WIN and sys.platform.startswith('win'): return
+        if  SKIP_DIAMONDIST_ON_WIN and sys.platform.startswith('win'): return
         cp = self.model.copy()
         self.assertAlmostEqual(self.model.diamonddist(cp), 0)
         # TODO non-trivial case
@@ -743,3 +746,55 @@ def tearDown(self):
 #    def test_randomize_with_unitary_raises(self):
 #        with self.assertRaises(AssertionError):
 #            self.model.randomize_with_unitary(1, rand_state=np.random.RandomState())  # scale shouldn't matter
+
+class ModelEquivalenceTester(BaseCase):
+
+    def setUp(self):
+        nQubits = 2
+        self.pspec_2Q = QubitProcessorSpec(nQubits, ('Gx', 'Gy', 'Gcnot'), geometry="line",
+                                           qubit_labels=['qb{}'.format(i) for i in range(nQubits)])
+
+    def check_model(self, mdl):
+        mcopy = mdl.copy()
+        self.assertFalse(mdl is mcopy)
+        self.assertTrue(mcopy.is_similar(mdl))
+        self.assertTrue(mcopy.is_equivalent(mdl))
+
+        if mdl.num_params > 0:
+            r = np.random.random(mdl.num_params)
+            if np.linalg.norm(r) < 1e6:  # just in case we randomly get all zeros!
+                r = 0.1 * np.ones(mcopy.num_params)
+            v_prime = mdl.to_vector() + r
+            mcopy.from_vector(v_prime)
+            self.assertTrue(mcopy.is_similar(mdl))
+            self.assertFalse(mcopy.is_equivalent(mdl))
+
+    def test_explicit_model_equal(self):
+        mdl_explicit = smq1Q_XYI.target_model()
+        self.check_model(mdl_explicit)
+
+    def test_local_model_equal(self):
+        mdl_local = create_crosstalk_free_model(self.pspec_2Q,
+                                                ideal_gate_type='H+S', ideal_spam_type='tensor product H+S',
+                                                independent_gates=False,
+                                                ensure_composed_gates=False)
+        self.check_model(mdl_local)
+
+        mdl_local = create_crosstalk_free_model(self.pspec_2Q,
+                                                ideal_gate_type='H+S', ideal_spam_type='computational',
+                                                independent_gates=True,
+                                                ensure_composed_gates=True)
+        self.check_model(mdl_local)
+
+    def test_cloud_model_equal(self):
+        mdl_cloud = create_cloud_crosstalk_model(self.pspec_2Q, depolarization_strengths={'Gx': 0.05},
+                                                 stochastic_error_probs={'Gy': (0.01, 0.02, 0.03)},
+                                                 lindblad_error_coeffs={'Gcnot': {('H','ZZ'): 0.07, ('S','XX'): 0.10}},
+                                                 independent_gates=False, independent_spam=True, verbosity=2)
+        self.check_model(mdl_cloud)
+
+        mdl_cloud = create_cloud_crosstalk_model(self.pspec_2Q, depolarization_strengths={'Gx': 0.05},
+                                                 stochastic_error_probs={'Gy': (0.01, 0.02, 0.03)},
+                                                 lindblad_error_coeffs={'Gcnot': {('H','ZZ'): 0.07, ('S','XX'): 0.10}},
+                                                 independent_gates=True, independent_spam=True, verbosity=2)
+        self.check_model(mdl_cloud)
\ No newline at end of file

From 3f6650463fd37c361b77cd587c083eb4aac3f1d8 Mon Sep 17 00:00:00 2001
From: Juan G Mendoza <jgonz19@sandia.gov>
Date: Wed, 9 Jul 2025 10:41:53 -0600
Subject: [PATCH 24/34] added __eq__ docstring

---
 pygsti/baseobjs/errorgenbasis.py | 14 ++++++++++++--
 pygsti/tools/fogitools.py        |  2 +-
 2 files changed, 13 insertions(+), 3 deletions(-)

diff --git a/pygsti/baseobjs/errorgenbasis.py b/pygsti/baseobjs/errorgenbasis.py
index 276316bbf..d737e8c8e 100644
--- a/pygsti/baseobjs/errorgenbasis.py
+++ b/pygsti/baseobjs/errorgenbasis.py
@@ -120,6 +120,16 @@ def __init__(self, state_space, labels, basis_1q=None):
         self._cached_supports = None
 
     def __eq__(self, other):
+        """Compare self with other. Return true only if they are identical, including the order of their labels.
+
+        Args:
+            other (ExplicitElementaryErrorgenBasis): Error generator basis to compare against
+
+        Returns:
+            Boolean: True if they are identical, False otherwise
+        """
+        if not isinstance(other, ExplicitElementaryErrorgenBasis):
+            return False
         return self.state_space.__eq__(other.state_space) and [label.__str__() for label in self.labels] == [label.__str__() for label in other.labels] and self._basis_1q.__eq__(other._basis_1q)
     
     def _to_nice_serialization(self):
@@ -802,8 +812,8 @@ def label_index(self, label, ok_if_missing=False, identity_label='I'):
         elif eetype in ('C', 'A'):
             assert(len(trivial_bel) == 1)  # assumes this is a single character
             nontrivial_inds = [i for i, letter in enumerate(bels[0]) if letter != trivial_bel]
-            left_support = tuple([label.sslbls[i] for i in nontrivial_inds])
-            #left_support = tuple([self.sslbls[i] for i in nontrivial_inds])
+            #left_support = tuple([label.sslbls[i] for i in nontrivial_inds])
+            left_support = tuple([self.sslbls[i] for i in nontrivial_inds])
             #This line above is supposed to be a bugfix, I want to verify this with
             #Corey before removing this comment
             if ok_if_missing and (support, left_support) not in self._offsets[eetype]:
diff --git a/pygsti/tools/fogitools.py b/pygsti/tools/fogitools.py
index e51c3add6..96093a690 100644
--- a/pygsti/tools/fogitools.py
+++ b/pygsti/tools/fogitools.py
@@ -70,7 +70,7 @@ def _embed(mx, target_labels, state_space):  # BOTTLENECK
         # - a full basis for gauge_action_deriv
         #global_row_space.add_labels(row_space.labels)  # labels would need to contain sslbls too
         action_row_labels = action_row_basis.labels
-        global_row_indices = elemgen_row_basis.label_indices(action_row_labels, ok_if_missing=True)
+        global_row_indices = elemgen_row_basis.label_indices(action_row_labels, ok_if_missing=True) 
 
         #DEBUG REMOVE
         #db_num_skipped = db_num_nonzero = 0

From 965c7856f8af28aec990163a686673b7f15b8305 Mon Sep 17 00:00:00 2001
From: Juan G Mendoza <jgonz19@sandia.gov>
Date: Wed, 9 Jul 2025 10:43:48 -0600
Subject: [PATCH 25/34] added another dosctring

---
 pygsti/baseobjs/errorgenspace.py | 9 +++++++++
 1 file changed, 9 insertions(+)

diff --git a/pygsti/baseobjs/errorgenspace.py b/pygsti/baseobjs/errorgenspace.py
index de949a221..aaf223bdc 100644
--- a/pygsti/baseobjs/errorgenspace.py
+++ b/pygsti/baseobjs/errorgenspace.py
@@ -94,7 +94,16 @@ def intersection(self, other_space, free_on_unspecified_space=False, use_nice_nu
         return ErrorgenSpace(intersection_vecs, common_basis)
 
     def __eq__(self, other):
+        """Compare self with other. Return true only if they are identical, including the order of their vectors.
 
+        Args:
+            other (ExplicitElementaryErrorgenBasis): Error generator basis to compare against
+
+        Returns:
+            Boolean: True if they are identical, False otherwise
+        """
+        if not isinstance(other, ErrorgenSpace):
+            return False
         return _np.allclose(self.vectors, other.vectors) and self.elemgen_basis.__eq__(other.elemgen_basis)
     def union(self, other_space):
         """

From b1078dcfdca71a2f1795e34e6833e7113f7b4f29 Mon Sep 17 00:00:00 2001
From: Juan G Mendoza <jgonz19@sandia.gov>
Date: Wed, 9 Jul 2025 10:46:12 -0600
Subject: [PATCH 26/34] added another docstring

---
 pygsti/baseobjs/errorgenspace.py |  2 +-
 pygsti/models/fogistore.py       | 10 +++++++++-
 2 files changed, 10 insertions(+), 2 deletions(-)

diff --git a/pygsti/baseobjs/errorgenspace.py b/pygsti/baseobjs/errorgenspace.py
index aaf223bdc..2bfe256c7 100644
--- a/pygsti/baseobjs/errorgenspace.py
+++ b/pygsti/baseobjs/errorgenspace.py
@@ -97,7 +97,7 @@ def __eq__(self, other):
         """Compare self with other. Return true only if they are identical, including the order of their vectors.
 
         Args:
-            other (ExplicitElementaryErrorgenBasis): Error generator basis to compare against
+            other (ErrorgenSpace): Error generator space to compare against
 
         Returns:
             Boolean: True if they are identical, False otherwise
diff --git a/pygsti/models/fogistore.py b/pygsti/models/fogistore.py
index 5483f95f8..2e5503e9d 100644
--- a/pygsti/models/fogistore.py
+++ b/pygsti/models/fogistore.py
@@ -395,8 +395,16 @@ def find_nice_fogiv_directions(self):
         pass
     
     def __eq__(self, other):
+        """Compare self with other. Return true only if they are identical, including the order of their attribute elements, such as labels.
 
-        assert isinstance(other, FirstOrderGaugeInvariantStore), 'Object provided is not of type FirstOrderGaugeInvariantStore'
+        Args:
+            other (FirstOrderGaugeInvariantStore): FOGIStore to compare against self
+
+        Returns:
+            Boolean: True if they are identical, False otherwise
+        """
+        if not  isinstance(other, FirstOrderGaugeInvariantStore):
+            return False
 
         if [label.__str__() for label in self.primitive_op_labels] != [label.__str__() for label in other.primitive_op_labels]:
             return False

From e392502fd03d066f5bc6b44f5635f0b88706292e Mon Sep 17 00:00:00 2001
From: Juan G Mendoza <jgonz19@sandia.gov>
Date: Wed, 9 Jul 2025 10:47:21 -0600
Subject: [PATCH 27/34] removed whitespace

---
 test/unit/objects/test_model.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/unit/objects/test_model.py b/test/unit/objects/test_model.py
index ee9369f80..7ad43e293 100644
--- a/test/unit/objects/test_model.py
+++ b/test/unit/objects/test_model.py
@@ -213,7 +213,7 @@ def test_jtracedist(self):
 
     @needs_cvxpy
     def test_diamonddist(self):
-        if  SKIP_DIAMONDIST_ON_WIN and sys.platform.startswith('win'): return
+        if SKIP_DIAMONDIST_ON_WIN and sys.platform.startswith('win'): return
         cp = self.model.copy()
         self.assertAlmostEqual(self.model.diamonddist(cp), 0)
         # TODO non-trivial case

From 67deafded4310174256abddaf159a6039a076f16 Mon Sep 17 00:00:00 2001
From: Juan G Mendoza <jgonz19@sandia.gov>
Date: Wed, 9 Jul 2025 11:37:16 -0600
Subject: [PATCH 28/34] improved label_index tests

---
 pygsti/baseobjs/errorgenbasis.py        |  4 +--
 test/unit/objects/test_errorgenbasis.py | 40 ++++++++++++++++++++++---
 2 files changed, 38 insertions(+), 6 deletions(-)

diff --git a/pygsti/baseobjs/errorgenbasis.py b/pygsti/baseobjs/errorgenbasis.py
index d737e8c8e..a5ae55c7c 100644
--- a/pygsti/baseobjs/errorgenbasis.py
+++ b/pygsti/baseobjs/errorgenbasis.py
@@ -812,8 +812,8 @@ def label_index(self, label, ok_if_missing=False, identity_label='I'):
         elif eetype in ('C', 'A'):
             assert(len(trivial_bel) == 1)  # assumes this is a single character
             nontrivial_inds = [i for i, letter in enumerate(bels[0]) if letter != trivial_bel]
-            #left_support = tuple([label.sslbls[i] for i in nontrivial_inds])
-            left_support = tuple([self.sslbls[i] for i in nontrivial_inds])
+            left_support = tuple([label.sslbls[i] for i in nontrivial_inds])
+            #left_support = tuple([self.sslbls[i] for i in nontrivial_inds])
             #This line above is supposed to be a bugfix, I want to verify this with
             #Corey before removing this comment
             if ok_if_missing and (support, left_support) not in self._offsets[eetype]:
diff --git a/test/unit/objects/test_errorgenbasis.py b/test/unit/objects/test_errorgenbasis.py
index 4d2f1da7f..e5cb55762 100644
--- a/test/unit/objects/test_errorgenbasis.py
+++ b/test/unit/objects/test_errorgenbasis.py
@@ -12,6 +12,7 @@ def setUp(self):
         
         #create a complete basis with default settings for reuse.
         self.complete_errorgen_basis_default_1Q = CompleteElementaryErrorgenBasis(self.basis_1q, self.state_space_1Q)
+        self.complete_errorgen_basis_default_2Q = CompleteElementaryErrorgenBasis(self.basis_1q, self.state_space_2Q)
 
     def test_default_construction(self):
         assert len(self.complete_errorgen_basis_default_1Q.labels) == 12
@@ -104,20 +105,51 @@ def test_elemgen_and_dual_construction(self):
         duals = self.complete_errorgen_basis_default_1Q.elemgen_dual_matrices
 
     def test_label_index(self):
+        
+        #1 qubit tests
         labels = self.complete_errorgen_basis_default_1Q.labels
+        test_eg = GlobalElementaryErrorgenLabel('A', ['X', 'Y'], (0,))
+        lbl_idx = self.complete_errorgen_basis_default_1Q.label_index(test_eg)
+        assert lbl_idx ==  labels.index(test_eg)
 
-        test_eg = GlobalElementaryErrorgenLabel('C', ['X', 'Y'], (0,))
-        test_eg_local = LocalElementaryErrorgenLabel('C', ['XI', 'YI'])
+        #Test missing label
         test_eg_missing = GlobalElementaryErrorgenLabel('C', ['X', 'Y'], (1,))
 
+        with self.assertRaises(KeyError):
+            self.complete_errorgen_basis_default_1Q.label_index(test_eg_missing)
+        assert self.complete_errorgen_basis_default_1Q.label_index(test_eg_missing, ok_if_missing=True) is None
+
+        #Test embedding 
+        test_eg = GlobalElementaryErrorgenLabel('C', ['X', 'Y'], (0,))
+        test_eg_local = LocalElementaryErrorgenLabel('C', ['XI', 'YI'])
+        
         lbl_idx = self.complete_errorgen_basis_default_1Q.label_index(test_eg)
         lbl_idx_1 = self.complete_errorgen_basis_default_1Q.label_index(test_eg_local)
         assert lbl_idx == lbl_idx_1
         assert lbl_idx ==  labels.index(test_eg)
 
+        # 2 qubit tests
+        labels = self.complete_errorgen_basis_default_2Q.labels
+
+        test_eg = GlobalElementaryErrorgenLabel('A', ['X', 'Y'], (0,))
+        lbl_idx = self.complete_errorgen_basis_default_2Q.label_index(test_eg)
+        assert lbl_idx ==  labels.index(test_eg)
+
+        #Test missing label
+        test_eg_missing = GlobalElementaryErrorgenLabel('C', ['X', 'Y'], (2,))
+
         with self.assertRaises(KeyError):
-            self.complete_errorgen_basis_default_1Q.label_index(test_eg_missing)
-        assert self.complete_errorgen_basis_default_1Q.label_index(test_eg_missing, ok_if_missing=True) is None
+            self.complete_errorgen_basis_default_2Q.label_index(test_eg_missing)
+        assert self.complete_errorgen_basis_default_2Q.label_index(test_eg_missing, ok_if_missing=True) is None
+
+        #Test embedding for 2 qubit labels
+        test_eg = GlobalElementaryErrorgenLabel('C', ['X', 'Y'], (1,))
+        test_eg_local = LocalElementaryErrorgenLabel('C', ['IX', 'IY'])
+
+        lbl_idx = self.complete_errorgen_basis_default_2Q.label_index(test_eg)
+        lbl_idx_1 = self.complete_errorgen_basis_default_2Q.label_index(test_eg_local)
+        assert lbl_idx == lbl_idx_1
+        assert lbl_idx ==  labels.index(test_eg)
 
     def test_create_subbasis(self):
         errorgen_basis = CompleteElementaryErrorgenBasis(self.basis_1q, self.state_space_2Q)

From 6784b870f7f70107c684a1d87fc4a9ab951faf7b Mon Sep 17 00:00:00 2001
From: Juan G Mendoza <jgonz19@sandia.gov>
Date: Wed, 9 Jul 2025 12:49:43 -0600
Subject: [PATCH 29/34] possible fix for when svd does not converge

---
 pygsti/tools/fogitools.py | 34 ++++++++++++++++++++++++++++++++--
 1 file changed, 32 insertions(+), 2 deletions(-)

diff --git a/pygsti/tools/fogitools.py b/pygsti/tools/fogitools.py
index 96093a690..e59ed6fa9 100644
--- a/pygsti/tools/fogitools.py
+++ b/pygsti/tools/fogitools.py
@@ -560,9 +560,39 @@ def resolve_norm_order(vecs_to_normalize, label_lists, given_norm_order):
                         # (A,B are faithful reps of gauge on intersection space, so pinv(ga_A) * ga_A
                         # restricted to intersection space is I:   int_spc.T * pinv(ga_A) * ga_A * int_spc == I
                         # (when int_spc vecs are orthonormal) and so the above redues to I - I = 0
+                        def scipy_pinv(A, rcond):
+                            # Step 1: Compute SVD
+                            U, S, VT = _spsl.svd(A, lapack_driver='gesvd')
 
-                        inv_diff_gauge_action = _np.concatenate((_np.linalg.pinv(gauge_action[0:n, :], rcond=1e-7),
-                                                                 -_np.linalg.pinv(gauge_action[n:, :], rcond=1e-7)),
+                            # Step 2: Construct Sigma
+                            Sigma = _np.zeros((U.shape[0], VT.shape[0]))
+                            _np.fill_diagonal(Sigma, S)
+
+                            # Step 3: Compute Sigma+
+                            S_plus = _np.zeros_like(Sigma.T)
+                            for i in range(len(S)):
+                                if _np.abs(S[i]) > rcond:
+                                    S_plus[i, i] = 1 / S[i]
+
+                            # Step 4: Compute the pseudoinverse
+                            return VT.T @ S_plus @ U.T
+
+                        try:
+                            inverse1 = _np.linalg.pinv(gauge_action[0:n, :], rcond=1e-7)
+                        except Exception as e:
+                            print('Pinverse failed with message: ', e, '\n Attempting gesvd algorithm for SVD')
+                            inverse1 = scipy_pinv(gauge_action[0:n, :], rcond=1e-7)
+                            print('success')
+                        
+                        try:
+                            inverse2 = -_np.linalg.pinv(gauge_action[n:, :], rcond=1e-7)
+                        except Exception as e:
+                            print('Pinverse failed with message: ', e, '\n Attempting gesvd algorithm for SVD')
+                            inverse2 = -scipy_pinv(gauge_action[n:, :], rcond=1e-7)
+                            print('success')
+
+                        inv_diff_gauge_action = _np.concatenate((inverse1,
+                                                                 inverse2),
                                                                 axis=1).T
 
                         #Equivalent:

From 2df848ae08ee779ea8615ac182b20a820c735ea5 Mon Sep 17 00:00:00 2001
From: Juan G Mendoza <jgonz19@sandia.gov>
Date: Mon, 14 Jul 2025 10:17:27 -0600
Subject: [PATCH 30/34] fixed import

---
 pygsti/tools/fogitools.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/pygsti/tools/fogitools.py b/pygsti/tools/fogitools.py
index e59ed6fa9..cdba367cf 100644
--- a/pygsti/tools/fogitools.py
+++ b/pygsti/tools/fogitools.py
@@ -13,7 +13,7 @@
 import numpy as _np
 import scipy.sparse as _sps
 import scipy.sparse.linalg as _spsl
-
+import scipy as _scp
 from . import matrixtools as _mt
 from . import optools as _ot
 
@@ -562,7 +562,7 @@ def resolve_norm_order(vecs_to_normalize, label_lists, given_norm_order):
                         # (when int_spc vecs are orthonormal) and so the above redues to I - I = 0
                         def scipy_pinv(A, rcond):
                             # Step 1: Compute SVD
-                            U, S, VT = _spsl.svd(A, lapack_driver='gesvd')
+                            U, S, VT = _scp.linalg.svd(A, lapack_driver='gesvd')
 
                             # Step 2: Construct Sigma
                             Sigma = _np.zeros((U.shape[0], VT.shape[0]))

From 3b401b21fe670270640eb53db44687c3fc37ccd1 Mon Sep 17 00:00:00 2001
From: Juan G Mendoza <jgonz19@sandia.gov>
Date: Wed, 16 Jul 2025 15:43:11 -0600
Subject: [PATCH 31/34] removed commented out code

---
 pygsti/models/fogistore.py | 17 -----------------
 1 file changed, 17 deletions(-)

diff --git a/pygsti/models/fogistore.py b/pygsti/models/fogistore.py
index 2e5503e9d..06035e7ad 100644
--- a/pygsti/models/fogistore.py
+++ b/pygsti/models/fogistore.py
@@ -459,24 +459,7 @@ def __eq__(self, other):
         if self.norm_order != other.norm_order or self._dependent_fogi_action != other._dependent_fogi_action:
             return False
         return True
-        
-        '''
-        self.primitive_op_labels = primitive_op_labels
-        self.gauge_space = gauge_space
-        self.elem_errorgen_labels_by_op = elem_errorgen_labels_by_op
-        self.op_errorgen_indices = op_errorgen_indices
-        self.fogi_directions = fogi_directions
-        self.fogi_metadata = fogi_metadata
-        self.dependent_dir_indices = dependent_dir_indices
-        self.fogv_directions = fogv_directions
-        self.allop_gauge_action = allop_gauge_action
-        self.gauge_space_directions = gauge_space_directions
-        self.norm_order = norm_order
-        self._dependent_fogi_action = dependent_fogi_action
 
-        self.errorgen_space_op_elem_labels = tuple([(op_label, elem_lbl) for op_label in self.primitive_op_labels
-                                                    for elem_lbl in self.elem_errorgen_labels_by_op[op_label]])
-        self.fogv_labels = ["%d gauge action" % i for i in range(self.fogv_directions.shape[1])]'''
     @property
     def errorgen_space_dim(self):
         return self.fogi_directions.shape[0]

From f17d24a14ec993b8ff6b9c4a82c39ce3b6d469cb Mon Sep 17 00:00:00 2001
From: Juan G Mendoza <jgonz19@sandia.gov>
Date: Fri, 18 Jul 2025 16:48:03 -0600
Subject: [PATCH 32/34] 1000x Murray update

---
 pygsti/tools/fogitools.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pygsti/tools/fogitools.py b/pygsti/tools/fogitools.py
index cdba367cf..058a7c793 100644
--- a/pygsti/tools/fogitools.py
+++ b/pygsti/tools/fogitools.py
@@ -712,7 +712,7 @@ def scipy_pinv(A, rcond):
                         new_fogi_dirs = new_fogi_dirs.tocsc()
 
                         # figure out which directions are independent
-                        indep_cols = _mt.independent_columns(new_fogi_dirs, fogi_dirs)
+                        indep_cols = _mt.independent_columns(new_fogi_dirs.toarray(), fogi_dirs.toarray())
                         #FOGI DEBUG print(" ==> %d independent columns" % len(indep_cols))
 
                         if dependent_fogi_action == "drop":

From 2171e3be28600a21762814cc88a3770706d9663a Mon Sep 17 00:00:00 2001
From: Corey Ostrove <ciostro@sandia.gov>
Date: Mon, 28 Jul 2025 15:15:24 -0600
Subject: [PATCH 33/34] Fix kwarg

Correct the kwarg on qr decomposition.
---
 pygsti/tools/matrixtools.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pygsti/tools/matrixtools.py b/pygsti/tools/matrixtools.py
index 8fcd29d3c..c6fbc29c9 100644
--- a/pygsti/tools/matrixtools.py
+++ b/pygsti/tools/matrixtools.py
@@ -420,7 +420,7 @@ def independent_columns(m, initial_independent_cols=None, tol=1e-7):
             # We assume initial_independent_cols is full column-rank.
             # This lets us use unpivoted QR instead of pivoted QR or SVD.
             assert initial_independent_cols.shape[0] == m.shape[0]
-            q = _spl.qr(initial_independent_cols, mode='econ')[0]
+            q = _spl.qr(initial_independent_cols, mode='economic')[0]
             # proj_m = (I - qq')m
             temp1 = q.T.conj() @ m
             temp2 = q @ temp1

From b68f18297c8a893c4d0ec8d8a4e1030066f06375 Mon Sep 17 00:00:00 2001
From: Juan G Mendoza <jgonz19@sandia.gov>
Date: Mon, 4 Aug 2025 10:27:44 -0600
Subject: [PATCH 34/34] removed unnecessary diagonal matrix as pointed out by
 Corey

---
 pygsti/tools/fogitools.py   | 10 +++-------
 pygsti/tools/matrixtools.py |  2 +-
 2 files changed, 4 insertions(+), 8 deletions(-)

diff --git a/pygsti/tools/fogitools.py b/pygsti/tools/fogitools.py
index 058a7c793..584cc331b 100644
--- a/pygsti/tools/fogitools.py
+++ b/pygsti/tools/fogitools.py
@@ -564,17 +564,13 @@ def scipy_pinv(A, rcond):
                             # Step 1: Compute SVD
                             U, S, VT = _scp.linalg.svd(A, lapack_driver='gesvd')
 
-                            # Step 2: Construct Sigma
-                            Sigma = _np.zeros((U.shape[0], VT.shape[0]))
-                            _np.fill_diagonal(Sigma, S)
-
-                            # Step 3: Compute Sigma+
-                            S_plus = _np.zeros_like(Sigma.T)
+                            # Step 2: Compute Sigma+
+                            S_plus = _np.zeros((VT.shape[0],U.shape[0]))
                             for i in range(len(S)):
                                 if _np.abs(S[i]) > rcond:
                                     S_plus[i, i] = 1 / S[i]
 
-                            # Step 4: Compute the pseudoinverse
+                            # Step 3: Compute the pseudoinverse
                             return VT.T @ S_plus @ U.T
 
                         try:
diff --git a/pygsti/tools/matrixtools.py b/pygsti/tools/matrixtools.py
index 8fcd29d3c..c6fbc29c9 100644
--- a/pygsti/tools/matrixtools.py
+++ b/pygsti/tools/matrixtools.py
@@ -420,7 +420,7 @@ def independent_columns(m, initial_independent_cols=None, tol=1e-7):
             # We assume initial_independent_cols is full column-rank.
             # This lets us use unpivoted QR instead of pivoted QR or SVD.
             assert initial_independent_cols.shape[0] == m.shape[0]
-            q = _spl.qr(initial_independent_cols, mode='econ')[0]
+            q = _spl.qr(initial_independent_cols, mode='economic')[0]
             # proj_m = (I - qq')m
             temp1 = q.T.conj() @ m
             temp2 = q @ temp1