diff --git a/.gitignore b/.gitignore index c366c2c6f..d983388e8 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,4 @@ +wip_notebook_sharing/case0_reports_* *~ *.tmp *.bak diff --git a/pygsti/algorithms/gaugeopt.py b/pygsti/algorithms/gaugeopt.py index c332553fa..a63d73822 100644 --- a/pygsti/algorithms/gaugeopt.py +++ b/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 diff --git a/pygsti/circuits/circuitparser/fastcircuitparser.pyx b/pygsti/circuits/circuitparser/fastcircuitparser.pyx index ad0767906..48671266b 100644 --- a/pygsti/circuits/circuitparser/fastcircuitparser.pyx +++ b/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': diff --git a/pygsti/circuits/circuitparser/slowcircuitparser.py b/pygsti/circuits/circuitparser/slowcircuitparser.py index c4a2d8bcb..5b3733808 100644 --- a/pygsti/circuits/circuitparser/slowcircuitparser.py +++ b/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: diff --git a/pygsti/data/datasetconstruction.py b/pygsti/data/datasetconstruction.py index bcf9b796b..7ecd95c57 100644 --- a/pygsti/data/datasetconstruction.py +++ b/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 diff --git a/pygsti/drivers/longsequence.py b/pygsti/drivers/longsequence.py index 5b2fc3503..f97240b9f 100644 --- a/pygsti/drivers/longsequence.py +++ b/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) diff --git a/pygsti/layouts/copalayout.py b/pygsti/layouts/copalayout.py index f2ea0d8e6..d1bda42dd 100644 --- a/pygsti/layouts/copalayout.py +++ b/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): diff --git a/pygsti/modelmembers/modelmembergraph.py b/pygsti/modelmembers/modelmembergraph.py index 5b2b2e743..99a7f4711 100644 --- a/pygsti/modelmembers/modelmembergraph.py +++ b/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 diff --git a/pygsti/modelmembers/povms/composedpovm.py b/pygsti/modelmembers/povms/composedpovm.py index 9f2f93bcd..968c4eb11 100644 --- a/pygsti/modelmembers/povms/composedpovm.py +++ b/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) """ diff --git a/pygsti/modelmembers/states/composedstate.py b/pygsti/modelmembers/states/composedstate.py index fd5af4220..1edeb8304 100644 --- a/pygsti/modelmembers/states/composedstate.py +++ b/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 diff --git a/pygsti/models/explicitmodel.py b/pygsti/models/explicitmodel.py index 5d276d9e8..cc1cb9aeb 100644 --- a/pygsti/models/explicitmodel.py +++ b/pygsti/models/explicitmodel.py @@ -239,6 +239,15 @@ def default_gauge_group(self, value): """ The default gauge group. """ + if isinstance(value, str): + value = value.lower() + if value == 'unitary': + value = _gg.UnitaryGaugeGroup(self.state_space, self.basis, self.evotype) + elif value == 'tp': + value = _gg.TPGaugeGroup(self.state_space, self.basis, self.evotypes) + else: + msg = f'string "{value}" cannot be used to set this model\'s gauge group.' + raise ValueError(msg) self._default_gauge_group = value @property @@ -983,6 +992,14 @@ def all_objects(self): self.factories.items()): yield (lbl, obj) + def members(self, include=('preps', 'operations', 'povms')): + out = dict() + for member_type in include: + member_dict = self.__getattribute__(member_type) + for k,v in member_dict.items(): + out[k] = v + return out + #TODO: how to handle these given possibility of different parameterizations... # -- maybe only allow these methods to be called when using a "full" parameterization? # -- or perhaps better to *move* them to the parameterization class @@ -1157,7 +1174,7 @@ def rotate(self, rotate=None, max_rotate=None, seed=None): newModel._clean_paramvec() # rotate may leave dirty members return newModel - def randomize_with_unitary(self, scale, seed=None, rand_state=None): + def randomize_with_unitary(self, scale, seed=None, rand_state=None, transform_spam=False, param_preserving=False): """ Create a new model with random unitary perturbations. @@ -1181,6 +1198,12 @@ def randomize_with_unitary(self, scale, seed=None, rand_state=None): across multiple random function calls but you don't want to bother with manually incrementing seeds between those calls. + param_preserving : bool + If False, then all members of the returned model will have "full" parameterization. + + If True, then we use ComposedOp (and ComposedState and ComposedPOVM) to define a + model with unitarily-transformed members while retaining the parameterization of `self`. + Returns ------- Model @@ -1198,18 +1221,55 @@ def randomize_with_unitary(self, scale, seed=None, rand_state=None): mdl_randomized = self.copy() - for opLabel, gate in self.operations.items(): - randMat = scale * (rndm.randn(unitary_dim, unitary_dim) - + 1j * rndm.randn(unitary_dim, unitary_dim)) - randMat = _np.transpose(_np.conjugate(randMat)) + randMat - # make randMat Hermetian: (A_dag + A)^dag = (A_dag + A) - randUnitary = _scipy.linalg.expm(-1j * randMat) - - randOp = _ot.unitary_to_superop(randUnitary, self.basis) + def rand_unitary_as_superop(): + rand_mat = rndm.randn(unitary_dim, unitary_dim) + 1j * rndm.randn(unitary_dim, unitary_dim) + rand_herm = rand_mat.T.conj() + rand_mat + rand_herm /= _scipy.linalg.norm(rand_herm) + rand_herm *= scale * _np.sqrt(unitary_dim) + rand_unitary = _scipy.linalg.expm(-1j * rand_herm) + rand_op = _op.StaticUnitaryOp(rand_unitary, self.basis) + return rand_op + + if param_preserving: + from pygsti.modelmembers.states import ComposedState + from pygsti.modelmembers.povms import ComposedPOVM + def transformed_gate(rand_op, gate): # type: ignore + ops_to_compose = gate.factorops if hasattr(gate, 'factorops') else [gate] + ops_to_compose.append(rand_op) + return _op.ComposedOp(ops_to_compose) + def transformed_stateprep(rand_op, rho): # type: ignore + return ComposedState(rho, rand_op) + def transformed_povm(rand_op, M): # type: ignore + return ComposedPOVM(rand_op, M) + else: + from pygsti.modelmembers.states import FullState + from pygsti.modelmembers.povms import create_from_dmvecs + def transformed_gate(rand_op, gate): + rand_op = rand_op.to_dense() + return _op.FullArbitraryOp(rand_op @ gate) + def transformed_stateprep(rand_op, rho): + rand_op = rand_op.to_dense() + return FullState(rand_op @ rho) + def transformed_povm(rand_op, M): + rand_op = rand_op.to_dense() + dmvecs = {elbl: rand_op @ e.to_dense() for elbl, e in M.items()} + return create_from_dmvecs(dmvecs, 'full') - mdl_randomized.operations[opLabel] = _op.FullArbitraryOp(_np.dot(randOp, gate.to_dense("HilbertSchmidt"))) + for opLabel, gate in self.operations.items(): + rand_op = rand_unitary_as_superop() + mdl_randomized.operations[opLabel] = transformed_gate(rand_op, gate) - #Note: this function does NOT randomize instruments + if transform_spam: + + for preplbl, rho in self.preps.items(): + rand_op = rand_unitary_as_superop() + mdl_randomized.preps[preplbl] = transformed_stateprep(rand_op, rho) + + for povmlbl, M in self.povms.items(): + rand_op = rand_unitary_as_superop() + mdl_randomized.povms[povmlbl] = transformed_povm(rand_op, M) + + # Note: this function does NOT randomize instruments return mdl_randomized diff --git a/pygsti/models/model.py b/pygsti/models/model.py index 5f6f8e940..cea0bcbf8 100644 --- a/pygsti/models/model.py +++ b/pygsti/models/model.py @@ -2320,7 +2320,14 @@ def copy(self): """ self._clean_paramvec() # ensure _paramvec is rebuilt if needed if OpModel._pcheck: self._check_paramvec() - ret = Model.copy(self) + try: + ret = Model.copy(self) + except: + # This can do a better job than `Model.copy` at resolving ambiguities + # in setting the `parent` field of Modelmembers in the copied model. + state = self.to_nice_serialization() + ret = type(self).from_nice_serialization(state) + if self._param_bounds is not None and self.parameter_labels is not None: ret._clean_paramvec() # will *always* rebuild paramvec; do now so we can preserve param bounds assert _np.all(self.parameter_labels == ret.parameter_labels) # ensure ordering is the same diff --git a/pygsti/models/modelconstruction.py b/pygsti/models/modelconstruction.py index 42c4763d2..3b8748c71 100644 --- a/pygsti/models/modelconstruction.py +++ b/pygsti/models/modelconstruction.py @@ -567,8 +567,26 @@ def _create_explicit_model_from_expressions(state_space, basis, if len(effect_vecs) > 0: # don't add POVMs with 0 effects ret.povms[povmLbl] = _povm.create_from_dmvecs(effect_vecs, povm_type, basis, evotype, state_space) - for (opLabel, opExpr) in zip(op_labels, op_expressions): - ret.operations[opLabel] = create_operation(opExpr, state_space, basis, gate_type, evotype) + from pygsti.circuits.circuitparser import parse_label + + parsed_op_labels = [parse_label(str(opLabel)) for opLabel in op_labels] + if len(set(parsed_op_labels)) != len(op_labels): + msg = f""" + There are fewer unique Label objects after parsing op_labels than + there are elements in op_labels. If we proceeded with the current + op_labels then we would not be able to return a model that could + be serialized and subequently deserialized (specifically, + deserialization would fail). Since all pyGSTi Model objects implement + the NicelySerializable API, we're raising an error. + + The initial op_labels are {op_labels}. + + The parsed op_labels are {parsed_op_labels}. + """ + raise ValueError(msg) + + for (lbl, opExpr) in zip(parsed_op_labels, op_expressions): + ret.operations[lbl] = create_operation(opExpr, state_space, basis, gate_type, evotype) if gate_type == "full": ret.default_gauge_group = _gg.FullGaugeGroup(ret.state_space, basis, evotype) diff --git a/pygsti/objectivefns/objectivefns.py b/pygsti/objectivefns/objectivefns.py index a88975f16..880079df4 100644 --- a/pygsti/objectivefns/objectivefns.py +++ b/pygsti/objectivefns/objectivefns.py @@ -15,6 +15,7 @@ import time as _time import pathlib as _pathlib import warnings as _warnings +from typing import Callable import numpy as _np @@ -29,6 +30,13 @@ from pygsti.models.model import OpModel as _OpModel from pygsti import SpaceT +def set_docstring(docstr): + def assign(fn): + fn.__doc__ = docstr + return fn + return assign + + def _objfn(objfn_cls, model, dataset, circuits=None, regularization=None, penalties=None, op_label_aliases=None, comm=None, mem_limit=None, method_names=None, array_types=None, @@ -143,8 +151,9 @@ class ObjectiveFunctionBuilder(_NicelySerializable): Penalty values (allowed keys depend on `cls_to_build`). """ - @classmethod - def cast(cls, obj): + # This was a classmethod, but I made it static because this class has no derived classes. + @staticmethod + def cast(obj): """ Cast `obj` to an `ObjectiveFunctionBuilder` instance. @@ -160,6 +169,7 @@ def cast(cls, obj): ------- ObjectiveFunctionBuilder """ + cls = ObjectiveFunctionBuilder if isinstance(obj, cls): return obj elif obj is None: return cls.create_from() elif isinstance(obj, str): return cls.create_from(objective=obj) @@ -167,15 +177,16 @@ def cast(cls, obj): elif isinstance(obj, (list, tuple)): return cls(*obj) else: raise ValueError("Cannot create an %s object from '%s'" % (cls.__name__, str(type(obj)))) - @classmethod - def create_from(cls, objective='logl', freq_weighted_chi2=False): + # This was a classmethod, but I made it static because this class has no derived classes. + @staticmethod + def create_from(objective='logl', freq_weighted_chi2=False, callable_penalty=None): """ Creates common :class:`ObjectiveFunctionBuilder` from a few arguments. Parameters ---------- - objective : {'logl', 'chi2'}, optional - The objective function type: log-likelihood or chi-squared. + objective : {'logl', 'chi2', 'tvd'}, optional + The objective function type: log-likelihood, chi-squared, or TVD. freq_weighted_chi2 : bool, optional Whether to use 1/frequency values as the weights in the `"chi2"` case. @@ -184,35 +195,63 @@ def create_from(cls, objective='logl', freq_weighted_chi2=False): ------- ObjectiveFunctionBuilder """ + if callable_penalty is None and isinstance(objective, tuple) and len(objective) >= 2 and isinstance(objective[1], Callable): + callable_penalty = objective[1] + # callable_penalty_scale = 1.0 if len(objective) == 2 else objective[2] + objective = objective[0] + + if objective == "chi2": if freq_weighted_chi2: - builder = FreqWeightedChi2Function.builder( + builder = ObjectiveFunctionBuilder(FreqWeightedChi2Function, name='fwchi2', description="Freq-weighted sum of Chi^2", - regularization={'min_freq_clip_for_weighting': 1e-4}) + regularization={'min_freq_clip_for_weighting': 1e-4}, + callable_penalty=callable_penalty + ) else: - builder = Chi2Function.builder( + builder = ObjectiveFunctionBuilder(Chi2Function, name='chi2', description="Sum of Chi^2", - regularization={'min_prob_clip_for_weighting': 1e-4}) + regularization={'min_prob_clip_for_weighting': 1e-4}, + callable_penalty=callable_penalty + ) elif objective == "logl": - builder = PoissonPicDeltaLogLFunction.builder( + builder = ObjectiveFunctionBuilder(PoissonPicDeltaLogLFunction, name='dlogl', description="2*Delta(log(L))", regularization={'min_prob_clip': 1e-4, 'radius': 1e-4}, penalties={'cptp_penalty_factor': 0, - 'spam_penalty_factor': 0}) - - elif objective == "tvd": - builder = TVDFunction.builder( - name='tvd', - description="Total Variational Distance (TVD)") + 'spam_penalty_factor': 0}, + callable_penalty=callable_penalty + ) + + elif 'tvd' in objective: + descr = "Total Variational Distance (TVD)" + if 'normalized' in objective: + descr = descr + ', normalized by circuit depth' + assert objective == 'normalized tvd' + else: + assert objective == 'tvd' + builder = ObjectiveFunctionBuilder( + TVDFunction, name=objective, description=descr, + callable_penalty=callable_penalty + ) + + elif isinstance(objective, tuple) and objective[0] == 'Lp^p': + power = objective[1] + builder = ObjectiveFunctionBuilder(LpNormToPowerP, + name='Lp^p', + description=f"L_{power} norm to the power {power}.", + power=objective[1], + callable_penalty=callable_penalty + ) else: raise ValueError("Invalid objective: %s" % objective) - assert(isinstance(builder, cls)), "This function should always return an ObjectiveFunctionBuilder!" + return builder def __init__(self, cls_to_build, name=None, description=None, regularization=None, penalties=None, **kwargs): @@ -555,6 +594,8 @@ def lsvec(self, probs, counts, total_counts, freqs, intermediates=None): """ return _np.sqrt(self.terms(probs, counts, total_counts, freqs, intermediates)) + # Infinite loop in evaluation of "dterms" and "dlsvec". + def dterms(self, probs, counts, total_counts, freqs, intermediates=None): """ Compute the derivatives of the terms of this objective function. @@ -589,10 +630,11 @@ def dterms(self, probs, counts, total_counts, freqs, intermediates=None): """ if intermediates is None: intermediates = self._intermediates(probs, counts, total_counts, freqs) - lsvec = self.lsvec(probs, counts, total_counts, freqs, intermediates) + u = self.lsvec(probs, counts, total_counts, freqs, intermediates) + v = self.dlsvec(probs, counts, total_counts, freqs, intermediates) # NOTE: this function is correct under the assumption that terms == lsvec**2, # independent of whether or not lsvec is nonnegative. - return 2 * lsvec * self.dlsvec(probs, counts, total_counts, freqs, intermediates) + return 2 * u * v def dlsvec(self, probs, counts, total_counts, freqs, intermediates=None): """ @@ -627,7 +669,8 @@ def dlsvec(self, probs, counts, total_counts, freqs, intermediates=None): A 1D array of length equal to that of each array argument. """ # lsvec = sqrt(terms) - # NOTE: ^ That's only correct if lsvec is >= 0, and some classes don't satisfy that. + # NOTE: ^ That's only correct if lsvec is >= 0. + # Any class that doesn't ensure lsvec >= 0 must override this function. # dlsvec = 0.5/lsvec * dterms # if intermediates is None: @@ -1095,9 +1138,9 @@ class MDCObjectiveFunction(ObjectiveFunction, EvaluatedModelDatasetCircuitsStore """ @classmethod - def create_from(cls, raw_objfn, model, dataset, circuits, resource_alloc=None, verbosity=0, array_types=()): + def create_from(cls, raw_objfn, model, dataset, circuits, resource_alloc=None, verbosity=0, array_types=(), **kwargs): mdc_store = ModelDatasetCircuitsStore(model, dataset, circuits, resource_alloc, array_types) - return cls(raw_objfn, mdc_store, verbosity) + return cls(raw_objfn, mdc_store, verbosity, **kwargs) @classmethod def _array_types_for_method(cls, method_name, fsim): @@ -1109,7 +1152,7 @@ def _array_types_for_method(cls, method_name, fsim): if method_name == 'dpercircuit': return cls._array_types_for_method('dterms', fsim) + ('cp',) return () - def __init__(self, raw_objfn, mdc_store, verbosity=0): + def __init__(self, raw_objfn, mdc_store, verbosity=0, **kwargs): """ Create a new MDCObjectiveFunction. @@ -1118,6 +1161,8 @@ def __init__(self, raw_objfn, mdc_store, verbosity=0): """ EvaluatedModelDatasetCircuitsStore.__init__(self, mdc_store, verbosity) self.raw_objfn = raw_objfn + self.callable_penalty = kwargs.get('callable_penalty', None) + self.callable_penalty_factor = 0.0 if self.callable_penalty is None else 1.0 @property def name(self): @@ -1356,7 +1401,7 @@ def dlsvec_percircuit(self, paramvec=None): # Note: don't need paramvec here since above call sets it return (0.5 / denom)[:, None] * self.dpercircuit() - def fn_local(self, paramvec=None): + def fn_local(self, paramvec=None, stateless=False): """ Evaluate the *local* value of this objective function. @@ -1375,9 +1420,15 @@ def fn_local(self, paramvec=None): ------- float """ - return _np.sum(self.terms(paramvec)) + if paramvec is None or not stateless: + return _np.sum(self.terms(paramvec)) + else: + old_paramvec = self.model.to_vector() + val = _np.sum(self.terms(paramvec)) + self.model.from_vector(old_paramvec) + return val - def fn(self, paramvec=None): + def fn(self, paramvec=None, stateless=False): """ Evaluate the value of this objective function. @@ -1392,13 +1443,20 @@ def fn(self, paramvec=None): float """ result, result_shm = _smt.create_shared_ndarray(self.resource_alloc, (1,), 'd') - local = _np.array([self.fn_local(paramvec)], 'd') + local = _np.array([self.fn_local(paramvec, stateless)], 'd') unit_ralloc = self.layout.resource_alloc('atom-processing') # proc group that computes same els self.resource_alloc.allreduce_sum(result, local, unit_ralloc) global_fnval = result[0] _smt.cleanup_shared_ndarray(result_shm) return global_fnval + def fn_from_model(self, model): + m = self.model + self.model = model + val = self.fn() + self.model = m + return val + def jacobian(self, paramvec=None): """ Compute the Jacobian of this objective function. @@ -2671,36 +2729,37 @@ def zero_freq_hterms(self, total_counts, probs): return 2 * _np.ones(len(probs)) -# The log(Likelihood) within the Poisson picture is: # noqa -# # noqa -# L = prod_{i,sl} lambda_{i,sl}^N_{i,sl} e^{-lambda_{i,sl}} / N_{i,sl}! # noqa -# # noqa -# Where lamba_{i,sl} := p_{i,sl}*N[i] is a rate, i indexes the operation sequence, # noqa -# and sl indexes the spam label. N[i] is the total counts for the i-th circuit, and # noqa -# so sum_{sl} N_{i,sl} == N[i]. We can ignore the p-independent N_j! and take the log: # noqa -# # noqa -# log L = sum_{i,sl} N_{i,sl} log(N[i]*p_{i,sl}) - N[i]*p_{i,sl} # noqa -# = sum_{i,sl} N_{i,sl} log(p_{i,sl}) - N[i]*p_{i,sl} (where we ignore the p-independent log(N[i]) terms) # noqa -# # noqa -# The objective function computes the negative log(Likelihood) as a vector of leastsq # noqa -# terms, where each term == sqrt( N_{i,sl} * -log(p_{i,sl}) + N[i] * p_{i,sl} ) # noqa -# # noqa -# See LikelihoodFunctions.py for details on patching # noqa -# The log(Likelihood) within the standard picture is: -# -# L = prod_{i,sl} p_{i,sl}^N_{i,sl} -# -# Where i indexes the operation sequence, and sl indexes the spam label. -# N[i] is the total counts for the i-th circuit, and -# so sum_{sl} N_{i,sl} == N[i]. We take the log: -# -# log L = sum_{i,sl} N_{i,sl} log(p_{i,sl}) -# -# The objective function computes the negative log(Likelihood) as a vector of leastsq -# terms, where each term == sqrt( N_{i,sl} * -log(p_{i,sl}) ) -# -# See LikelihoodFunction.py for details on patching - +""" +The log(Likelihood) within the Poisson picture is: # noqa + # noqa +L = prod_{i,sl} lambda_{i,sl}^N_{i,sl} e^{-lambda_{i,sl}} / N_{i,sl}! # noqa + # noqa +Where lamba_{i,sl} := p_{i,sl}*N[i] is a rate, i indexes the operation sequence, # noqa + and sl indexes the spam label. N[i] is the total counts for the i-th circuit, and # noqa + so sum_{sl} N_{i,sl} == N[i]. We can ignore the p-independent N_j! and take the log: # noqa + # noqa +log L = sum_{i,sl} N_{i,sl} log(N[i]*p_{i,sl}) - N[i]*p_{i,sl} # noqa + = sum_{i,sl} N_{i,sl} log(p_{i,sl}) - N[i]*p_{i,sl} (where we ignore the p-independent log(N[i]) terms) # noqa + # noqa +The objective function computes the negative log(Likelihood) as a vector of leastsq # noqa + terms, where each term == sqrt( N_{i,sl} * -log(p_{i,sl}) + N[i] * p_{i,sl} ) # noqa + # noqa +See LikelihoodFunctions.py for details on patching # noqa +The log(Likelihood) within the standard picture is: + +L = prod_{i,sl} p_{i,sl}^N_{i,sl} + +Where i indexes the operation sequence, and sl indexes the spam label. + N[i] is the total counts for the i-th circuit, and + so sum_{sl} N_{i,sl} == N[i]. We take the log: + +log L = sum_{i,sl} N_{i,sl} log(p_{i,sl}) + +The objective function computes the negative log(Likelihood) as a vector of leastsq + terms, where each term == sqrt( N_{i,sl} * -log(p_{i,sl}) ) + +See LikelihoodFunction.py for details on patching +""" class RawPoissonPicDeltaLogLFunction(RawObjectiveFunction): """ @@ -3938,6 +3997,9 @@ def __init__(self, regularization=None, resource_alloc=None, name='tvd', description="Total Variational Distance (TVD)", verbosity=0): super().__init__(regularization, resource_alloc, name, description, verbosity) + def chi2k_distributed_qty(self, objective_function_value): + return -1 + def terms(self, probs, counts, total_counts, freqs, intermediates=None): """ Compute the terms of the objective function. @@ -4003,7 +4065,10 @@ def dterms(self, probs, counts, total_counts, freqs, intermediates=None): numpy.ndarray A 1D array of length equal to that of each array argument. """ - raise NotImplementedError("Derivatives not implemented for TVD yet!") + t = probs - freqs + d = 0.5*_np.ones_like(t) + d[t < 0] *= -1 + return d def hterms(self, probs, counts, total_counts, freqs, intermediates=None): """ @@ -4113,6 +4178,30 @@ def zero_freq_hterms(self, total_counts, probs): raise NotImplementedError("Derivatives not implemented for TVD yet!") +class RawAbsPower(RawObjectiveFunction): + + def __init__(self, power: float, regularization=None, resource_alloc=None, + name='Lp^p', description="Elementwise absolute value and raising to a power.", verbosity=0): + super().__init__(regularization, resource_alloc, name, description, verbosity) + assert power >= 1 + self.power = power + + def chi2k_distributed_qty(self, objective_function_value): + return -1 + + def terms(self, probs, counts, total_counts, freqs, intermediates=None): + return _np.abs(probs - freqs) ** self.power + + def dterms(self, probs, counts, total_counts, freqs, intermediates=None): + t = probs - freqs + d = self.power * _np.abs(t) ** (self.power - 1) + d[t < 0] *= -1 + return d + + def zero_freq_terms(self, total_counts, probs): + return _np.abs(probs) ** self.power + + ###################################################### # # Start MDCObjectiveFunction subclasses @@ -4121,15 +4210,23 @@ def zero_freq_hterms(self, total_counts, probs): class TimeIndependentMDCObjectiveFunction(MDCObjectiveFunction): + + TEMPLATE_FIELDS = ( """ A time-independent model-based (:class:`MDCObjectiveFunction`-derived) objective function. - - Parameters - ---------- + """, + """ raw_objfn : RawObjectiveFunction The raw objective function - specifies how probability and count values are turned into objective function values. + """, "" + ) + DOCSTR_TEMPLATE = \ + """ + %s + Parameters + ----------%s mdl : Model The model - specifies how parameter values are turned into probabilities for each circuit outcome. @@ -4165,33 +4262,9 @@ class TimeIndependentMDCObjectiveFunction(MDCObjectiveFunction): Whether hessian calculations are allowed. If `True` then more resources are needed. If `False`, calls to hessian-requiring function will result in an error. + %s """ - @classmethod - def builder(cls, name=None, description=None, regularization=None, penalties=None, **kwargs): - """ - Create an :class:`ObjectiveFunctionBuilder` that builds an objective function of this type. - - Parameters - ---------- - name : str, optional - A name for the built objective function (can be anything). - - description : str, optional - A description for the built objective function (can be anything) - - regularization : dict, optional - Regularization values. - - penalties : dict, optional - Penalty values. - - Returns - ------- - ObjectiveFunctionBuilder - """ - return ObjectiveFunctionBuilder(cls, name, description, regularization, penalties, **kwargs) - @classmethod def _create_mdc_store(cls, model, dataset, circuits, resource_alloc, method_names=('fn',), array_types=(), verbosity=0): @@ -4208,11 +4281,10 @@ def _create_mdc_store(cls, model, dataset, circuits, resource_alloc, return ModelDatasetCircuitsStore(model, dataset, circuits, resource_alloc, array_types, None, verbosity) @classmethod - def create_from(cls, raw_objfn, model, dataset, circuits, resource_alloc=None, penalties=None, - verbosity=0, method_names=('fn',), array_types=()): - mdc_store = cls._create_mdc_store(model, dataset, circuits, resource_alloc, method_names, - array_types, verbosity) - return cls(raw_objfn, mdc_store, penalties, verbosity) + def create_from(cls, model, dataset, circuits, regularization=None, penalties=None, resource_alloc=None, + name=None, description=None, verbosity=0, method_names=('fn',), array_types=(), **kwargs): + mdc_store = cls._create_mdc_store(model, dataset, circuits, resource_alloc, method_names, array_types, verbosity) + return cls(mdc_store, regularization, penalties, name, description, verbosity, **kwargs) @classmethod def _array_types_for_method(cls, method_name, fsim): @@ -4243,9 +4315,10 @@ def compute_array_types(cls, method_names, fsim): return array_types - def __init__(self, raw_objfn, mdc_store, penalties=None, verbosity=0): + @set_docstring(DOCSTR_TEMPLATE % TEMPLATE_FIELDS) + def __init__(self, raw_objfn, mdc_store, penalties=None, verbosity=0, **kwargs): - super().__init__(raw_objfn, mdc_store, verbosity) + super().__init__(raw_objfn, mdc_store, verbosity, **kwargs) if penalties is None: penalties = {} self.ex = self.set_penalties(**penalties) @@ -4363,7 +4436,11 @@ def set_penalties(self, regularize_factor=0, cptp_penalty_factor=0, spam_penalty if self.cptp_penalty_factor != 0: ex += _cptp_penalty_size(self.model) if self.spam_penalty_factor != 0: ex += _spam_penalty_size(self.model) if self.errorgen_penalty_factor != 0: ex += _errorgen_penalty_size(self.model) - + if self.callable_penalty is not None: ex += 1 + # ^ We want to be able to change the penalty factor without having to rebuild + # this objective function. So we allocate space for the callable penalty + # if it exists at all, rather than requiring self.callable_penalty_factor > 0. + # return ex def terms(self, paramvec=None, oob_check=False, profiler_str="TERMS OBJECTIVE"): @@ -4515,6 +4592,31 @@ def _dterms_fill_penalty(self, paramvec, terms_jac): terms_jac[off:off+n, :] *= 2*errorgen_penalty[:, None] off += n + if self.callable_penalty is not None: + # + # Ideally we'd have put in the effort to do this earlier when doing finite-difference to get + # the Jacobian of terms. But we start with a simple implementation. + # + # --> TODO: change so that if callable_penalty has a grad method then we'll call it instead + # of relying on finite-differences over all model parameters. This would help when + # the penalty is constant w.r.t. some params (like those that only appear in SPAM). + # + terms_jac[off:off+1, :] = 0.0 + if self.callable_penalty_factor: + vec0 = self.model.to_vector() + val0 = self.callable_penalty(self.model) + eps = 1e-7 + for i in range(vec0.size): + vec0[i] += eps + self.model.from_vector(vec0) + val = self.callable_penalty(self.model) + dval = (val - val0) / eps + terms_jac[off, i] = dval + vec0[i] -= eps + self.model.from_vector(vec0) + off += 1 + pass + assert(off == self.local_ex) return @@ -4543,6 +4645,11 @@ def _terms_penalty(self, paramvec): errorgen_penalty = _errorgen_penalty(self.model, self.errorgen_penalty_factor) ** 2 blocks.append(errorgen_penalty) + if self.callable_penalty is not None: + f = self.callable_penalty_factor + val = f * self.callable_penalty(self.model) if f > 0 else 0.0 + blocks.append(_np.array([val])) + return _np.concatenate(blocks) def _clip_probs(self): @@ -4752,231 +4859,82 @@ def _hessian_from_block(self, hprobs, dprobs12, probs, element_slice, counts, to class Chi2Function(TimeIndependentMDCObjectiveFunction): - """ - Model-based chi-squared function: `N(p-f)^2 / p` - - Parameters - ---------- - mdl : Model - The model - specifies how parameter values are turned into probabilities - for each circuit outcome. - - dataset : DataSet - The data set - specifies how counts and total_counts are obtained for each - circuit outcome. - - circuits : list or CircuitList - The circuit list - specifies what probabilities and counts this objective - function compares. If `None`, then the keys of `dataset` are used. - - regularization : dict, optional - Regularization values. - - penalties : dict, optional - Penalty values. Penalties usually add additional (penalty) terms to the sum - of per-circuit-outcome contributions that evaluate to the objective function. - - resource_alloc : ResourceAllocation, optional - Available resources and how they should be allocated for computations. - - name : str, optional - A name for this objective function (can be anything). - - description : str, optional - A description for this objective function (can be anything) - - verbosity : int, optional - Level of detail to print to stdout. - enable_hessian : bool, optional - Whether hessian calculations are allowed. If `True` then more resources are - needed. If `False`, calls to hessian-requiring function will result in an - error. + TEMPLATE_FIELDS = ( """ - @classmethod - def create_from(cls, model, dataset, circuits, regularization=None, penalties=None, resource_alloc=None, - name=None, description=None, verbosity=0, method_names=('fn',), array_types=()): - mdc_store = cls._create_mdc_store(model, dataset, circuits, resource_alloc, method_names, - array_types, verbosity) - return cls(mdc_store, regularization, penalties, name, description, verbosity) + Model-based chi-squared function: `N(p-f)^2 / p` + """, "", "" + ) - def __init__(self, mdc_store, regularization=None, penalties=None, name=None, description=None, verbosity=0): + @set_docstring(TimeIndependentMDCObjectiveFunction.DOCSTR_TEMPLATE % TEMPLATE_FIELDS) + def __init__(self, mdc_store, regularization=None, penalties=None, name=None, description=None, verbosity=0, **kwargs): raw_objfn = RawChi2Function(regularization, mdc_store.resource_alloc, name, description, verbosity) - super().__init__(raw_objfn, mdc_store, penalties, verbosity) + super().__init__(raw_objfn, mdc_store, penalties, verbosity, **kwargs) class ChiAlphaFunction(TimeIndependentMDCObjectiveFunction): + + TEMPLATE_FIELDS = ( """ Model-based chi-alpha function: `N[x + 1/(alpha * x^alpha) - (1 + 1/alpha)]` where `x := p/f`. - - Parameters - ---------- - mdl : Model - The model - specifies how parameter values are turned into probabilities - for each circuit outcome. - - dataset : DataSet - The data set - specifies how counts and total_counts are obtained for each - circuit outcome. - - circuits : list or CircuitList - The circuit list - specifies what probabilities and counts this objective - function compares. If `None`, then the keys of `dataset` are used. - - regularization : dict, optional - Regularization values. - - penalties : dict, optional - Penalty values. Penalties usually add additional (penalty) terms to the sum - of per-circuit-outcome contributions that evaluate to the objective function. - - resource_alloc : ResourceAllocation, optional - Available resources and how they should be allocated for computations. - - name : str, optional - A name for this objective function (can be anything). - - description : str, optional - A description for this objective function (can be anything) - - verbosity : int, optional - Level of detail to print to stdout. - - enable_hessian : bool, optional - Whether hessian calculations are allowed. If `True` then more resources are - needed. If `False`, calls to hessian-requiring function will result in an - error. - + """, + "", # no custom leading arguments. + """ alpha : float, optional The alpha parameter, which lies in the interval (0,1]. - """ + """ # one custom trailing argument. + ) @classmethod def create_from(cls, model, dataset, circuits, regularization=None, penalties=None, resource_alloc=None, name=None, description=None, verbosity=0, - method_names=('fn',), array_types=(), alpha=1): - mdc_store = cls._create_mdc_store(model, dataset, circuits, resource_alloc, method_names, - array_types, verbosity) - return cls(mdc_store, regularization, penalties, name, description, verbosity, alpha) + method_names=('fn',), array_types=(), alpha=1, **kwargs): + mdc_store = cls._create_mdc_store(model, dataset, circuits, resource_alloc, method_names, array_types, verbosity) + return cls(mdc_store, regularization, penalties, name, description, verbosity, alpha, **kwargs) + @set_docstring(TimeIndependentMDCObjectiveFunction.DOCSTR_TEMPLATE % TEMPLATE_FIELDS) def __init__(self, mdc_store, regularization=None, penalties=None, name=None, description=None, verbosity=0, - alpha=1): + alpha=1, **kwargs): raw_objfn = RawChiAlphaFunction(regularization, mdc_store.resource_alloc, name, description, verbosity, alpha) - super().__init__(raw_objfn, mdc_store, penalties, verbosity) + super().__init__(raw_objfn, mdc_store, penalties, verbosity, **kwargs) class FreqWeightedChi2Function(TimeIndependentMDCObjectiveFunction): - """ - Model-based frequency-weighted chi-squared function: `N(p-f)^2 / f` - - Parameters - ---------- - mdl : Model - The model - specifies how parameter values are turned into probabilities - for each circuit outcome. - - dataset : DataSet - The data set - specifies how counts and total_counts are obtained for each - circuit outcome. - - circuits : list or CircuitList - The circuit list - specifies what probabilities and counts this objective - function compares. If `None`, then the keys of `dataset` are used. - - regularization : dict, optional - Regularization values. - - penalties : dict, optional - Penalty values. Penalties usually add additional (penalty) terms to the sum - of per-circuit-outcome contributions that evaluate to the objective function. - - resource_alloc : ResourceAllocation, optional - Available resources and how they should be allocated for computations. - - name : str, optional - A name for this objective function (can be anything). - - description : str, optional - A description for this objective function (can be anything) - - verbosity : int, optional - Level of detail to print to stdout. - enable_hessian : bool, optional - Whether hessian calculations are allowed. If `True` then more resources are - needed. If `False`, calls to hessian-requiring function will result in an - error. + TEMPLATE_FIELDS = ( """ + Model-based frequency-weighted chi-squared function: `N(p-f)^2 / f` + """, "", "" + ) - @classmethod - def create_from(cls, model, dataset, circuits, regularization=None, penalties=None, - resource_alloc=None, name=None, description=None, verbosity=0, - method_names=('fn',), array_types=()): - mdc_store = cls._create_mdc_store(model, dataset, circuits, resource_alloc, method_names, - array_types, verbosity) - return cls(mdc_store, regularization, penalties, name, description, verbosity) - - def __init__(self, mdc_store, regularization=None, penalties=None, name=None, description=None, verbosity=0): + @set_docstring(TimeIndependentMDCObjectiveFunction.DOCSTR_TEMPLATE % TEMPLATE_FIELDS) + def __init__(self, mdc_store, regularization=None, penalties=None, name=None, description=None, verbosity=0, **kwargs): raw_objfn = RawFreqWeightedChi2Function(regularization, mdc_store.resource_alloc, name, description, verbosity) - super().__init__(raw_objfn, mdc_store, penalties, verbosity) + super().__init__(raw_objfn, mdc_store, penalties, verbosity, **kwargs) class CustomWeightedChi2Function(TimeIndependentMDCObjectiveFunction): + + TEMPLATE_FIELDS = ( """ Model-based custom-weighted chi-squared function: `cw^2 (p-f)^2` - - Parameters - ---------- - mdl : Model - The model - specifies how parameter values are turned into probabilities - for each circuit outcome. - - dataset : DataSet - The data set - specifies how counts and total_counts are obtained for each - circuit outcome. - - circuits : list or CircuitList - The circuit list - specifies what probabilities and counts this objective - function compares. If `None`, then the keys of `dataset` are used. - - regularization : dict, optional - Regularization values. - - penalties : dict, optional - Penalty values. Penalties usually add additional (penalty) terms to the sum - of per-circuit-outcome contributions that evaluate to the objective function. - - resource_alloc : ResourceAllocation, optional - Available resources and how they should be allocated for computations. - - name : str, optional - A name for this objective function (can be anything). - - description : str, optional - A description for this objective function (can be anything) - - verbosity : int, optional - Level of detail to print to stdout. - - enable_hessian : bool, optional - Whether hessian calculations are allowed. If `True` then more resources are - needed. If `False`, calls to hessian-requiring function will result in an - error. - + """, "", + """ custom_weights : numpy.ndarray, optional One-dimensional array of the custom weights, which linearly multiply the *least-squares* terms, i.e. `(p - f)`. If `None`, then unit weights are used and the objective function computes the sum of unweighted squares. """ + ) @classmethod def create_from(cls, model, dataset, circuits, regularization=None, penalties=None, resource_alloc=None, name=None, description=None, verbosity=0, - method_names=('fn',), array_types=(), custom_weights=None): - mdc_store = cls._create_mdc_store(model, dataset, circuits, resource_alloc, method_names, - array_types, verbosity) - return cls(mdc_store, regularization, penalties, name, description, verbosity, custom_weights) + method_names=('fn',), array_types=(), custom_weights=None, **kwargs): + mdc_store = cls._create_mdc_store(model, dataset, circuits, resource_alloc, method_names, array_types, verbosity) + return cls(mdc_store, regularization, penalties, name, description, verbosity, custom_weights, **kwargs) + @set_docstring(TimeIndependentMDCObjectiveFunction.DOCSTR_TEMPLATE % TEMPLATE_FIELDS) def __init__(self, mdc_store, regularization=None, penalties=None, name=None, description=None, verbosity=0, custom_weights=None): raw_objfn = RawCustomWeightedChi2Function(regularization, mdc_store.resource_alloc, name, description, @@ -4985,230 +4943,215 @@ def __init__(self, mdc_store, regularization=None, penalties=None, name=None, de class PoissonPicDeltaLogLFunction(TimeIndependentMDCObjectiveFunction): - """ - Model-based poisson-picture delta log-likelihood function: `N*f*log(f/p) - N*(f-p)`. - - Parameters - ---------- - mdl : Model - The model - specifies how parameter values are turned into probabilities - for each circuit outcome. - dataset : DataSet - The data set - specifies how counts and total_counts are obtained for each - circuit outcome. - - circuits : list or CircuitList - The circuit list - specifies what probabilities and counts this objective - function compares. If `None`, then the keys of `dataset` are used. - - regularization : dict, optional - Regularization values. - - penalties : dict, optional - Penalty values. Penalties usually add additional (penalty) terms to the sum - of per-circuit-outcome contributions that evaluate to the objective function. - - resource_alloc : ResourceAllocation, optional - Available resources and how they should be allocated for computations. - - name : str, optional - A name for this objective function (can be anything). - - description : str, optional - A description for this objective function (can be anything) - - verbosity : int, optional - Level of detail to print to stdout. - - enable_hessian : bool, optional - Whether hessian calculations are allowed. If `True` then more resources are - needed. If `False`, calls to hessian-requiring function will result in an - error. + TEMPLATE_FIELDS = ( """ + Model-based poisson-picture delta log-likelihood function: `N*f*log(f/p) - N*(f-p)`. + """, "", "" + ) - @classmethod - def create_from(cls, model, dataset, circuits, regularization=None, penalties=None, - resource_alloc=None, name=None, description=None, verbosity=0, - method_names=('fn',), array_types=()): - mdc_store = cls._create_mdc_store(model, dataset, circuits, resource_alloc, method_names, - array_types, verbosity) - return cls(mdc_store, regularization, penalties, name, description, verbosity) - - def __init__(self, mdc_store, regularization=None, penalties=None, name=None, description=None, verbosity=0): + @set_docstring(TimeIndependentMDCObjectiveFunction.DOCSTR_TEMPLATE % TEMPLATE_FIELDS) + def __init__(self, mdc_store, regularization=None, penalties=None, name=None, description=None, verbosity=0, **kwargs): raw_objfn = RawPoissonPicDeltaLogLFunction(regularization, mdc_store.resource_alloc, name, description, verbosity) - super().__init__(raw_objfn, mdc_store, penalties, verbosity) + super().__init__(raw_objfn, mdc_store, penalties, verbosity, **kwargs) class DeltaLogLFunction(TimeIndependentMDCObjectiveFunction): + + TEMPLATE_FIELDS = ( """ Model-based delta log-likelihood function: `N*f*log(f/p)`. + """, "", "" + ) - Parameters - ---------- - mdl : Model - The model - specifies how parameter values are turned into probabilities - for each circuit outcome. - - dataset : DataSet - The data set - specifies how counts and total_counts are obtained for each - circuit outcome. - - circuits : list or CircuitList - The circuit list - specifies what probabilities and counts this objective - function compares. If `None`, then the keys of `dataset` are used. - - regularization : dict, optional - Regularization values. - - penalties : dict, optional - Penalty values. Penalties usually add additional (penalty) terms to the sum - of per-circuit-outcome contributions that evaluate to the objective function. - - resource_alloc : ResourceAllocation, optional - Available resources and how they should be allocated for computations. - - name : str, optional - A name for this objective function (can be anything). + @set_docstring(TimeIndependentMDCObjectiveFunction.DOCSTR_TEMPLATE % TEMPLATE_FIELDS) + def __init__(self, mdc_store, regularization=None, penalties=None, name=None, description=None, verbosity=0, **kwargs): + raw_objfn = RawDeltaLogLFunction(regularization, mdc_store.resource_alloc, name, description, verbosity) + super().__init__(raw_objfn, mdc_store, penalties, verbosity, **kwargs) - description : str, optional - A description for this objective function (can be anything) - verbosity : int, optional - Level of detail to print to stdout. +class MaxLogLFunction(TimeIndependentMDCObjectiveFunction): - enable_hessian : bool, optional - Whether hessian calculations are allowed. If `True` then more resources are - needed. If `False`, calls to hessian-requiring function will result in an - error. + TEMPLATE_FIELDS = ( """ + Model-based maximum-model log-likelihood function: `N*f*log(f)` + """, "", "" + ) @classmethod def create_from(cls, model, dataset, circuits, regularization=None, penalties=None, resource_alloc=None, name=None, description=None, verbosity=0, - method_names=('fn',), array_types=()): - mdc_store = cls._create_mdc_store(model, dataset, circuits, resource_alloc, method_names, - array_types, verbosity) - return cls(mdc_store, regularization, penalties, name, description, verbosity) + method_names=('fn',), array_types=(), poisson_picture=True, **kwargs): + mdc_store = cls._create_mdc_store(model, dataset, circuits, resource_alloc, method_names, array_types, verbosity) + return cls(mdc_store, regularization, penalties, name, description, verbosity, poisson_picture, **kwargs) - def __init__(self, mdc_store, regularization=None, penalties=None, name=None, description=None, verbosity=0): - raw_objfn = RawDeltaLogLFunction(regularization, mdc_store.resource_alloc, name, description, verbosity) - super().__init__(raw_objfn, mdc_store, penalties, verbosity) + @set_docstring(TimeIndependentMDCObjectiveFunction.DOCSTR_TEMPLATE % TEMPLATE_FIELDS) + def __init__(self, mdc_store, regularization=None, penalties=None, name=None, description=None, verbosity=0, + poisson_picture=True, **kwargs): + raw_objfn = RawMaxLogLFunction(regularization, mdc_store.resource_alloc, name, description, verbosity, + poisson_picture) + super().__init__(raw_objfn, mdc_store, penalties, verbosity, **kwargs) -class MaxLogLFunction(TimeIndependentMDCObjectiveFunction): +class TermWeighted(TimeIndependentMDCObjectiveFunction): """ - Model-based maximum-model log-likelihood function: `N*f*log(f)` - - Parameters - ---------- - mdl : Model - The model - specifies how parameter values are turned into probabilities - for each circuit outcome. - - dataset : DataSet - The data set - specifies how counts and total_counts are obtained for each - circuit outcome. + Represents an objective whose term function takes the form + f(params) = w * g(params), + where w is constant and g(params) is a black-box term function. - circuits : list or CircuitList - The circuit list - specifies what probabilities and counts this objective - function compares. If `None`, then the keys of `dataset` are used. + This class' implementation of terms() and dterms() has a substantial amount of code + duplication with the implementations in TimeIndependentMDCObjectiveFunction + """ - regularization : dict, optional - Regularization values. + @property + def terms_weights(self) -> _np.ndarray: + if not hasattr(self, '_terms_weights'): + self._update_terms_weights() + return self._terms_weights - penalties : dict, optional - Penalty values. Penalties usually add additional (penalty) terms to the sum - of per-circuit-outcome contributions that evaluate to the objective function. + def _update_terms_weights(self): + self._terms_weights = _np.ones(self.layout.num_elements) + return - resource_alloc : ResourceAllocation, optional - Available resources and how they should be allocated for computations. + def terms(self, paramvec=None, oob_check=False, profiler_str="TERMS OBJECTIVE"): + tm = _time.time() + if paramvec is None: + paramvec = self.model.to_vector() + else: + self.model.from_vector(paramvec) + terms = self.obj.view() - name : str, optional - A name for this objective function (can be anything). + unit_ralloc = self.layout.resource_alloc('atom-processing') + shared_mem_leader = unit_ralloc.is_host_leader - description : str, optional - A description for this objective function (can be anything) + with self.resource_alloc.temporarily_track_memory(self.nelements): # 'e' (terms) + self.model.sim.bulk_fill_probs(self.probs, self.layout) + self._clip_probs() + + if oob_check: # Only used for termgap cases + if not self.model.sim.bulk_test_if_paths_are_sufficient(self.layout, self.probs, verbosity=1): + raise ValueError("Out of bounds!") # signals LM optimizer - verbosity : int, optional - Level of detail to print to stdout. + if shared_mem_leader: + terms_no_penalty = self.raw_objfn.terms(self.probs, self.counts, self.total_counts, self.freqs) + terms[:self.nelements] = terms_no_penalty + if self._process_penalties: + terms[self.nelements:] = self._terms_penalty(paramvec) - enable_hessian : bool, optional - Whether hessian calculations are allowed. If `True` then more resources are - needed. If `False`, calls to hessian-requiring function will result in an - error. - """ + if self.firsts is not None and shared_mem_leader: + omitted_probs = 1.0 - _np.array([self.probs[self.layout.indices_for_index(i)].sum() for i in self.indicesOfCircuitsWithOmittedData]) + omitted_probs_firsts_terms = self.raw_objfn.zero_freq_terms(self.total_counts[self.firsts], omitted_probs) + terms[self.firsts] += omitted_probs_firsts_terms + + self._update_terms_weights() + terms[:self.nelements] *= self.terms_weights + + unit_ralloc.host_comm_barrier() - @classmethod - def create_from(cls, model, dataset, circuits, regularization=None, penalties=None, - resource_alloc=None, name=None, description=None, verbosity=0, - method_names=('fn',), array_types=(), poisson_picture=True): - mdc_store = cls._create_mdc_store(model, dataset, circuits, resource_alloc, method_names, - array_types, verbosity) - return cls(mdc_store, regularization, penalties, name, description, verbosity, poisson_picture) + self.raw_objfn.resource_alloc.profiler.add_time(profiler_str, tm) + assert(terms.shape == (self.nelements + self.local_ex,)) + terms *= self.terms_weights + return terms - def __init__(self, mdc_store, regularization=None, penalties=None, name=None, description=None, verbosity=0, - poisson_picture=True): - raw_objfn = RawMaxLogLFunction(regularization, mdc_store.resource_alloc, name, description, verbosity, - poisson_picture) - super().__init__(raw_objfn, mdc_store, penalties, verbosity) + def dterms(self, paramvec=None): + tm = _time.time() + unit_ralloc = self.layout.resource_alloc('param-processing') + shared_mem_leader = unit_ralloc.is_host_leader + if paramvec is None: + paramvec = self.model.to_vector() + else: + self.model.from_vector(paramvec) + + dprobs = self.jac[0:self.nelements, :] + dprobs.shape = (self.nelements, self.nparams) -class TVDFunction(TimeIndependentMDCObjectiveFunction): - """ - Model-based TVD function: `0.5 * |p-f|`. + with self.resource_alloc.temporarily_track_memory(2 * self.nelements): + self.model.sim.bulk_fill_dprobs(dprobs, self.layout, self.probs) + self._clip_probs() + if shared_mem_leader: + if self.firsts is not None: + for ii, i in enumerate(self.indicesOfCircuitsWithOmittedData): + self.dprobs_omitted_rowsum[ii, :] = _np.sum(dprobs[self.layout.indices_for_index(i), :], axis=0) + dg_probs = self.raw_objfn.dterms(self.probs, self.counts, self.total_counts, self.freqs) + dprobs *= dg_probs[:, None] + + if shared_mem_leader and self.firsts is not None: + total_counts_firsts = self.total_counts[self.firsts] + omitted_probs = 1.0 - _np.array([_np.sum(self.probs[self.layout.indices_for_index(i)]) for i in self.indicesOfCircuitsWithOmittedData]) + omitted_dprobs_firsts_dterms = self.raw_objfn.zero_freq_dterms(total_counts_firsts, omitted_probs) + dprobs[self.firsts] -= omitted_dprobs_firsts_dterms[:, None] * self.dprobs_omitted_rowsum + + self._update_terms_weights() + dprobs[:self.nelements] *= self.terms_weights[:, None] - Parameters - ---------- - mdl : Model - The model - specifies how parameter values are turned into probabilities - for each circuit outcome. + if shared_mem_leader and self._process_penalties: + self._dterms_fill_penalty(paramvec, self.jac[self.nelements:, :]) - dataset : DataSet - The data set - specifies how counts and total_counts are obtained for each - circuit outcome. + unit_ralloc.host_comm_barrier() + self.raw_objfn.resource_alloc.profiler.add_time("JACOBIAN", tm) + return self.jac - circuits : list or CircuitList - The circuit list - specifies what probabilities and counts this objective - function compares. If `None`, then the keys of `dataset` are used. - regularization : dict, optional - Regularization values. +class TVDFunction(TermWeighted): - penalties : dict, optional - Penalty values. Penalties usually add additional (penalty) terms to the sum - of per-circuit-outcome contributions that evaluate to the objective function. + TEMPLATE_FIELDS = ( + """ + Model-based TVD function: `0.5 * w * |p-f|`, where w is a vector of weights. - resource_alloc : ResourceAllocation, optional - Available resources and how they should be allocated for computations. + If `name == 'normalized tvd'`, then `w[i]` will be 1/(length of circuit associated with i). + Otherwise, w[i] will equal 1. + """, "", "" + ) - name : str, optional - A name for this objective function (can be anything). + @set_docstring(TermWeighted.DOCSTR_TEMPLATE % TEMPLATE_FIELDS) + def __init__(self, mdc_store, regularization=None, penalties=None, name=None, description=None, verbosity=0, **kwargs): + raw_objfn = RawTVDFunction(regularization, mdc_store.resource_alloc, name, description, verbosity) + super().__init__(raw_objfn, mdc_store, penalties, verbosity, **kwargs) + self.normalized = name == 'normalized tvd' + self._terms_weights = None + self._update_terms_weights() + + def _update_terms_weights(self): + if self.normalized: + num_elements = self.layout.num_elements + circuit_sizes = _np.zeros((num_elements,)) + for elind, circuit, _ in self.layout: + circuit_sizes[elind] = 1 + circuit.num_gates + assert _np.all(circuit_sizes > 0) + self._terms_weights = 1 / circuit_sizes + self._terms_weights /= _np.mean(self.terms_weights) + # ^ Make the weights mean-1 to avoid unintended changes to Levenberg-Marquardt + # stopping criteria. This has an undesirable side effect of affecting how + # one would interpret individual circuits' contributions to the terms vector. + # However, this function is first and foremost a sum of terms, and the sum + # never has an interesting interpretation. + # + elif self._terms_weights is None: + self._terms_weights = _np.ones(self.layout.num_elements) + return - description : str, optional - A description for this objective function (can be anything) - verbosity : int, optional - Level of detail to print to stdout. +class LpNormToPowerP(TermWeighted): - enable_hessian : bool, optional - Whether hessian calculations are allowed. If `True` then more resources are - needed. If `False`, calls to hessian-requiring function will result in an - error. + TEMPLATE_FIELDS = ( """ + Model-based loss function: `0.5 * |p-f|^power`. + """, "", + """ + power : float, optonal + Must be >= 1. + """ + ) - @classmethod - def create_from(cls, model, dataset, circuits, regularization=None, penalties=None, - resource_alloc=None, name=None, description=None, verbosity=0, - method_names=('fn',), array_types=()): - mdc_store = cls._create_mdc_store(model, dataset, circuits, resource_alloc, method_names, - array_types, verbosity) - return cls(mdc_store, regularization, penalties, name, description, verbosity) - - def __init__(self, mdc_store, regularization=None, penalties=None, name=None, description=None, verbosity=0): - raw_objfn = RawTVDFunction(regularization, mdc_store.resource_alloc, name, description, verbosity) - super().__init__(raw_objfn, mdc_store, penalties, verbosity) + @set_docstring(TermWeighted.DOCSTR_TEMPLATE % TEMPLATE_FIELDS) + def __init__(self, mdc_store, regularization=None, penalties=None, name=None, description=None, + verbosity=0, power=2, **kwargs): + raw_objfn = RawAbsPower(power, regularization, mdc_store.resource_alloc, name, description, verbosity) + super().__init__(raw_objfn, mdc_store, penalties, verbosity, **kwargs) + self.power = raw_objfn.power + self._update_terms_weights() class TimeDependentMDCObjectiveFunction(MDCObjectiveFunction): @@ -5246,31 +5189,6 @@ class TimeDependentMDCObjectiveFunction(MDCObjectiveFunction): Level of detail to print to stdout. """ - @classmethod - def builder(cls, name=None, description=None, regularization=None, penalties=None, **kwargs): - """ - Create an :class:`ObjectiveFunctionBuilder` that builds an objective function of this type. - - Parameters - ---------- - name : str, optional - A name for the built objective function (can be anything). - - description : str, optional - A description for the built objective function (can be anything) - - regularization : dict, optional - Regularization values. - - penalties : dict, optional - Penalty values. - - Returns - ------- - ObjectiveFunctionBuilder - """ - return ObjectiveFunctionBuilder(cls, name, description, regularization, penalties, **kwargs) - #This objective function can handle time-dependent circuits - that is, circuits are treated as # potentially time-dependent and mdl as well. For now, we don't allow any regularization or penalization # in this case. @@ -5278,7 +5196,7 @@ def builder(cls, name=None, description=None, regularization=None, penalties=Non @classmethod def create_from(cls, model, dataset, circuits, regularization=None, penalties=None, resource_alloc=None, name=None, description=None, verbosity=0, - method_names=('fn',), array_types=()): + method_names=('fn',), array_types=(), **kwargs): #Array types are used to construct memory estimates (as a function of element number, etc) for layout creation. # They account for memory used in: # 1) an optimization method (if present), @@ -5286,7 +5204,7 @@ def create_from(cls, model, dataset, circuits, regularization=None, penalties=No # 2b) intermediate memory allocated by methods of the created object (possibly an objective function) array_types += cls.compute_array_types(method_names, model.sim) mdc_store = ModelDatasetCircuitsStore(model, dataset, circuits, resource_alloc, array_types) - return cls(mdc_store, regularization, penalties, name, description, verbosity) + return cls(mdc_store, regularization, penalties, name, description, verbosity, **kwargs) @classmethod def compute_array_types(cls, method_names, fsim): diff --git a/pygsti/optimize/simplerlm.py b/pygsti/optimize/simplerlm.py index a027bfd53..9b7512b0f 100644 --- a/pygsti/optimize/simplerlm.py +++ b/pygsti/optimize/simplerlm.py @@ -287,6 +287,8 @@ def run(self, objective: TimeIndependentMDCObjectiveFunction, profiler, printer) objective_func, jacobian, x0, max_iter=self.maxiter, num_fd_iters=self.fditer, + # NOTE: I think the fallback values below will NEVER be triggered. + # Should probably remove them. See __init__ instead. f_norm2_tol=self.tol.get('f', 1.0), jac_norm_tol=self.tol.get('jac', 1e-6), rel_ftol=self.tol.get('relf', 1e-6), @@ -583,7 +585,7 @@ def simplish_leastsq( if norm_JTf < jac_norm_tol: if oob_check_interval <= 1: - msg = "norm(jacobian) is at most %g" % jac_norm_tol + msg = "norm(J'f) is at most %g" % jac_norm_tol converged = True break else: diff --git a/pygsti/protocols/estimate.py b/pygsti/protocols/estimate.py index d759dd758..2b81b1669 100644 --- a/pygsti/protocols/estimate.py +++ b/pygsti/protocols/estimate.py @@ -163,7 +163,7 @@ def __init__(self, parent, models=None, parameters=None, extra_parameters=None): self.profiler = parameters.get('profiler', None) self._final_mdc_store = parameters.get('final_mdc_store', None) self._final_objfn_cache = parameters.get('final_objfn_cache', None) - self.final_objfn_builder = parameters.get('final_objfn_builder', _objfns.PoissonPicDeltaLogLFunction.builder()) + self.final_objfn_builder = parameters.get('final_objfn_builder', _objfns.ObjectiveFunctionBuilder(_objfns.PoissonPicDeltaLogLFunction)) self._final_objfn = parameters.get('final_objfn', None) self._per_iter_mdc_store = parameters.get('per_iter_mdc_store', None) diff --git a/pygsti/protocols/gst.py b/pygsti/protocols/gst.py index 83680912f..c6fd21c56 100644 --- a/pygsti/protocols/gst.py +++ b/pygsti/protocols/gst.py @@ -735,8 +735,9 @@ class GSTObjFnBuilders(_NicelySerializable): on the final GST iteration. """ - @classmethod - def cast(cls, obj): + # This used to be a class method, but this class has no derived classes. + @staticmethod + def cast(obj): """ Cast `obj` to a :class:`GSTObjFnBuilders` object. @@ -751,20 +752,22 @@ def cast(cls, obj): ------- GSTObjFnBuilders """ + cls = GSTObjFnBuilders if isinstance(obj, cls): return obj elif obj is None: return cls.create_from() elif isinstance(obj, dict): return cls.create_from(**obj) elif isinstance(obj, (list, tuple)): return cls(*obj) else: raise ValueError("Cannot create an %s object from '%s'" % (cls.__name__, str(type(obj)))) - @classmethod - def create_from(cls, objective='logl', freq_weighted_chi2=False, always_perform_mle=False, only_perform_mle=False): + # This used to be a class method, but this class has no derived classes. + @staticmethod + def create_from(objective='logl', freq_weighted_chi2=False, always_perform_mle=False, only_perform_mle=False): """ Creates a common :class:`GSTObjFnBuilders` object from several arguments. Parameters ---------- - objective : {'logl', 'chi2'}, optional + objective : {'logl', 'chi2', 'tvd'}, optional Whether to create builders for maximum-likelihood or minimum-chi-squared GST. freq_weighted_chi2 : bool, optional @@ -786,10 +789,15 @@ def create_from(cls, objective='logl', freq_weighted_chi2=False, always_perform_ chi2_builder = _objfns.ObjectiveFunctionBuilder.create_from('chi2', freq_weighted_chi2) mle_builder = _objfns.ObjectiveFunctionBuilder.create_from('logl') - if objective == "chi2": + if not isinstance(objective, str): + import warnings + warnings.warn(f'Trying to create an objective function from non-string specification, "{objective}". \ + \nSupport for this kind of specification is experimental!') + iteration_builders = [chi2_builder] + final_builders = [_objfns.ObjectiveFunctionBuilder.create_from(objective)] + elif objective == "chi2": iteration_builders = [chi2_builder] final_builders = [] - elif objective == "logl": if always_perform_mle: iteration_builders = [mle_builder] if only_perform_mle else [chi2_builder, mle_builder] @@ -798,8 +806,13 @@ def create_from(cls, objective='logl', freq_weighted_chi2=False, always_perform_ iteration_builders = [chi2_builder] final_builders = [mle_builder] else: - raise ValueError("Invalid objective: %s" % objective) - return cls(iteration_builders, final_builders) + iteration_builders = [chi2_builder] + try: + final_builders = [_objfns.ObjectiveFunctionBuilder.create_from(objective)] + except Exception as e: + raise ValueError("Invalid objective: %s" % objective) + + return GSTObjFnBuilders(iteration_builders, final_builders) def __init__(self, iteration_builders, final_builders=()): super().__init__() @@ -1662,7 +1675,7 @@ def run(self, data, memlimit=None, comm=None): parameters['protocol'] = self # Estimates can hold sub-Protocols <=> sub-results parameters['profiler'] = profiler parameters['final_mdc_store'] = final_store - parameters['final_objfn_builder'] = _objfns.PoissonPicDeltaLogLFunction.builder() + parameters['final_objfn_builder'] = _objfns.ObjectiveFunctionBuilder(_objfns.PoissonPicDeltaLogLFunction) # just set final objective function as default logl objective (for ease of later comparison) ret = ModelEstimateResults(data, self) @@ -3136,6 +3149,9 @@ def add_estimate(self, estimate, estimate_key='default'): _warnings.warn("Re-initializing the %s estimate" % estimate_key + " of this Results object! Usually you don't" + " want to do this.") + + if estimate.parent is None: + estimate.set_parent(self) self.estimates[estimate_key] = estimate diff --git a/pygsti/protocols/locking.py b/pygsti/protocols/locking.py new file mode 100644 index 000000000..9f9bf764e --- /dev/null +++ b/pygsti/protocols/locking.py @@ -0,0 +1,130 @@ +# *************************************************************************************************** +# Copyright 2025 National Technology & Engineering Solutions of Sandia, LLC (NTESS). +# Under the terms of Contract DE-NA0003525 with NTESS, the U.S. Government retains certain rights +# in this software. +# Licensed under the Apache License, Version 2.0 (the "License"); you may not use this file except +# in compliance with the License. You may obtain a copy of the License at +# http://www.apache.org/licenses/LICENSE-2.0 or in the LICENSE file in the root pyGSTi directory. +# *************************************************************************************************** +import copy as _copy +import numpy as _np +import itertools as _itertools +import pathlib as _pathlib +import warnings as _warnings + +from pygsti.protocols.treenode import TreeNode as _TreeNode +from pygsti import io as _io +from collections.abc import Container as _Container +from pygsti.circuits import ( + Circuit as _Circuit, + CircuitList as _CircuitList, + CircuitPlaquette as _CircuitPlaquette +) +from pygsti import data as _data +from pygsti.tools import NamedDict as _NamedDict +from pygsti.tools import listtools as _lt +from pygsti.tools.dataframetools import _process_dataframe +from pygsti.baseobjs.label import ( + LabelStr as _LabelStr +) +from pygsti.baseobjs.mongoserializable import MongoSerializable as _MongoSerializable +from pygsti.baseobjs.nicelyserializable import NicelySerializable as _NicelySerializable +from pygsti.protocols.protocol import Protocol as _Protocol, CircuitListsDesign as _CircuitListsDesign + +from typing import Union, Literal + +BinningStrategy = Union[ + int, Literal['auto-int'], Literal['auto'], Literal['fd'], Literal['doane'], Literal['scott'], Literal['stone'], Literal['sturges'] +] +LengthTransformer = Union[ + Literal['log'], Union[Literal['none'], None], _np.ufunc +] + +def histonested_circuitlists( + circuits: Union[_CircuitList, list[_Circuit]], + bins: BinningStrategy = 'auto-int', + trans: LengthTransformer = 'log' + ) -> list[_Circuit]: + """ + This is a helper function for building CircuitListsDesign objects with certain + nested structures. If `clists` is the output of this function, then the induced + design is canonically + + d = CircuitListsDesign(clists nested=True, remove_duplicates=True). + + If `circuits` contained no duplicates, then we'll have + + set(circuits) == set(d.all_circuits_needing_data). + """ + assert len(circuits) > 0 + lengths = _np.array([len(c) + 1 for c in circuits]) + if isinstance(bins, str) and 'auto' in bins and 'int' in bins: # type: ignore + bins = int(_np.log2(_np.max(lengths))) + if isinstance(trans, _np.ufunc): + lengths = trans(lengths) + elif trans == 'log': + lengths = _np.log2(lengths) + elif (trans != 'none') and (trans is not None): + raise ValueError(f'Argument `trans` had unsupported value, {trans}.') + # get bin edges from numpy, then drop empty bins. + counts, edges = _np.histogram(lengths, bins) + edges = _np.concatenate([[edges[0]], edges[1:][counts > 0]]) + assignments = _np.digitize(lengths, edges) + assignments -= 1 + # edges[ assignments[j] ] <= lengths[j] < edges[ assignments[j]+1 ] + num_bins = edges.size - 1 + circuit_lists = [list() for _ in range(num_bins)] + for j, c in zip(assignments, circuits): + for i in range(j, num_bins): + circuit_lists[i].append(c) + """ + # The following approch to building circuit_lists is (in theory) less efficient than + # the approach above, but it has the advantage of avoiding the nested for-loop. + circuit_lists = [] + last_size = 0 + edges = _np.histogram_bin_edges(lengths, bins) + for upperbound in edges[1:]: + cur_inds = _np.where(lengths < upperbound)[0] + if cur_inds.size == last_size: + continue # empty bin + last_size = cur_inds.size + circuit_lists.append([circuits[j] for j in cur_inds]) + """ + return circuit_lists # type: ignore + + +def logspaced_prefix_circuits( + c: _Circuit, + povms_to_keep: _Container[_LabelStr] = (_LabelStr('Mdefault'),), + base: Union[int, float]=2, + editable = False + ) -> list[_Circuit]: + + if len(c) > 0 and c[-1] in povms_to_keep: + M = c[-1] + if not isinstance(M, _LabelStr): + M = _LabelStr(M) + c = c[:-1] + circuits = logspaced_prefix_circuits(c, tuple(), base, editable=True) + for c in circuits: + c._labels.append(M) # type: ignore + c.done_editing() + return circuits + + if editable and not c._static: + c = c.copy(editable=True) # type: ignore + + circuits = [c] + assert base > 1 + next_len = int(len(c) // base) + while next_len > 0: + c = c[:next_len] + circuits.append(c) + next_len = int(len(c) // base) + + if not editable: + for c in circuits: + c.done_editing() + + return circuits + diff --git a/pygsti/report/factory.py b/pygsti/report/factory.py index c8781f53f..02f9266a6 100644 --- a/pygsti/report/factory.py +++ b/pygsti/report/factory.py @@ -1199,6 +1199,7 @@ def construct_standard_report(results, title="auto", in this particular report. Strings will be cast to lowercase, stripped of white space, and then mapped to omitted Section classes as follows + { 'summary' : SummarySection, 'goodness' : GoodnessSection, diff --git a/pygsti/report/section/summary.py b/pygsti/report/section/summary.py index 4d6764763..d4e522bad 100644 --- a/pygsti/report/section/summary.py +++ b/pygsti/report/section/summary.py @@ -28,8 +28,9 @@ def final_model_fit_progress_bar_plot_sum(workspace, switchboard=None, max_lengt def final_model_fit_histogram(workspace, switchboard=None, linlog_percentile=5, comm=None, bgcolor='white', **kwargs): return workspace.ColorBoxPlot( - switchboard.objfn_builder, switchboard.circuits_final, - switchboard.modvi_ds, switchboard.mdl_final_modvi, + switchboard.objfn_builder_modvi, + switchboard.circuits_final, + switchboard.modvi_ds, switchboard.mdl_current_modvi, linlg_pcntle=linlog_percentile / 100, typ='histogram', comm=comm, bgcolor=bgcolor, mdc_store= switchboard.final_mdc_store diff --git a/pygsti/report/workspaceplots.py b/pygsti/report/workspaceplots.py index dadd179b1..92a43d22c 100644 --- a/pygsti/report/workspaceplots.py +++ b/pygsti/report/workspaceplots.py @@ -4234,11 +4234,14 @@ def _create(self, x_names, circuits_by_x, model_by_x, dataset_by_x, objfn_builde if mdc_stores is None: mdc_stores = [None]*len(model_by_x) - for X, mdl, circuits, dataset, Np, mdc_store in zip(x_names, model_by_x, circuits_by_x, - dataset_by_x, np_by_x, mdc_stores): + for X, mdl, circuits, dataset, Np, mdc_store in zip( + x_names, model_by_x, circuits_by_x, dataset_by_x, np_by_x, mdc_stores + ): if circuits is None or mdl is None: Nsig, rating = _np.nan, 5 else: + # from pygsti.baseobjs.smartcache import _get_fn_name_key + # self.ws.smartCache.ineffective.add(_get_fn_name_key(_ph.rated_n_sigma)) Nsig, rating, _, _, _, _ = self._ccompute(_ph.rated_n_sigma, dataset, mdl, circuits, objfn_builder, Np, wildcard, return_all=True, comm=comm, mdc_store=mdc_store) diff --git a/pygsti/tools/optools.py b/pygsti/tools/optools.py index 5dab09e8a..a3dd4c00c 100644 --- a/pygsti/tools/optools.py +++ b/pygsti/tools/optools.py @@ -40,7 +40,18 @@ from typing import Union, Optional, Any BasisLike = Union[str, _Basis] -IMAG_TOL = 1e-7 # tolerance for imaginary part being considered zero + +__SCALAR_TOL_EXPONENT__ = 0.75 +"""^ +__SCALAR_TOL_EXPONENT__ is used when checking properties of matrices and vectors. +It's intended only when we can check the property without incurring rounding +errors from some reduction (like a sum or a matrix-vector product). If we're +working with a matrix whose dtype is `d`, then we set + + __SCALAR_TOL__ = d.eps ** __SCALAR_TOL_EXPONENT__, + +or a modest multiple thereof. +""" def _flat_mut_blks(i, j, block_dims): @@ -82,7 +93,7 @@ def fidelity(a, b): float The resulting fidelity. """ - __SCALAR_TOL__ = _np.finfo(a.dtype).eps ** 0.75 + __SCALAR_TOL__ = _np.finfo(a.dtype).eps ** __SCALAR_TOL_EXPONENT__ # ^ use for checks that have no dimensional dependence; about 1e-12 for double precision. __VECTOR_TOL__ = (a.shape[0] ** 0.5) * __SCALAR_TOL__ # ^ use for checks that do have dimensional dependence (will naturally increase for larger matrices) @@ -99,6 +110,18 @@ def assert_hermitian(mat): assert_hermitian(a) assert_hermitian(b) + def check_unit_trace(mat): + tr = _np.trace(mat) + if abs(tr - 1) > __VECTOR_TOL__: + message = f""" + The input matrix is trace {tr}, which deviates from 1 by more than {__VECTOR_TOL__}. + Beware result! + """ + _warnings.warn(message) + + check_unit_trace(a) + check_unit_trace(b) + def check_rank_one_density(mat): """ mat is Hermitian of order n. This function uses an O(n^2) time randomized algorithm to @@ -120,7 +143,11 @@ def check_rank_one_density(mat): """ n = mat.shape[0] - if _np.linalg.norm(mat) < __VECTOR_TOL__: + ZERO_THRESHOLD = n**0.5 * _np.finfo(mat.dtype).eps**0.75 + # ^ This value affects correctness, not just when we raise an error or a warning. + # Don't change from the value given here. + + if _np.linalg.norm(mat) < ZERO_THRESHOLD: # We prefer to return the zero vector instead of None to simplify how we handle # this function's output. return 0, _np.zeros(n, dtype=complex) @@ -134,7 +161,7 @@ def check_rank_one_density(mat): alpha = _np.real(candidate_v.conj() @ mat @ candidate_v) reconstruction = alpha * _np.outer(candidate_v, candidate_v.conj()) - if _np.linalg.norm(mat - reconstruction) > __VECTOR_TOL__: + if _np.linalg.norm(mat - reconstruction) > ZERO_THRESHOLD: # We can't certify that mat is rank-1. return 2, None @@ -144,9 +171,9 @@ def check_rank_one_density(mat): return 0, _np.zeros(n) if abs(alpha - 1) > __SCALAR_TOL__: - message = f"The input matrix is not trace-1 up to tolerance {__SCALAR_TOL__}. Beware result!" + message = f"The input matrix is not trace-1 up to tolerance {__SCALAR_TOL__}." _warnings.warn(message) - candidate_v *= _np.sqrt(alpha) + candidate_v *= _np.sqrt(alpha) return 1, candidate_v @@ -171,6 +198,7 @@ def psd_square_root(mat): if _np.min(evals) < -__SCALAR_TOL__: message = f""" Input matrix is not PSD up to tolerance {__SCALAR_TOL__}. + The negative eigenvalues are {evals[evals < 0]}. We'll project out the bad eigenspaces to only work with the PSD part. """ _warnings.warn(message) diff --git a/pygsti/tools/rbtheory.py b/pygsti/tools/rbtheory.py index bb0cd03ed..57558e936 100644 --- a/pygsti/tools/rbtheory.py +++ b/pygsti/tools/rbtheory.py @@ -326,6 +326,14 @@ def L_matrix(model, target_model, weights=None): # noqa N802 weights[key] = 1. normalizer = _np.sum(_np.array([weights[key] for key in list(target_model.operations.keys())])) + # TODO: improve efficiency + # + # 1. Accumuate the summands in this list comprehension in-place. (Might already happen but that's non-obvious) + # 2. Use the fact that target gates are unitary and so their superoperator representation inverses should + # be their transposes. + # 3. Have the option to return this matrix as an implicit abstract linear operator, so that anyone who wants + # eigenvalue info can try to get it from an iterative method instead of a full eigendecomposition. + # L_matrix = (1 / normalizer) * _np.sum( weights[key] * _np.kron( model.operations[key].to_dense("HilbertSchmidt").T, diff --git a/test/performance/mpi_2D_scaling/run_me_with_mpirun.py b/test/performance/mpi_2D_scaling/run_me_with_mpirun.py index b4b966c52..9f1514550 100644 --- a/test/performance/mpi_2D_scaling/run_me_with_mpirun.py +++ b/test/performance/mpi_2D_scaling/run_me_with_mpirun.py @@ -31,10 +31,14 @@ ds = ds_ref MINCLIP = 1e-4 -chi2_builder = pygsti.objectivefns.Chi2Function.builder( - 'chi2', regularization={'min_prob_clip_for_weighting': MINCLIP}, penalties={'cptp_penalty_factor': 0.0}) -mle_builder = pygsti.objectivefns.PoissonPicDeltaLogLFunction.builder( - 'logl', regularization={'min_prob_clip': MINCLIP, 'radius': MINCLIP}) +chi2_builder = pygsti.objectivefns.ObjectiveFunctionBuilder( + pygsti.objectivefns.Chi2Function, + 'chi2', regularization={'min_prob_clip_for_weighting': MINCLIP}, penalties={'cptp_penalty_factor': 0.0} +) +mle_builder = pygsti.objectivefns.ObjectiveFunctionBuilder( + pygsti.objectivefns.PoissonPicDeltaLogLFunction, + 'logl', regularization={'min_prob_clip': MINCLIP, 'radius': MINCLIP} +) iteration_builders = [chi2_builder]; final_builders = [mle_builder] builders = pygsti.protocols.GSTObjFnBuilders(iteration_builders, final_builders) diff --git a/test/test_packages/drivers/test_timedep.py b/test/test_packages/drivers/test_timedep.py index 7e7830cb9..8c726a6b1 100644 --- a/test/test_packages/drivers/test_timedep.py +++ b/test/test_packages/drivers/test_timedep.py @@ -1,4 +1,6 @@ import logging + +import pygsti.objectivefns mpl_logger = logging.getLogger('matplotlib') mpl_logger.setLevel(logging.WARNING) @@ -109,7 +111,8 @@ def test_time_dependent_gst_staticdata(self): target_model.sim = pygsti.forwardsims.MapForwardSimulator(max_cache_size=0) # No caching allowed for time-dependent calcs self.assertEqual(ds.degrees_of_freedom(aggregate_times=False), 57) - builders = pygsti.protocols.GSTObjFnBuilders([pygsti.objectivefns.TimeDependentPoissonPicLogLFunction.builder()], []) + builder = pygsti.objectivefns.ObjectiveFunctionBuilder(pygsti.objectivefns.TimeDependentPoissonPicLogLFunction) + builders = pygsti.protocols.GSTObjFnBuilders([builder], []) gst = pygsti.protocols.GateSetTomography(target_model, gaugeopt_suite=None, objfn_builders=builders, optimizer={'maxiter':2,'tol': 1e-4}) @@ -162,7 +165,8 @@ def test_time_dependent_gst(self): target_model.operations['Gi',0] = MyTimeDependentIdle(0) # start assuming no time dependent decay target_model.sim = pygsti.forwardsims.MapForwardSimulator(max_cache_size=0) # No caching allowed for time-dependent calcs - builders = pygsti.protocols.GSTObjFnBuilders([pygsti.objectivefns.TimeDependentPoissonPicLogLFunction.builder()], []) + builder = pygsti.objectivefns.ObjectiveFunctionBuilder(pygsti.objectivefns.TimeDependentPoissonPicLogLFunction) + builders = pygsti.protocols.GSTObjFnBuilders([builder], []) gst = pygsti.protocols.GateSetTomography(target_model, gaugeopt_suite=None, objfn_builders=builders, optimizer={'maxiter':10,'tol': 1e-4}) data = pygsti.protocols.ProtocolData(edesign, ds) diff --git a/test/test_packages/objects/test_hessian.py b/test/test_packages/objects/test_hessian.py index a22331a7f..8f3f24fa6 100644 --- a/test/test_packages/objects/test_hessian.py +++ b/test/test_packages/objects/test_hessian.py @@ -130,7 +130,7 @@ def test_confidenceRegion(self): res = proto.ModelEstimateResults(data, proto.StandardGST(modes="full TP")) #Add estimate for hessian-based CI -------------------------------------------------- - builder = pygsti.objectivefns.PoissonPicDeltaLogLFunction.builder() + builder = pygsti.objectivefns.ObjectiveFunctionBuilder(pygsti.objectivefns.PoissonPicDeltaLogLFunction) res.add_estimate( proto.estimate.Estimate.create_gst_estimate( res, smq1Q_XY.target_model(), smq1Q_XY.target_model(), diff --git a/test/unit/algorithms/test_core.py b/test/unit/algorithms/test_core.py index 737a5f7f8..af7c06765 100644 --- a/test/unit/algorithms/test_core.py +++ b/test/unit/algorithms/test_core.py @@ -6,7 +6,7 @@ from pygsti.baseobjs import Label from pygsti.circuits import Circuit, CircuitList from pygsti.objectivefns import Chi2Function, FreqWeightedChi2Function, \ - PoissonPicDeltaLogLFunction + PoissonPicDeltaLogLFunction, ObjectiveFunctionBuilder from . import fixtures from ..util import BaseCase @@ -117,7 +117,7 @@ def test_do_mc2gst(self): # TODO assert correctness def test_do_mc2gst_regularize_factor(self): - obj_builder = Chi2Function.builder( + obj_builder = ObjectiveFunctionBuilder(Chi2Function, name='chi2', description="Sum of chi^2", regularization={'min_prob_clip_for_weighting': 1e-4}, @@ -130,7 +130,7 @@ def test_do_mc2gst_regularize_factor(self): # TODO assert correctness def test_do_mc2gst_CPTP_penalty_factor(self): - obj_builder = Chi2Function.builder( + obj_builder = ObjectiveFunctionBuilder(Chi2Function, name='chi2', description="Sum of chi^2", regularization={'min_prob_clip_for_weighting': 1e-4}, @@ -143,7 +143,7 @@ def test_do_mc2gst_CPTP_penalty_factor(self): # TODO assert correctness def test_do_mc2gst_SPAM_penalty_factor(self): - obj_builder = Chi2Function.builder( + obj_builder = ObjectiveFunctionBuilder(Chi2Function, name='chi2', description="Sum of chi^2", regularization={'min_prob_clip_for_weighting': 1e-4}, @@ -156,7 +156,7 @@ def test_do_mc2gst_SPAM_penalty_factor(self): # TODO assert correctness def test_do_mc2gst_CPTP_SPAM_penalty_factor(self): - obj_builder = Chi2Function.builder( + obj_builder = ObjectiveFunctionBuilder(Chi2Function, name='chi2', description="Sum of chi^2", regularization={'min_prob_clip_for_weighting': 1e-4}, @@ -202,7 +202,7 @@ def test_do_iterative_mc2gst(self): # TODO assert correctness def test_do_iterative_mc2gst_regularize_factor(self): - obj_builder = Chi2Function.builder( + obj_builder = ObjectiveFunctionBuilder(Chi2Function, name='chi2', description="Sum of chi^2", regularization={'min_prob_clip_for_weighting': 1e-4}, @@ -218,7 +218,7 @@ def test_do_iterative_mc2gst_regularize_factor(self): # TODO assert correctness def test_do_iterative_mc2gst_use_freq_weighted_chi2(self): - obj_builder = FreqWeightedChi2Function.builder( + obj_builder = ObjectiveFunctionBuilder(FreqWeightedChi2Function, name='freq-weighted-chi2', description="Sum of chi^2", regularization={'min_freq_clip_for_weighting': 1e-4} @@ -269,7 +269,7 @@ def test_do_mlgst(self): # TODO assert correctness def test_do_mlgst_CPTP_penalty_factor(self): - obj_builder = PoissonPicDeltaLogLFunction.builder( + obj_builder = ObjectiveFunctionBuilder(PoissonPicDeltaLogLFunction, name='logl', description="2*DeltaLogL", regularization={'min_prob_clip': 1e-4}, @@ -283,7 +283,7 @@ def test_do_mlgst_CPTP_penalty_factor(self): # TODO assert correctness def test_do_mlgst_SPAM_penalty_factor(self): - obj_builder = PoissonPicDeltaLogLFunction.builder( + obj_builder = ObjectiveFunctionBuilder(PoissonPicDeltaLogLFunction, name='logl', description="2*DeltaLogL", regularization={'min_prob_clip': 1e-4}, @@ -302,7 +302,7 @@ def test_do_mlgst_CPTP_SPAM_penalty_factor(self): # FUTURE: see what we can do in custom LM about scaling large # jacobians... #self.skipTest("Ignore for now.") - obj_builder = PoissonPicDeltaLogLFunction.builder( + obj_builder = ObjectiveFunctionBuilder(PoissonPicDeltaLogLFunction, name='logl', description="2*DeltaLogL", regularization={'min_prob_clip': 1e-4}, @@ -356,7 +356,7 @@ def test_do_iterative_mlgst(self): # ) def test_do_iterative_mlgst_use_freq_weighted_chi2(self): - obj_builder = FreqWeightedChi2Function.builder( + obj_builder = ObjectiveFunctionBuilder(FreqWeightedChi2Function, name='freq-weighted-chi2', description="Sum of chi^2", regularization={'min_freq_clip_for_weighting': 1e-4} @@ -404,7 +404,7 @@ def test_do_mlgst_raises_on_out_of_memory(self): # XXX if this function needs explicit coverage, it should be public! def test_do_mlgst_base_forcefn_grad(self): forcefn_grad = np.ones((1, self.mdl_clgst.num_params), 'd') - obj_builder = PoissonPicDeltaLogLFunction.builder( + obj_builder = ObjectiveFunctionBuilder(PoissonPicDeltaLogLFunction, name='logl', description="2*DeltaLogL", regularization={'min_prob_clip': 1e-4}, diff --git a/test/unit/algorithms/test_gaugeopt_correctness.py b/test/unit/algorithms/test_gaugeopt_correctness.py index af4b39c70..3012bc7a7 100644 --- a/test/unit/algorithms/test_gaugeopt_correctness.py +++ b/test/unit/algorithms/test_gaugeopt_correctness.py @@ -41,15 +41,15 @@ def check_gate_metrics_are_nontrivial(metrics, tol): for lbl, val in metrics['infids'].items(): if lbl == idle_label: continue - assert val > tol, f"{val} is at most {tol}, failure for gate {lbl} w.r.t. infidelity." + assert val > tol, f"{val} is <= {tol}, failure for gate {lbl} w.r.t. infidelity." for lbl, val in metrics['frodists'].items(): if lbl == idle_label: continue - assert val > tol, f"{val} is at most {tol}, failure for gate {lbl} w.r.t. frobenius distance." + assert val > tol, f"{val} is <= {tol}, failure for gate {lbl} w.r.t. frobenius distance." for lbl, val in metrics['tracedists'].items(): if lbl == idle_label: continue - assert val > tol, f"{val} is at most {tol}, failure for gate {lbl} w.r.r. trace distance." + assert val > tol, f"{val} is <= {tol}, failure for gate {lbl} w.r.r. trace distance." return @@ -111,18 +111,20 @@ def _gauge_transform_model(self, seed): self.model = self.target.copy() self.model.default_gauge_group = self.default_gauge_group np.random.seed(seed) - self.U = la.expm(np.random.randn()/2 * -1j * (pgbc.sigmax + pgbc.sigmaz)/np.sqrt(2)) + strength = (np.random.rand() + 1)/100 + self.U = la.expm(strength * -1j * (pgbc.sigmax + pgbc.sigmaz)/np.sqrt(2)) self.gauge_grp_el = self.gauge_grp_el_class(pgo.unitary_to_pauligate(self.U)) self.model.transform_inplace(self.gauge_grp_el) self.metrics_before = gate_metrics_dict(self.model, self.target) - check_gate_metrics_are_nontrivial(self.metrics_before, tol=1e-2) + check_gate_metrics_are_nontrivial(self.metrics_before, tol=1e-4) return def _main_tester(self, method, seed, test_tol, alg_tol, gop_objective): assert gop_objective in {'frobenius', 'frobenius_squared', 'tracedist', 'fidelity'} self._gauge_transform_model(seed) tic = time.time() - newmodel = gop.gaugeopt_to_target(self.model, self.target, method=method, tol=alg_tol, spam_metric=gop_objective, gates_metric=gop_objective) + verbosity = 0 + newmodel = gop.gaugeopt_to_target(self.model, self.target, method=method, tol=alg_tol, spam_metric=gop_objective, gates_metric=gop_objective, verbosity=verbosity) toc = time.time() dt = toc - tic metrics_after = gate_metrics_dict(newmodel, self.target) @@ -133,7 +135,7 @@ def test_lbfgs_frodist(self): self.setUp() times = [] for seed in range(self.N_REPS): - dt = self._main_tester('L-BFGS-B', seed, test_tol=1e-6, alg_tol=1e-8, gop_objective='frobenius') + dt = self._main_tester('L-BFGS-B', seed, test_tol=1e-5, alg_tol=1e-7, gop_objective='frobenius') times.append(dt) print(f'GaugeOpt over {self.gauge_space_name} w.r.t. Frobenius dist, using L-BFGS: {times}.') return @@ -142,7 +144,7 @@ def test_lbfgs_tracedist(self): self.setUp() times = [] for seed in range(self.N_REPS): - dt = self._main_tester('L-BFGS-B', seed, test_tol=1e-6, alg_tol=1e-8, gop_objective='tracedist') + dt = self._main_tester('L-BFGS-B', seed, test_tol=1e-5, alg_tol=1e-7, gop_objective='tracedist') times.append(dt) print(f'GaugeOpt over {self.gauge_space_name} w.r.t. trace dist, using L-BFGS: {times}.') return diff --git a/test/unit/objects/fixtures.py b/test/unit/objects/fixtures.py index 1bf1eab6a..44b3b80f3 100644 --- a/test/unit/objects/fixtures.py +++ b/test/unit/objects/fixtures.py @@ -65,9 +65,12 @@ def lsgstStructs(self): @ns.memo def mdl_lsgst(self): - chi2_builder = pygsti.objectivefns.Chi2Function.builder( + from pygsti.objectivefns.objectivefns import ObjectiveFunctionBuilder, Chi2Function + chi2_builder = ObjectiveFunctionBuilder( + Chi2Function, regularization={'min_prob_clip_for_weighting': 1e-6}, - penalties={'prob_clip_interval': (-1e6, 1e6)}) + penalties={'prob_clip_interval': (-1e6, 1e6)} + ) models, _, _, _ = pygsti.algorithms.core.run_iterative_gst( self.dataset, self.mdl_clgst, self.lsgstStrings, optimizer=None, diff --git a/test/unit/objects/test_objectivefns.py b/test/unit/objects/test_objectivefns.py index df5873705..d1e3650a1 100644 --- a/test/unit/objects/test_objectivefns.py +++ b/test/unit/objects/test_objectivefns.py @@ -270,12 +270,6 @@ def setUp(self): super().setUp() self.objfns = self.build_objfns() - def test_builder(self): - #All objective function should be of the same type - cls = self.objfns[0].__class__ - builder = cls.builder("test_name", "test_description") - self.assertTrue(isinstance(builder, _objfns.ObjectiveFunctionBuilder)) - def test_value(self): for objfn in self.objfns: terms = objfn.terms().copy() diff --git a/test/unit/protocols/test_gst.py b/test/unit/protocols/test_gst.py index 2b712f868..9365c61dd 100644 --- a/test/unit/protocols/test_gst.py +++ b/test/unit/protocols/test_gst.py @@ -2,7 +2,7 @@ from pygsti.forwardsims.mapforwardsim import MapForwardSimulator from pygsti.modelpacks import smq1Q_XYI from pygsti.modelpacks.legacy import std1Q_XYI, std2Q_XYICNOT -from pygsti.objectivefns.objectivefns import PoissonPicDeltaLogLFunction +from pygsti.objectivefns.objectivefns import PoissonPicDeltaLogLFunction, ObjectiveFunctionBuilder from pygsti.models.gaugegroup import TrivialGaugeGroup from pygsti.objectivefns import FreqWeightedChi2Function from pygsti.optimize.simplerlm import SimplerLMOptimizer @@ -64,7 +64,7 @@ def test_gaugeopt_suite_raises_on_bad_suite(self): GSTGaugeOptSuite("foobar").to_dictionary(model_1Q, verbosity=1) def test_add_badfit_estimates(self): - builder = PoissonPicDeltaLogLFunction.builder() + builder = ObjectiveFunctionBuilder(PoissonPicDeltaLogLFunction) opt = SimplerLMOptimizer() badfit_opts = gst.GSTBadFitOptions(threshold=-10, actions=("robust", "Robust", "robust+", "Robust+", "wildcard", "do nothing")) diff --git a/test/unit/tools/fixtures.py b/test/unit/tools/fixtures.py index d24c77c70..6792fc723 100644 --- a/test/unit/tools/fixtures.py +++ b/test/unit/tools/fixtures.py @@ -58,9 +58,12 @@ def lsgstStrings(self): @ns.memo def mdl_lsgst(self): - chi2_builder = pygsti.objectivefns.Chi2Function.builder( + from pygsti.objectivefns.objectivefns import ObjectiveFunctionBuilder, Chi2Function + chi2_builder = ObjectiveFunctionBuilder( + Chi2Function, regularization={'min_prob_clip_for_weighting': 1e-6}, - penalties={'prob_clip_interval': (-1e6, 1e6)}) + penalties={'prob_clip_interval': (-1e6, 1e6)} + ) models, _, _, _ = pygsti.algorithms.core.run_iterative_gst( self.dataset, self.mdl_clgst, self.lsgstStrings, optimizer=None, diff --git a/test/unit/util.py b/test/unit/util.py index 54c930632..6bf0ef5d8 100644 --- a/test/unit/util.py +++ b/test/unit/util.py @@ -6,18 +6,25 @@ from contextlib import contextmanager from pathlib import Path from tempfile import TemporaryDirectory +import importlib.util import numpy as np +__cvxpy_not_importable__ = importlib.util.find_spec('cvxpy') is None + +__deap_not_importable__ = importlib.util.find_spec('deap') is None + def needs_cvxpy(fn): """Shortcut decorator for skipping tests that require CVXPY""" - return unittest.skipIf('SKIP_CVXPY' in os.environ, "skipping cvxpy tests")(fn) + cond = __cvxpy_not_importable__ or ('SKIP_CVXPY' in os.environ) + return unittest.skipIf(cond, "skipping cvxpy tests")(fn) def needs_deap(fn): """Shortcut decorator for skipping tests that require deap""" - return unittest.skipIf('SKIP_DEAP' in os.environ, "skipping deap tests")(fn) + cond = __deap_not_importable__ or ('SKIP_DEAP' in os.environ) + return unittest.skipIf(cond, "skipping deap tests")(fn) def needs_matplotlib(fn): diff --git a/wip_notebook_sharing/README.md b/wip_notebook_sharing/README.md new file mode 100644 index 000000000..6a233cd35 --- /dev/null +++ b/wip_notebook_sharing/README.md @@ -0,0 +1,4 @@ +# This folder + +Contents of this folder should be deleted before merge. + diff --git a/wip_notebook_sharing/case0-gst-with-outliers.ipynb b/wip_notebook_sharing/case0-gst-with-outliers.ipynb new file mode 100644 index 000000000..f3e993ad8 --- /dev/null +++ b/wip_notebook_sharing/case0-gst-with-outliers.ipynb @@ -0,0 +1,318 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "%load_ext autoreload\n", + "%autoreload 2\n", + "\n", + "import pygsti\n", + "from pygsti.modelpacks import smq1Q_XYI\n", + "import numpy as np\n", + "from pprint import pprint\n", + "import copy\n", + "\n", + "\n", + "def make_depolarized_dataset(modelpack, depol_level=0.01, max_max_len=128):\n", + " ideal_model = modelpack.target_model()\n", + " prep_fids = modelpack.prep_fiducials()\n", + " meas_fids = modelpack.meas_fiducials()\n", + " germs = modelpack.germs()\n", + " max_lens = [2**p for p in range(1+int(np.log2(max_max_len)))]\n", + " lsgst_circuit_lists = pygsti.circuits.create_lsgst_circuit_lists(ideal_model, prep_fids, meas_fids, germs, max_lens)\n", + " all_circuits = lsgst_circuit_lists[-1]\n", + " shots_per_circuit = 1000\n", + " rng_state = np.random.default_rng(0)\n", + " depol_model = ideal_model.depolarize(op_noise=depol_level)\n", + " ds = pygsti.data.simulate_data(depol_model, all_circuits, shots_per_circuit, rand_state=rng_state)\n", + " return ds\n", + "\n", + "\n", + "def corrupt_dataset(ds, prop_corrupt, rng=0):\n", + " dsc = ds.copy_nonstatic()\n", + " rng = np.random.default_rng(rng)\n", + " num_circs = len(dsc)\n", + " selected = rng.choice(np.arange(num_circs), size=int(num_circs*prop_corrupt), replace=False)\n", + " circuits = list(dsc.keys())\n", + " selected = [circuits[i] for i in selected]\n", + " for c in selected:\n", + " num_shots = dsc[c].total\n", + " old_row = dsc[c].to_dict()\n", + " distn = rng.random(len(old_row))\n", + " distn /= np.sum(distn)\n", + " new_row = {k: num_shots * distn[i] for i,k in enumerate(old_row.keys())}\n", + " dsc[c] = new_row\n", + " dsc.comment = 'corrupt'\n", + " return dsc, selected\n", + "\n", + "\n", + "def cptp_gst(ds, fids, germs, target_model, final_objective: str, verbosity: int, mode='CPTPLND'):\n", + " \"\"\"\n", + " In the context of this notebook, `ds` is produced by either make_depolarized_dataset or corrupt_dataset.\n", + " final_objective can be 'tvd', 'chi2', or 'logl'.\n", + "\n", + " This function wraps up three steps of a GST pipeline.\n", + "\n", + " 1. Construct a StandardGSTDesign based on (target_model, ds, fids, germs).\n", + " * processor_spec is the value returned from target_model.create_processor_spec.\n", + " * max_lens list is all powers of two that are <= the depth of the longest circuit in ds.\n", + " * circuits in the design are filtered to only include circuits that appeared in ds.\n", + "\n", + " 2. Construct a StandardGST protocol object based on (final_objective, mode, verbosity).\n", + " * The gauge optimization suite is 'stdgaugeopt', minus the TPSpam optimization step.\n", + " * objfn_builders, optimizer, and badfit_options are all set so the final \n", + " iteration's objective function is based on final_objective.\n", + "\n", + " 3. Run GST with checkpointing turned off. \n", + " We dot NOT save the results to disk! The calling function is responsible for that.\n", + " \"\"\"\n", + " max_exp = int(np.log2(np.max([len(c) for c in ds.keys()])))\n", + " max_lens = [2**p for p in range(1 + max_exp)]\n", + " prep_fids, meas_fids = fids\n", + "\n", + " target_model = target_model.copy()\n", + " target_model.default_gauge_group = 'unitary'\n", + "\n", + " gos = pygsti.protocols.gst.GSTGaugeOptSuite.cast('stdgaugeopt')\n", + " gop_params = gos.to_dictionary(target_model)\n", + " # ^ a dict with one key, 'stdgaugeopt', whose corresponding value is a list of dicts.\n", + " # The internal dicts will indicate Frobenius-based losses for gates and SPAM,\n", + " # along with varying weights. Additional elements can be added to any one of these\n", + " # internal dicts to be passed to gaugeopt_to_target.\n", + " gop_params['stdgaugeopt'] = gop_params['stdgaugeopt'][:-1]\n", + " # ^ drop the 1-dimensional TPSpam gauge optimization step.\n", + "\n", + " exp_design = pygsti.protocols.StandardGSTDesign(\n", + " target_model.create_processor_spec(),\n", + " prep_fids, meas_fids, germs, max_lens,\n", + " None, # germ_length_limits\n", + " None, 1, None, # fidPairs, keepFraction, keepSeed\n", + " True, True, # include_lgst, nested_circuit_lists\n", + " None, # string_manipulation_rules\n", + " None, # op_label_aliases\n", + " ds, 'drop', verbosity=verbosity\n", + " )\n", + " data = pygsti.protocols.ProtocolData(exp_design, ds)\n", + " \n", + " from pygsti.drivers.longsequence import _get_gst_builders, _get_optimizer, _get_badfit_options\n", + " advanced_options = {'objective': final_objective}\n", + " proto = pygsti.protocols.StandardGST(\n", + " (mode,), gop_params, target_model, None, \n", + " objfn_builders = _get_gst_builders(advanced_options),\n", + " optimizer = _get_optimizer(advanced_options, target_model),\n", + " badfit_options = _get_badfit_options(advanced_options),\n", + " verbosity = verbosity\n", + " )\n", + " results = proto.run(data, disable_checkpointing=True)\n", + " return results\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "mp = smq1Q_XYI\n", + "target = mp.target_model()\n", + "fids = (mp.prep_fiducials(), mp.meas_fiducials())\n", + "germs = mp.germs()\n", + "ds = make_depolarized_dataset(mp, depol_level=0.01, max_max_len=128)" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "--- Circuit Creation ---\n", + "-- Std Practice: [--------------------------------------------------] 0.0% (full TP) --\r" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/rjmurr/Documents/pygsti-tvd/pygsti/pygsti/objectivefns/objectivefns.py:4483: RuntimeWarning: divide by zero encountered in divide\n", + " p5over_lsvec = 0.5/lsvec\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "-- Std Practice: [##################################################] 100.0% (full TP) --\n", + "--- Circuit Creation ---\n", + "-- Std Practice: [--------------------------------------------------] 0.0% (full TP) --\r" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/rjmurr/Documents/pygsti-tvd/pygsti/pygsti/objectivefns/objectivefns.py:4014: UserWarning: This derivative is discontinuous and does not return a full subgradient.\n", + " _warnings.warn('This derivative is discontinuous and does not return a full subgradient.')\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "-- Std Practice: [##################################################] 100.0% (full TP) --\n", + "--- Circuit Creation ---\n", + "-- Std Practice: [--------------------------------------------------] 0.0% (full TP) --\r" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/rjmurr/Documents/pygsti-tvd/pygsti/pygsti/objectivefns/objectivefns.py:4483: RuntimeWarning: divide by zero encountered in divide\n", + " p5over_lsvec = 0.5/lsvec\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "-- Std Practice: [##################################################] 100.0% (full TP) --\n", + "--- Circuit Creation ---\n", + "-- Std Practice: [##################################################] 100.0% (full TP) --\n" + ] + } + ], + "source": [ + "fit_mode = 'full TP'\n", + "\n", + "results = cptp_gst(ds, fids, germs, target, 'logl', verbosity=1, mode=fit_mode)\n", + "results.estimates['fit-true-logl'] = results.estimates.pop(fit_mode)\n", + "tvd_res = cptp_gst(ds, fids, germs, target, 'tvd', verbosity=1, mode=fit_mode)\n", + "results.estimates['fit-true-tvd'] = tvd_res.estimates.pop(fit_mode)\n", + "\n", + "dsc, selected = corrupt_dataset(ds, prop_corrupt=0.025)\n", + "logl_res = cptp_gst(dsc, fids, germs, target, 'logl', verbosity=1, mode=fit_mode)\n", + "results.estimates['fit-corrupt-logl'] = logl_res.estimates.pop(fit_mode)\n", + "tvd_res = cptp_gst(dsc, fids, germs, target, 'tvd', verbosity=1, mode=fit_mode)\n", + "results.estimates['fit-corrupt-tvd'] = tvd_res.estimates.pop(fit_mode)" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Running idle tomography\n", + "Computing switchable properties\n", + "Found standard clifford compilation from smq1Q_XYI\n", + "Found standard clifford compilation from smq1Q_XYI\n", + "Found standard clifford compilation from smq1Q_XYI\n", + "Found standard clifford compilation from smq1Q_XYI\n", + "Found standard clifford compilation from smq1Q_XYI\n", + "Found standard clifford compilation from smq1Q_XYI\n", + "Found standard clifford compilation from smq1Q_XYI\n", + "Found standard clifford compilation from smq1Q_XYI\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/Users/rjmurr/Documents/pygsti-tvd/pygsti/pygsti/tools/optools.py:164: UserWarning:\n", + "\n", + "\n", + " Input matrix is not PSD up to tolerance 1.8189894035458565e-12.\n", + " We'll project out the bad eigenspaces to only work with the PSD part.\n", + " \n", + "\n", + "/Users/rjmurr/Documents/pygsti-tvd/pygsti/pygsti/tools/optools.py:172: UserWarning:\n", + "\n", + "\n", + " The PSD part of the input matrix is not trace-1 up to tolerance 3.637978807091713e-12.\n", + " Beware result!\n", + " \n", + "\n", + "/Users/rjmurr/Documents/pygsti-tvd/pygsti/pygsti/forwardsims/mapforwardsim.py:732: UserWarning:\n", + "\n", + "Generating dense process matrix representations of circuits or gates \n", + "can be inefficient and should be avoided for the purposes of forward \n", + "simulation/calculation of circuit outcome probability distributions \n", + "when using the MapForwardSimulator.\n", + "\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Statistical hypothesis tests did NOT find inconsistency between the data at 5.00% significance.\n", + "The data are INCONSISTENT at 5.00% significance.\n", + " - Details:\n", + " - The aggregate log-likelihood ratio test is significant at 167.04 standard deviations.\n", + " - The aggregate log-likelihood ratio test standard deviations signficance threshold is 2.00\n", + " - The number of sequences with data that is inconsistent is 21\n", + " - The maximum SSTVD over all sequences is 0.78\n", + " - The maximum SSTVD was observed for Qubit 0 ---|Gxpi2|-|Gxpi2|-|Gypi2|-|Gypi2|-|Gypi2|-|Gypi2|---\n", + "\n", + "The data are INCONSISTENT at 5.00% significance.\n", + " - Details:\n", + " - The aggregate log-likelihood ratio test is significant at 167.04 standard deviations.\n", + " - The aggregate log-likelihood ratio test standard deviations signficance threshold is 2.00\n", + " - The number of sequences with data that is inconsistent is 21\n", + " - The maximum SSTVD over all sequences is 0.78\n", + " - The maximum SSTVD was observed for Qubit 0 ---|Gxpi2|-|Gxpi2|-|Gypi2|-|Gypi2|-|Gypi2|-|Gypi2|---\n", + "\n", + "Statistical hypothesis tests did NOT find inconsistency between the data at 5.00% significance.\n" + ] + } + ], + "source": [ + "moar_results = results.copy()\n", + "dsc.comment = 'corrupt'\n", + "moar_results.data.dataset = dsc\n", + "report = pygsti.report.construct_standard_report(\n", + " {'eval-true' : results,\n", + " 'eval-corrupt' : moar_results\n", + " },\n", + " advanced_options={'skip_sections': ('colorbox',)},\n", + " title=\"GST Example Report\", verbosity=2\n", + ")\n", + "# NOTE: can reach in and change the entry in report.switchboard.objfn_builder_modvi,\n", + "# or anything else in the switchboard, according to my whims.\n", + "report.write_html(\"case0_reports_250217/exampleReport\", auto_open=True, verbosity=1)\n" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "rogst", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.11" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +}