WIP: TVD and other objectives in GST

.gitignore

@@ -1,3 +1,4 @@
+wip_notebook_sharing/case0_reports_*
 *~
 *.tmp
 *.bak

pygsti/algorithms/gaugeopt.py

@@ -335,12 +335,21 @@ def _call_jacobian_fn(gauge_group_el_vec):
         if bToStdout and (comm is None or comm.Get_rank() == 0):
             print_obj_func = _opt.create_objfn_printer(_call_objective_fn)  # only ever prints to stdout!
             # print_obj_func(x0) #print initial point (can be a large vector though)
-        else: print_obj_func = None
+        else:
+            print_obj_func = None
+
+        minimize_kwargs = {'method': method, 'maxiter': maxiter, 'maxfev': maxfev, 'tol': tol}
+        minimize_kwargs['jac'] = '3-point' if _call_jacobian_fn is None else _call_jacobian_fn
 
-        minSol = _opt.minimize(_call_objective_fn, x0,
-                               method=method, maxiter=maxiter, maxfev=maxfev,
-                               tol=tol, jac=_call_jacobian_fn,
-                               callback=print_obj_func)
+        minSol = _opt.minimize(_call_objective_fn, x0, callback=print_obj_func, **minimize_kwargs)
+        if not minSol.success:
+            msg = f"""\n
+            Gauge optimization failed to converge! Algorithm parameters were
+            {minimize_kwargs}.
+
+            The optimizer returned with message "{minSol.message}".
+            """
+            _warnings.warn(msg)
         solnX = minSol.x
         solnF = minSol.fun
 

pygsti/circuits/circuitparser/slowcircuitparser.py

@@ -65,7 +65,7 @@ def parse_circuit(code, create_subcircuits=True, integerize_sslbls=True):
     return tuple(result), labels, occurrence_id, compilable_indices
 
 
-def parse_label(code, integerize_sslbls=True):
+def parse_label(code: str, integerize_sslbls=True) -> _lbl.Label:
     create_subcircuits = False
     segment = 0  # segment for gates/instruments vs. preps vs. povms: 0 = *any*
     interlayer_marker = u''  # matches nothing - no interlayer markerg

pygsti/drivers/longsequence.py

@@ -25,7 +25,7 @@
 from pygsti.models.modelconstruction import _create_explicit_model, create_explicit_model
 from pygsti.protocols.gst import _load_pspec_or_model
 from pygsti.forwardsims import ForwardSimulator
-from typing import Optional
+from typing import Optional, Any
 
 ROBUST_SUFFIX_LIST = [".robust", ".Robust", ".robust+", ".Robust+"]
 DEFAULT_BAD_FIT_THRESHOLD = 2.0
@@ -315,7 +315,7 @@ def run_linear_gst(data_filename_or_set, target_model_filename_or_object,
 def run_long_sequence_gst(data_filename_or_set, target_model_filename_or_object,
                           prep_fiducial_list_or_filename, meas_fiducial_list_or_filename,
                           germs_list_or_filename, max_lengths, gauge_opt_params=None,
-                          advanced_options=None, comm=None, mem_limit=None,
+                          advanced_options: Optional[dict[str,Any]]=None, comm=None, mem_limit=None,
                           output_pkl=None, verbosity=2, checkpoint=None, checkpoint_path=None,
                           disable_checkpointing=False,
                           simulator: Optional[ForwardSimulator.Castable]=None,
@@ -386,7 +386,12 @@ def run_long_sequence_gst(data_filename_or_set, target_model_filename_or_object,
         Specifies advanced options most of which deal with numerical details of
         the objective function or expert-level functionality.  The allowed keys
         and values include:
-        - objective = {'chi2', 'logl'}
+
+        HISTORICAL NOTE: "XX" indicates that we've at least _intended_ for the
+        keyword argument to be removed.
+
+        - objective = typically, a string in {'chi2', 'logl', 'tvd'}. But this can
+          be anything accepted by the `ObjectiveFunctionBuilder.create_from` function.
         - op_labels = list of strings
         - circuit_weights = dict or None
         - starting_point = "LGST-if-possible" (default), "LGST", or "target"
@@ -402,19 +407,16 @@ def run_long_sequence_gst(data_filename_or_set, target_model_filename_or_object,
         - prob_clip_interval = tuple (default == (-1e6,1e6)
         - radius = float (default == 1e-4)
         - use_freq_weighted_chi2 = True / False (default)
-        - XX nested_circuit_lists = True (default) / False
-        - XX include_lgst = True / False (default is True)
+        - XX nested_circuit_lists = True (default) / False  -- passed to StandardGSTDesign; seems to be functional
+        - XX include_lgst = True / False (default is True)  -- passed to StandardGSTDesign; seems to be functional
         - distribute_method = "default", "circuits" or "deriv"
         - profile = int (default == 1)
         - check = True / False (default)
-        - XX op_label_aliases = dict (default = None)
+        - XX op_label_aliases = dict (default = None) -- passed to StandardGSTDesign and used to set a GateSetTomography protocol object's oplabel_aliases field.
         - always_perform_mle = bool (default = False)
         - only_perform_mle = bool (default = False)
-        - XX truncScheme = "whole germ powers" (default) or "truncated germ powers" or "length as exponent"
         - appendTo = Results (default = None)
-        - estimateLabel = str (default = "default")
-        - XX missingDataAction = {'drop','raise'} (default = 'drop')
-        - XX string_manipulation_rules = list of (find,replace) tuples
+        - XX string_manipulation_rules = list of (find,replace) tuples -- passed to _proto.StandardGSTDesign construct as its "circuit_rules" argument.
         - germ_length_limits = dict of form {germ: maxlength}
         - record_output = bool (default = True)
         - timeDependent = bool (default = False)

pygsti/layouts/copalayout.py

@@ -202,13 +202,14 @@ def __init__(self, circuits, unique_circuits, to_unique, elindex_outcome_tuples,
             "Inconsistency: %d distinct indices but max index + 1 is %d!" % (len(indices), self._size)
 
         self._outcomes = dict()
-        self._element_indices = dict()
+        self._element_indices : dict[tuple, slice] = dict()
         sort_idx_func = lambda x: x[0]
         for i_unique, tuples in elindex_outcome_tuples.items():
             sorted_tuples = sorted(tuples, key=sort_idx_func)  # sort by element index
             elindices, outcomes = zip(*sorted_tuples)  # sorted by elindex so we make slices whenever possible
             self._outcomes[i_unique] = tuple(outcomes)
-            self._element_indices[i_unique] = _slct.list_to_slice(elindices, array_ok=True)
+            s = _slct.list_to_slice(elindices, array_ok=True)
+            self._element_indices[i_unique] = s # type: ignore
 
 #    def hotswap_circuits(self, circuits, unique_complete_circuits=None):
 #        self.circuits = circuits if isinstance(circuits, _CircuitList) else _CircuitList(circuits)
@@ -740,7 +741,9 @@ def indices_and_outcomes_for_index(self, index):
 
     def __iter__(self):
         for circuit, i in self._unique_circuit_index.items():
-            for element_index, outcome in zip(self._element_indices[i], self._outcomes[i]):
+            indices = _slct.to_array(self._element_indices[i])
+            iterator = zip(indices, self._outcomes[i])
+            for element_index, outcome in iterator:
                 yield element_index, circuit, outcome
 
     def iter_unique_circuits(self):

pygsti/modelmembers/modelmember.py

@@ -309,9 +309,20 @@ def relink_parent(self, parent):
         """
         for subm in self.submembers():
             subm.relink_parent(parent)
-
-        if self._parent is parent: return  # OK to relink multiple times
-        assert(self._parent is None), "Cannot relink parent: parent is not None!"
+        #
+        # We need to either set self._parent or raise an error if it's set incorrectly.
+        # This is delicate, because ...
+        #
+        #   1. The `_parent`` attribute might be missing from self.__dict__.
+        #
+        #   2. type(self) might have its own __getattr__ method, and the logic of
+        #      that method might break if its called when the _parent attribute
+        #      is missing.
+        #
+        if '_parent' in self.__dict__:
+            if self._parent is parent:
+                return  # OK to relink multiple times
+            assert(self._parent is None), "Cannot relink parent: current parent is not None!"
         self._parent = parent  # assume no dependent objects
 
     def unlink_parent(self, force=False):

pygsti/modelmembers/modelmembergraph.py

@@ -40,14 +40,37 @@ def load_modelmembers_from_serialization_dict(sdict, parent_model):
 
         for mm_node_serialized_id_str, mm_node_dict in sdict.items():
             mm_node_serialized_id = int(mm_node_serialized_id_str)  # convert string keys back to integers
-            mm_class = ModelMember._state_class(mm_node_dict)  # checks this is a ModelMember-derived class
+            mm_class = ModelMember._state_class(mm_node_dict)       # checks this is a ModelMember-derived class
             mm_serial[mm_node_serialized_id] = mm_class.from_memoized_dict(mm_node_dict, mm_serial, parent_model)
-            if 'memberdict_types' in mm_node_dict and 'memberdict_labels' in mm_node_dict:
-                for mm_type, lbl_str in zip(mm_node_dict['memberdict_types'], mm_node_dict['memberdict_labels']):
-                    lbl = _parse_label(lbl_str)
-                    if mm_type not in mm_nodes:
-                        mm_nodes[mm_type] = {}
-                    mm_nodes[mm_type][lbl] = mm_serial[mm_node_serialized_id]
+
+            md_types  = mm_node_dict.get('memberdict_types',  dict())
+            md_labels = mm_node_dict.get('memberdict_labels', dict())
+
+            assert len(md_types) == len(md_labels)
+
+            for mm_type, lbl_str in zip(md_types, md_labels):
+                lbl = _parse_label(lbl_str)
+                if mm_type not in mm_nodes:
+                    mm_nodes[mm_type] = {}
+                d = mm_nodes[mm_type]
+                if lbl in d:
+                    msg = f"""
+                    We've encountered a collision during deserialization. The 
+                    string-valued {mm_type} label "{lbl_str}" was mapped to the
+                    Label object `{lbl}`, but data has already been recorded for
+                    a(n) {mm_type} with this Label.
+
+                    This can happen when a string-valued modelmember label doesn't
+                    follow formatting rules required by the parse_labels function
+                    in pygsti.circuits.circuitparser.
+
+                    We're raising an error because we can't gaurantee correctness
+                    of the result if we proceeded.
+                    """
+                    raise RuntimeError(msg)
+                    # ^ Another approach would be to raise an error when `lbl_str not in lbl`
+                    #   if `lbl = parse_label(lbl_str)` is a LabelStr.
+                d[lbl] = mm_serial[mm_node_serialized_id]
 
         return mm_nodes
 

pygsti/modelmembers/povms/composedpovm.py

@@ -22,8 +22,8 @@
 
 class ComposedPOVM(_POVM):
     """
-    TODO: update docstring
-    A POVM that is effectively a *single* Lindblad-parameterized gate followed by a computational-basis POVM.
+    A parameterized POVM that is effectively a *single* parameterized gate followed by
+    a static (usually, computational basis) POVM.
 
     Parameters
     ----------
@@ -37,7 +37,7 @@ class ComposedPOVM(_POVM):
     povm : POVM, optional
         A sub-POVM which supplies the set of "reference" effect vectors
         that `errormap` acts on to produce the final effect vectors of
-        this LindbladPOVM.  This POVM must be *static*
+        this ComposedPOVM.  This POVM must be *static*
         (have zero parameters) and its evolution type must match that of
         `errormap`.  If None, then a :class:`ComputationalBasisPOVM` is
         used on the number of qubits appropriate to `errormap`'s dimension.
@@ -50,59 +50,34 @@ class ComposedPOVM(_POVM):
     """
 
     def __init__(self, errormap, povm=None, mx_basis=None):
-        """
-        Creates a new LindbladPOVM object.
-
-        Parameters
-        ----------
-        errormap : MapOperator
-            The error generator action and parameterization, encapsulated in
-            a gate object.  Usually a :class:`LindbladOp`
-            or :class:`ComposedOp` object.  (This argument is *not* copied,
-            to allow ComposedPOVMEffects to share error generator
-            parameters with other gates and spam vectors.)
-
-        povm : POVM, optional
-            A sub-POVM which supplies the set of "reference" effect vectors
-            that `errormap` acts on to produce the final effect vectors of
-            this LindbladPOVM.  This POVM must be *static*
-            (have zero parameters) and its evolution type must match that of
-            `errormap`.  If None, then a :class:`ComputationalBasisPOVM` is
-            used on the number of qubits appropriate to `errormap`'s dimension.
-
-        mx_basis : {'std', 'gm', 'pp', 'qt'} or Basis object
-            The basis for this spam vector. Allowed values are Matrix-unit (std),
-            Gell-Mann (gm), Pauli-product (pp), and Qutrit (qt) (or a custom
-            basis object).  If None, then this is extracted (if possible) from
-            `errormap`.
-        """
         self.error_map = errormap
-        state_space = self.error_map.state_space
+        state_space = errormap.state_space
+        evotype = errormap._evotype
+
+        if povm is None:
+            povm = _ComputationalBasisPOVM(state_space.num_qubits, evotype)
+        elif isinstance(povm, ComposedPOVM):
+            from pygsti.modelmembers.operations import ComposedOp
+            errormap = ComposedOp([povm.errormap, errormap])
+            povm = povm.base_povm
+
+        assert(povm.evotype == evotype), \
+            ("Evolution type of `povm` (%s) must match that of "
+                "`errormap` (%s)!") % (povm.evotype, evotype)
+        assert(povm.num_params == 0), \
+            "Given `povm` must be static (have 0 parameters)!"
+
+        self.base_povm = povm
 
         if mx_basis is None:
-            if (isinstance(errormap, (_op.ExpErrorgenOp, _op.IdentityPlusErrorgenOp))
-                and isinstance(errormap.errorgen, _op.LindbladErrorgen)):
-                mx_basis = errormap.errorgen.matrix_basis
+            if hasattr(errormap, 'errorgen') and isinstance(errormap.errorgen, _op.LindbladErrorgen): # type: ignore
+                mx_basis = errormap.errorgen.matrix_basis # type: ignore
             else:
-                raise ValueError("Cannot extract a matrix-basis from `errormap` (type %s)"
-                                 % str(type(errormap)))
+                raise ValueError(f"Cannot extract a matrix-basis from `errormap` (type {type(errormap)})")
 
         self.matrix_basis = _Basis.cast(mx_basis, state_space)
-        evotype = self.error_map._evotype
-
-        if povm is None:
-            assert(state_space.num_qubits >= 0), \
-                ("A default computational-basis POVM can only be used with an"
-                 " integral number of qubits!")
-            povm = _ComputationalBasisPOVM(state_space.num_qubits, evotype)
-        else:
-            assert(povm.evotype == evotype), \
-                ("Evolution type of `povm` (%s) must match that of "
-                 "`errormap` (%s)!") % (povm.evotype, evotype)
-            assert(povm.num_params == 0), \
-                "Given `povm` must be static (have 0 parameters)!"
-        self.base_povm = povm
 
+        self.errormap = errormap
         items = []  # init as empty (lazy creation of members)
         try:
             rep = evotype.create_composed_povm_rep(self.error_map._rep, self.base_povm._rep, state_space)
@@ -188,7 +163,7 @@ def __getitem__(self, key):
             assert(self.parent is None or num_new == 0)  # ensure effect inds are already allocated to current model
             _collections.OrderedDict.__setitem__(self, key, effect)
             return effect
-        else: raise KeyError("%s is not an outcome label of this LindbladPOVM" % key)
+        else: raise KeyError("%s is not an outcome label of this ComposedPOVM" % key)
 
     def __reduce__(self):
         """ Needed for OrderedDict-derived classes (to set dict items) """

pygsti/modelmembers/states/composedstate.py

@@ -19,21 +19,25 @@
 from pygsti.modelmembers.errorgencontainer import ErrorMapContainer as _ErrorMapContainer
 from pygsti import SpaceT
 
-class ComposedState(_State):  # , _ErrorMapContainer
+class ComposedState(_State):
     """
-    TODO: update docstring
-    A Lindblad-parameterized State (that is also expandable into terms).
+    A representation of the superket `errormap @ state_vec` induced by a superoperator
+    `errormap` and a superket `state_vec`, where `state_vec.num_params == 0`. 
+
 
     Parameters
     ----------
-    pure_vec : numpy array or State
+    static_state : numpy array or StaticState or ComposedState
         An array or State in the *full* density-matrix space (this
-        vector will have dimension 4 in the case of a single qubit) which
-        represents a pure-state preparation or projection.  This is used as
-        the "base" preparation or projection that is followed or preceded
-        by, respectively, the parameterized Lindblad-form error generator.
-        (This argument is *not* copied if it is a State.  A numpy array
-        is converted to a new StaticState.)
+        vector will have dimension 4 in the case of a single qubit).
+
+        This argument is *not* copied if it is a State. If it's an
+        ndarray then we convert it to a new StaticState. If it's a
+        ComposedState, then we use a ComposedOp to pack `errormap`
+        and `static_state.errormap` into a single superoperator, and
+        we replace `static_state = static_state.state_vec`.
+
+        After possible conversions, we set `self.state_vec = static_state`.
 
     errormap : MapOperator
         The error generator action and parameterization, encapsulated in
@@ -43,21 +47,23 @@ class ComposedState(_State):  # , _ErrorMapContainer
         parameters with other gates and spam vectors.)
     """
 
+    STATIC_STATE_ERROR_MSG = "`static_state` 'reference' must have *zero* parameters!"
+
     def __init__(self, static_state, errormap):
         evotype = errormap._evotype
-        #from .operation import LindbladOp as _LPGMap
-        #assert(evotype in ("densitymx", "svterm", "cterm")), \
-        #    "Invalid evotype: %s for %s" % (evotype, self.__class__.__name__)
-
         if not isinstance(static_state, _State):
             # UNSPECIFIED BASIS - change None to static_state.basis once we have a std attribute
             static_state = _StaticState(static_state, None, evotype)  # assume spamvec is just a vector
 
+        if isinstance(static_state, ComposedState):
+            from pygsti.modelmembers.operations import ComposedOp
+            errormap = ComposedOp([static_state.error_map, errormap])
+            static_state = static_state.state_vec
+
         assert(static_state._evotype == evotype), \
             "`static_state` evotype must match `errormap` ('%s' != '%s')" % (static_state._evotype, evotype)
-        assert(static_state.num_params == 0), "`static_state` 'reference' must have *zero* parameters!"
+        assert(static_state.num_params == 0), ComposedState.STATIC_STATE_ERROR_MSG
 
-        #d2 = static_state.dim
         self.state_vec = static_state
         self.error_map = errormap
         self.terms = {}
@@ -67,7 +73,7 @@ def __init__(self, static_state, errormap):
         rep = evotype.create_composed_state_rep(self.state_vec._rep, self.error_map._rep, static_state.state_space)
 
         _State.__init__(self, rep, evotype)
-        _ErrorMapContainer.__init__(self, self.error_map)
+        _ErrorMapContainer.__init__(self, self.error_map) # type: ignore
         self.init_gpindices()  # initialize our gpindices based on sub-members
 
     @classmethod