diff --git a/pygsti/drivers/longsequence.py b/pygsti/drivers/longsequence.py index e6de76803..9785dff8b 100644 --- a/pygsti/drivers/longsequence.py +++ b/pygsti/drivers/longsequence.py @@ -385,7 +385,13 @@ 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. I've removed documentation for options + that we never reference in the code (truncScheme, estimate_label, and + missingDataAction). + + - objective = {'chi2', 'logl', 'tvd'} - op_labels = list of strings - circuit_weights = dict or None - starting_point = "LGST-if-possible" (default), "LGST", or "target" @@ -401,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 34e907d8f..a74ff89b5 100644 --- a/pygsti/layouts/copalayout.py +++ b/pygsti/layouts/copalayout.py @@ -207,12 +207,16 @@ def __init__(self, circuits, unique_circuits, to_unique, elindex_outcome_tuples, self._outcomes = dict() self._element_indices = dict() + self._element_inditer = dict() + ind_array = _np.arange(self._size) 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 + self._element_inditer[i_unique] = ind_array[s] # def hotswap_circuits(self, circuits, unique_complete_circuits=None): # self.circuits = circuits if isinstance(circuits, _CircuitList) else _CircuitList(circuits) @@ -744,7 +748,11 @@ 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]): + try: + iterator = zip(self._element_indices[i], self._outcomes[i]) + except TypeError: + iterator = zip(self._element_inditer[i], self._outcomes[i]) + for element_index, outcome in iterator: yield element_index, circuit, outcome def iter_unique_circuits(self): diff --git a/pygsti/models/explicitmodel.py b/pygsti/models/explicitmodel.py index b1f601298..75dc6e906 100644 --- a/pygsti/models/explicitmodel.py +++ b/pygsti/models/explicitmodel.py @@ -234,6 +234,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 diff --git a/pygsti/objectivefns/objectivefns.py b/pygsti/objectivefns/objectivefns.py index fd24655c1..a9d645584 100644 --- a/pygsti/objectivefns/objectivefns.py +++ b/pygsti/objectivefns/objectivefns.py @@ -555,6 +555,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 +591,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 +630,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: @@ -2671,36 +2675,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 +3943,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 +4011,11 @@ 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!") + _warnings.warn('This derivative is discontinuous and does not return a full subgradient.') + 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): """ @@ -5208,6 +5220,109 @@ def create_from(cls, model, dataset, circuits, regularization=None, penalties=No 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) + self.num_circuits = 0 + self.terms_weights = _np.array([]) + self._update_terms_weights() # <-- properly computes the two properties above. + + def _update_terms_weights(self): + # if self.num_circuits == self.layout.num_circuits: + # # Assume there's nothing to do. + # return + # else, assume that self.layout has been updated. + self.num_circuits = self.layout.num_circuits + 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. + return + + # QUESTION: could I just call the base class's implementation of terms and then scale the + # result? If so then I can avoid a ton of code duplication. + 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() + + unit_ralloc = self.layout.resource_alloc('atom-processing') + shared_mem_leader = unit_ralloc.is_host_leader + + 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 + + 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) + + 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() + + 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 + + # QUESTION: could I just call the base class's implementation of dterms and then scale the + # leading rows of the result? If so then I can avoid a ton of code duplication. + 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) + + 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] + + if shared_mem_leader and self._process_penalties: + self._dterms_fill_penalty(paramvec, self.jac[self.nelements:, :]) + + unit_ralloc.host_comm_barrier() + self.raw_objfn.resource_alloc.profiler.add_time("JACOBIAN", tm) + return self.jac class TimeDependentMDCObjectiveFunction(MDCObjectiveFunction): diff --git a/pygsti/protocols/gst.py b/pygsti/protocols/gst.py index dba21e8e6..b708589b1 100644 --- a/pygsti/protocols/gst.py +++ b/pygsti/protocols/gst.py @@ -785,7 +785,6 @@ def create_from(cls, objective='logl', freq_weighted_chi2=False, always_perform_ if 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] @@ -793,6 +792,9 @@ def create_from(cls, objective='logl', freq_weighted_chi2=False, always_perform_ else: iteration_builders = [chi2_builder] final_builders = [mle_builder] + elif objective == "tvd": + iteration_builders = [chi2_builder] + final_builders = [_objfns.ObjectiveFunctionBuilder.create_from('tvd')] else: raise ValueError("Invalid objective: %s" % objective) return cls(iteration_builders, final_builders) diff --git a/pygsti/report/factory.py b/pygsti/report/factory.py index 8d2f675d7..15f5670e4 100644 --- a/pygsti/report/factory.py +++ b/pygsti/report/factory.py @@ -184,8 +184,7 @@ def _create_master_switchboard(ws, results_dict, confidence_level, Ls = None for results in results_dict.values(): - est_labels = _add_new_estimate_labels(est_labels, results.estimates, - combine_robust) + est_labels = _add_new_estimate_labels(est_labels, results.estimates, combine_robust) loc_Ls = results.circuit_lists['final'].xs \ if isinstance(results.circuit_lists['final'], _PlaquetteGridCircuitStructure) else [0] Ls = _add_new_labels(Ls, loc_Ls) @@ -305,10 +304,9 @@ def _create_master_switchboard(ws, results_dict, confidence_level, else: est_modvi = est - switchBd.objfn_builder[d, i] = est.parameters.get( - 'final_objfn_builder', _objfns.ObjectiveFunctionBuilder.create_from('logl')) - switchBd.objfn_builder_modvi[d, i] = est_modvi.parameters.get( - 'final_objfn_builder', _objfns.ObjectiveFunctionBuilder.create_from('logl')) + switchBd.objfn_builder[d, i] = _objfns.ObjectiveFunctionBuilder.create_from('logl') # est.parameters.get('final_objfn_builder', _objfns.ObjectiveFunctionBuilder.create_from('logl')) + switchBd.objfn_builder_modvi[d, i] = _objfns.ObjectiveFunctionBuilder.create_from('logl') + # est_modvi.parameters.get('final_objfn_builder', _objfns.ObjectiveFunctionBuilder.create_from('logl')) switchBd.params[d, i] = est.parameters switchBd.clifford_compilation[d, i] = est.parameters.get("clifford compilation", 'auto') @@ -1172,6 +1170,32 @@ def construct_standard_report(results, title="auto", - idt_idle_oplabel : Label, optional The label identifying the idle gate (for use with idle tomography). + - skip_sections : tuple[str], optional + Contains names of standard report sections that should be skipped + 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, + 'colorbox' : GoodnessColorBoxPlotSection, + 'invariantgates' : GaugeInvariantsGatesSection, + 'invariantgerms' : GaugeInvariantsGermsSection, + 'variant' : GaugeVariantSection, + 'variantraw' : GaugeVariantsRawSection, + 'variantdecomp' : GaugeVariantsDecompSection, + 'varianterrorgen' : GaugeVariantsErrorGenSection, + 'input' : InputSection, + 'meta' : MetaSection, + 'help' : HelpSection + } + + A KeyError will be raised if skip_sections contains a string + that is not in the keys of the above dict (after casting to + lower case and stripping white space). + + verbosity : int, optional How much detail to send to stdout. @@ -1240,20 +1264,30 @@ def construct_standard_report(results, title="auto", flags.add('CombineRobust') # build section list - sections = [ - _section.SummarySection(), - _section.GoodnessSection(), - _section.GoodnessColorBoxPlotSection(), - _section.GaugeInvariantsGatesSection(), - _section.GaugeInvariantsGermsSection(), - _section.GaugeVariantSection(), - _section.GaugeVariantsRawSection(), - _section.GaugeVariantsDecompSection(), - _section.GaugeVariantsErrorGenSection(), - _section.InputSection(), - _section.MetaSection(), - _section.HelpSection() - ] + possible_sections = { + 'summary' : _section.SummarySection(), + 'goodness' : _section.GoodnessSection(), + 'colorbox' : _section.GoodnessColorBoxPlotSection(), + 'invariantgates' : _section.GaugeInvariantsGatesSection(), + 'invariantgerms' : _section.GaugeInvariantsGermsSection(), + 'variant' : _section.GaugeVariantSection(), + 'variantraw' : _section.GaugeVariantsRawSection(), + 'variantdecomp' : _section.GaugeVariantsDecompSection(), + 'varianterrorgen' : _section.GaugeVariantsErrorGenSection(), + 'input' : _section.InputSection(), + 'meta' : _section.MetaSection(), + 'help' : _section.HelpSection() + } + + skip_sections = advanced_options.get('skip_sections', tuple()) + if skip_sections: + if isinstance(skip_sections, str): + skip_sections = [skip_sections] + skip_sections = [s.lower().replace(' ','') for s in skip_sections] + for s in skip_sections: + possible_sections.pop(s) + sections = list(possible_sections.values()) + # ^ This whole process won't affect ordering of objects in "sections". if 'ShowScaling' in flags: sections.append(_section.GoodnessScalingSection()) diff --git a/pygsti/report/section/summary.py b/pygsti/report/section/summary.py index 9ac7263f4..90a47b74a 100644 --- a/pygsti/report/section/summary.py +++ b/pygsti/report/section/summary.py @@ -28,7 +28,8 @@ 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.objfn_builder_modvi, # NOTE: this should objfun_builder_modvi + switchboard.circuits_final, switchboard.modvi_ds, switchboard.mdl_current_modvi, linlg_pcntle=linlog_percentile / 100, typ='histogram', comm=comm, bgcolor=bgcolor 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 +}