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/fastcircuitparser.pyx

@@ -236,6 +236,11 @@ cdef get_next_simple_lbl(unicode s, INT start, INT end, bool integerize_sslbls,
     else:
         while i < end:
             c = s[i]
+            if c != 'G':
+                # We need to convert to lowercase in case there are strings like
+                # "GXI" "GXX", "GIGY", or "GCNOT". We don't convert 'G' because
+                # that's a special character for our purposes.
+                c = c.lower()
             if u'a' <= c <= u'z' or u'0' <= c <= u'9' or c == u'_':
                 i += 1
             else:
@@ -250,7 +255,11 @@ cdef get_next_simple_lbl(unicode s, INT start, INT end, bool integerize_sslbls,
             c = s[i]
             if u'0' <= c <= u'9' or c == u'.' or c == u'-' or c == u'e': #last case for scientific notation
                 i += 1
-            elif u'a' <= c <= u'z' or c == u'_' or c == u'Q' or c == u'/':
+                continue
+            if c != 'G':
+                # We convert to lowercase here for the same reason as above.
+                c = c.lower()
+            if u'a' <= c <= u'z' or c == u'_' or c == u'Q' or c == u'/':
                 i += 1; is_float = False
             else:
                 break
@@ -262,6 +271,9 @@ cdef get_next_simple_lbl(unicode s, INT start, INT end, bool integerize_sslbls,
         last = i; is_int = True
         while i < end:
             c = s[i]
+            if c != 'G':
+                # We convert to lowercase here for the same reason as above.
+                c = c.lower()
             if u'0' <= c <= u'9':
                 i += 1
             elif u'a' <= c <= u'z' or c == u'_' or c == u'Q':

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
@@ -171,6 +171,11 @@ def _get_next_simple_lbl(s, start, end, integerize_sslbls, segment):
     else:
         while i < end:
             c = s[i]
+            if c != 'G':
+                # We need to convert to lowercase in case there are strings like
+                # "GXI" "GXX", "GIGY", or "GCNOT". We don't convert 'G' because
+                # that's a special character for our purposes.
+                c = c.lower()
             if 'a' <= c <= 'z' or '0' <= c <= '9' or c == '_':
                 i += 1
             else:
@@ -183,6 +188,9 @@ def _get_next_simple_lbl(s, start, end, integerize_sslbls, segment):
         last = i
         while i < end:
             c = s[i]
+            if c != 'G':
+                # We convert to lowercase here for the same reason as above.
+                c = c.lower()
             if 'a' <= c <= 'z' or '0' <= c <= '9' or c == '_' or c == 'Q' or c == '.' or c == '/' or c == '-':
                 i += 1
             else:
@@ -199,6 +207,9 @@ def _get_next_simple_lbl(s, start, end, integerize_sslbls, segment):
         last = i
         while i < end:
             c = s[i]
+            if c != 'G':
+                # We convert to lowercase here for the same reason as above.
+                c = c.lower()
             if 'a' <= c <= 'z' or '0' <= c <= '9' or c == '_' or c == 'Q':
                 i += 1
             else:

pygsti/data/datasetconstruction.py

@@ -209,6 +209,39 @@ def simulate_data(model_or_dataset, circuit_list, num_samples,
     return dataset
 
 
+def mix_datasets(dsa, dsb, p, integral=True, choose_machine=False, seed=None):
+    dsc = dsa.copy_nonstatic()
+    # arr = _np.array(dsc.repData).ravel()
+    # print((arr, arr.size))
+    # print((dsb.repData, dsb.repData.size))
+    num_circuits = len(dsb)
+    if choose_machine:
+        if seed is None:
+            _warnings.warn('Set the random seed! Using 42.')
+            seed = 42
+        rngstate = _np.random.default_rng(seed)
+        interp_weights = rngstate.uniform(low=0, high=1, size=num_circuits)
+        interp_weights[interp_weights < p] = 0.0
+        interp_weights[interp_weights > 0] = 1.0 
+    else:
+        interp_weights = p * _np.ones(num_circuits)
+
+    for i, (_, dsrow) in enumerate(dsb.items()):
+        p_i = interp_weights[i]
+        interpolated = p_i * dsc.repData[i] + (1-p_i) * dsrow.reps
+        if integral and (not choose_machine):
+            assert interpolated.size == 2
+            total = int(_np.ceil((_np.sum(interpolated))))
+            j = _np.argmin(interpolated)
+            interpolated[j] = _np.ceil(interpolated[j])
+            interpolated[1 - j] = total - interpolated[j]
+        dsc.repData[i][:] = interpolated
+    dsc.done_adding_data()
+    # arr = np.array(dsc.repData).ravel()
+    # print((arr, arr.size))
+    return dsc
+
+
 def _adjust_probabilities_inbounds(ps, tol):
     #Adjust to probabilities if needed (and warn if not close to in-bounds)
     # ps is a dict w/keys = outcome labels and values = probabilities

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/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) """