Now that we've covered what circuits, models, and data sets are, let's see what we can do with them! This tutorial is intended to be an overview of the things that pyGSTi is able to do, with links to more detailed explanations and demonstrations as is appropriate. We begin with the simpler applications and proceed to the more complex ones. Here's a table of contents to give you a sense of what's here and so you can skip around if you'd like. Each of the sections here can more-or-less stand on its own.
We'll begin by setting up a Workspace
so we can display pretty interactive figures inline (see the intro to Workspaces tutorial for more details).
import pygsti
import numpy as np
ws = pygsti.report.Workspace()
ws.init_notebook_mode(autodisplay=True)
One of the simplest uses of pyGSTi is to construct a Model
and use it to compute the outcome probabilities of one or more Circuit
objects. This is generally accomplished using the .probs
method of a Model
object as shown below (this was also demonstrated in the essential objects tutorial). The real work is in constructing the Circuit
and Model
objects, which is covered in more detail in the circuits tutorial and in the explicit-model (best for 1-2 qubits) and implicit-model (best for 3+ qubits) tutorials.
PyGSTi refers to the process of computing circuit-outcome probabilities as forward simulation, and there are several methods of forward simulation available in pyGSTi. For example, one method multiplies together dense process matrices while another performs sparse matrix-vector products. A Model
is constructed for a single type of forward simulation, and it stores this within its .simtype
member. For more information on using different types of forward simulation see the forward simulation types tutorial.
mdl = pygsti.construction.build_explicit_model((0,1),
[(), ('Gx',0), ('Gy',0), ('Gx',1), ('Gy',1), ('Gcnot',0,1)],
["I(0,1)","X(pi/2,0)", "Y(pi/2,0)", "X(pi/2,1)", "Y(pi/2,1)", "CNOT(0,1)"])
c = pygsti.objects.Circuit([('Gx',0),('Gcnot',0,1),('Gy',1)] , line_labels=[0,1])
print("mdl will simulate probabilities using the '%s' forward-simulation method." % mdl.simtype)
mdl.probs(c) # Compute the outcome probabilities of circuit `c`
mdl will simulate probabilities using the 'matrix' forward-simulation method.
OutcomeLabelDict([(('00',), 0.25000000000000033), (('01',), 0.25000000000000006), (('10',), 0.25), (('11',), 0.25000000000000006)])
Only slightly more complex than computing circuit-outcome probabilities is generating simulated data (outcome counts rather than probabilities). This is performed by the generate_fake_data
function, which just samples the circuit-outcome probability distribution. You supply a list of Circuits
, a number of samples, and often times a seed to initialize the random sampler. This is an easy way to create a DataSet
to test other pyGSTi functions or to use independently.
The default behavior is to sample the multinomial distribution associated with the given outcome probabilities and number of samples. It's possible to turn off finite-sample error altogether and make the data-set counts exactly equal the probability values multiplied by the number of samples by setting sampleError='none'
, as demonstrated below.
circuit_list = pygsti.construction.circuit_list([ (),
(('Gx',0),),
(('Gx',0),('Gy',1)),
(('Gx',0),)*4,
(('Gx',0),('Gcnot',0,1)) ], line_labels=(0,1))
ds_fake = pygsti.construction.generate_fake_data(mdl, circuit_list, nSamples=100,
sampleError='multinomial', seed=1234)
print("Normal:")
print(ds_fake)
ds_nosampleerr = pygsti.construction.generate_fake_data(mdl, circuit_list, nSamples=100,
sampleError='none')
print("Without any sample error:")
print(ds_nosampleerr)
Normal: {} : {('00',): 100.0} Gx:0@(0,1) : {('00',): 47.0, ('10',): 53.0} Gx:0Gy:1@(0,1) : {('00',): 16.0, ('01',): 31.0, ('10',): 24.0, ('11',): 29.0} Gx:0Gx:0Gx:0Gx:0@(0,1) : {('00',): 100.0} Gx:0Gcnot:0:1@(0,1) : {('00',): 56.0, ('11',): 44.0} Without any sample error: {} : {('00',): 100.0} Gx:0@(0,1) : {('00',): 50.0, ('10',): 50.0} Gx:0Gy:1@(0,1) : {('00',): 25.0, ('01',): 25.0, ('10',): 25.0, ('11',): 25.0} Gx:0Gx:0Gx:0Gx:0@(0,1) : {('00',): 100.0, ('10',): -2.220446049250313e-14} Gx:0Gcnot:0:1@(0,1) : {('00',): 50.0, ('11',): 50.0}
The above section showed how the circuit-outcome probabilities computed by a Model
object can be used to generate data. We can also compare these probabilities with the outcome counts in an existing DataSet
, that is, ask the question: "For each circuit, how well do the frequencies of the outcomes (in the data) align with the probabilities predicted by the model?". There are several common statistics for this purpose; the two used most often in pyGSTi are the $\chi^2$ and log-likelihood ($\log\mathcal{L}$) statistics. If you're not sure what these are, the Methods section of this paper provides some details and here are a few practical considerations:
Here's how we compute the $\chi^2$ and $2\Delta\log\mathcal{L}$ between some data and a model:
import numpy as np
dataset_txt = \
"""## Columns = 00 count, 01 count, 10 count, 11 count
{} 100 0 0 0
Gx:0 55 5 40 0
Gx:0Gy:1 20 27 23 30
Gx:0^4 85 3 10 2
Gx:0Gcnot:0:1 45 1 4 50
"""
with open("tutorial_files/Example_Dataset.txt","w") as f:
f.write(dataset_txt)
ds = pygsti.io.load_dataset("tutorial_files/Example_Dataset.txt")
def compare(prefix, model, dataset):
chi2 = pygsti.tools.chi2(model, dataset)
logl = pygsti.tools.logl(model, dataset)
max_logl = pygsti.tools.logl_max(model, dataset)
k = dataset.get_degrees_of_freedom()
Nsigma = (2*(max_logl-logl) - k)/np.sqrt(2*k)
print(prefix, "\n chi^2 = ",chi2,"\n 2DeltaLogL = ", 2*(max_logl-logl),
"\n #std-deviations away from expected (%g) = " % k,Nsigma,"\n")
print("\nModel compared with:")
compare("1. Hand-chosen data (doesn't agree): ",mdl,ds)
compare("2. Model-generated data (agrees): ",mdl,ds_fake)
compare("3. Model-generated data w/no sample err (agrees *exactly*): ",mdl,ds_nosampleerr)
Loading tutorial_files/Example_Dataset.txt: 100% Model compared with: 1. Hand-chosen data (doesn't agree): chi^2 = 15507.570000000045 2DeltaLogL = 575.9735817498567 #std-deviations away from expected (15) = 102.41929496295496 2. Model-generated data (agrees): chi^2 = 7.160000000000015 2DeltaLogL = 7.508315185184529 #std-deviations away from expected (15) = -1.3677882555977823 3. Model-generated data w/no sample err (agrees *exactly*): chi^2 = 1.013809522725441e-28 2DeltaLogL = -9.094947017729282e-13 #std-deviations away from expected (15) = -2.7386127875259967
We can also look at these values on a per-circuit basis:
logl_percircuit = pygsti.tools.logl_terms(mdl, ds)
max_logl_percircuit = pygsti.tools.logl_max_terms(mdl, ds)
print("2DeltaLogL per circuit = ", 2*(max_logl_percircuit - logl_percircuit))
#ws.ColorBoxPlot('logl', pygsti.obj.LsGermsSerialStructure([0],)) TODO
2DeltaLogL per circuit = [ 0. 115.83041852 2.33389357 349.09795747 108.7113122 ]
It's possible to display model testing results within figure and HTML reports too. For more information on model testing, especially alongside GST, see the tutorial on model testing.
PyGSTi is able to perform two types of Randomized Benchmarking (RB). First, there is the standard Clifford-circuit-based RB protocol first defined by Magesan et al. Second, there is "Direct RB", which is particularly suited to multi-qubit benchmarking. More more details on using these protocols (e.g. how to generate a set of RB sequences) see the separate Clifford RB tutorial and Direct RB tutorial notebooks. Below we demonstrate how easy it is to analyze RB data (the same using either protocol). For further details and options on analyzing RB data, see the RB data analysis tutorial.
from pygsti.extras import rb
rbresults_txt = \
"""# Results from a DRB simulation
# Number of qubits
5
# RB length // Success counts // Total counts // Circuit depth // Circuit two-qubit gate count
0 37 50 50 25
10 36 50 51 33
20 38 50 74 56
30 35 50 78 62
50 28 50 102 87
100 25 50 139 163
200 20 50 241 288
400 10 50 446 546
"""
with open("tutorial_files/Example_RBData.txt","w") as f:
f.write(rbresults_txt)
rbdata = rb.io.import_rb_summary_data('tutorial_files/Example_RBData.txt')
rbresults = rb.analysis.std_practice_analysis(rbdata)
ws.RandomizedBenchmarkingPlot(rbresults)
Importing tutorial_files/Example_RBData.txt...Complete.
<pygsti.report.workspaceplots.RandomizedBenchmarkingPlot at 0x11f76c320>
The Robust Phase Estimation (RPE) protocol is designed to efficiently estimate a few specific parameters of certain single-qubit models. Below we demonstrate how to run RPE with the single-qubit model containing $X(\pi/2)$ and $Y(\pi/2)$ gates. The list of requisite circuits is given by make_rpe_angle_string_list_dict
and simulated noisy data is analyzed using analyze_rpe_data
. For more information on running RPE see the RPE tutorial.
from pygsti.extras import rpe
from pygsti.construction import std1Q_XY
import numpy as np
#Declare the particular RPE instance we are interested in (X and Y pi/2 rotations)
# Note: Prep and measurement are for the |0> state.
rpeconfig_inst = rpe.rpeconfig_GxPi2_GyPi2_00
stringListsRPE = rpe.rpeconstruction.make_rpe_angle_string_list_dict(10,rpeconfig_inst)
angleList = ['alpha','epsilon','theta']
numStrsD = {'RPE' : [6*i for i in np.arange(1,12)] }
#Create fake noisy model
mdl_real = std1Q_XY.target_model().randomize_with_unitary(.01,seed=0)
ds_rpe = pygsti.construction.generate_fake_data(mdl_real,stringListsRPE['totalStrList'],
1000,sampleError='binomial',seed=1)
#Run RPE protocol
resultsRPE = rpe.analyze_rpe_data(ds_rpe,mdl_real,stringListsRPE,rpeconfig_inst)
print('alpha_true - alpha_est_final =',resultsRPE['alphaErrorList'][-1])
print('epsilon_true - epsilon_est_final =',resultsRPE['epsilonErrorList'][-1])
print('theta_true - theta_est_final =',resultsRPE['thetaErrorList'][-1])
alpha_true - alpha_est_final = 4.136221313455479e-05 epsilon_true - epsilon_est_final = 4.843096931717028e-05 theta_true - theta_est_final = 0.012088866000938565
The DataComparator
object is designed to check multiple DataSet
objects for consistency. This procedure essentially answers the question: "Is it better to describe DataSet
s $A$ and $B$ as having been generated by the same set of probabilities or different sets?". This quick test is useful for detecting drift in experimental setups from one round of data-taking to the next, and doesn't require constructing any Model
objects. Below, we generate three DataSet
objects - two from the same underlying model and one from a different model - and show that we can detect this difference. For more information, see the tutorial on data set comparison.
from pygsti.construction import std1Q_XYI
#Generate data from two different models
mdlA = std1Q_XYI.target_model().randomize_with_unitary(.01,seed=0)
mdlB = std1Q_XYI.target_model().randomize_with_unitary(.01,seed=1)
circuits = pygsti.construction.make_lsgst_experiment_list(std1Q_XYI.target_model(),std1Q_XYI.fiducials,
std1Q_XYI.fiducials,std1Q_XYI.germs,[1,2,4,8])
#Generate the data for the two datasets, using the same model, and one with a different model
dsA1 = pygsti.construction.generate_fake_data(mdlA,circuits,100,'binomial',seed=10)
dsA2 = pygsti.construction.generate_fake_data(mdlA,circuits,100,'binomial',seed=20)
dsB = pygsti.construction.generate_fake_data(mdlB,circuits,100,'binomial',seed=30)
#Let's compare the two datasets.
print("Compare two *consistent* DataSets (generated from the same underlying model)")
comparator_A1_A2 = pygsti.objects.DataComparator([dsA1,dsA2])
comparator_A1_A2.implement(significance=0.05)
print("\nCompare two *inconsistent* DataSets (generated from different model)")
comparator_A1_B = pygsti.objects.DataComparator([dsA1,dsB])
comparator_A1_B.implement(significance=0.05)
#Plots of consistent (top) and inconsistent (bottom) cases
ws.DatasetComparisonHistogramPlot(comparator_A1_A2, log=True, display='pvalue', scale=0.8)
ws.DatasetComparisonHistogramPlot(comparator_A1_B, log=True, display='pvalue', scale=0.8)
Compare two *consistent* DataSets (generated from the same underlying model) Statistical hypothesis tests did NOT find inconsistency between the datasets at 5.00% significance. Compare two *inconsistent* DataSets (generated from different model) The datasets are INCONSISTENT at 5.00% significance. - Details: - The aggregate log-likelihood ratio test is significant at 17.93 standard deviations. - The aggregate log-likelihood ratio test standard deviations signficance threshold is 2.01 - The number of sequences with data that is inconsistent is 8 - The maximum SSTVD over all sequences is 0.42 - The maximum SSTVD was observed for Qubit * ---|Gx|-|Gx|-|Gx|-|Gi|-|Gi|-|Gi|-|Gi|-|Gi|-|Gi|-|Gi|-|Gi|-|Gy|-|Gy|-|Gy|---
<pygsti.report.workspaceplots.DatasetComparisonHistogramPlot at 0x1202f0cf8>
Gate set tomography (GST) is a protocol designed to solve the inverse of "use this model to simulate observed data"; its goal is to infer a model based on actual observed data. From a functional perspective, GST can be viewed as an inverse of the generate_fake_data
function we've used a bunch above: it takes a DataSet
and produces a Model
.
Because this inverse problem traverses some technical challenges, GST also requires a structured set of Circuits
to work reliably and efficiently. Here enters the concepts of "fiducial" and "germ" circuits, as well as a list of "maximum-repeated-germ-lengths" or just "max-lengths". For details, see the tutorial on the structure of GST circuits and the tutorial on fiducial and germ selection. The important takeaway is that the GST circuits are described below by the 4 variables: prep_fiducials
, meas_fiducials
, germs
, and maxLengths
.
Below, we generate a set of GST circuits and simulate them using a noisy (slightly depolarized) model of the some ideal operations to get a DataSet
.
construction is explained in a tutorial on We explain what these circuits are and how to build described above, that is, an inverse to . whatis the most complex protocol
from pygsti.construction import std1Q_XYI
# 1) get the target Model
mdl_ideal = std1Q_XYI.target_model()
# 2) get the building blocks needed to specify which circuits are needed
prep_fiducials, meas_fiducials = std1Q_XYI.prepStrs, std1Q_XYI.effectStrs
germs = std1Q_XYI.germs
maxLengths = [1,2,4] # roughly gives the length of the sequences used by GST
# 3) generate "fake" data from a depolarized version of mdl_ideal
mdl_true = mdl_ideal.depolarize(op_noise=0.01, spam_noise=0.001)
listOfExperiments = pygsti.construction.make_lsgst_experiment_list(
mdl_ideal, prep_fiducials, meas_fiducials, germs, maxLengths)
ds = pygsti.construction.generate_fake_data(mdl_true, listOfExperiments, nSamples=1000,
sampleError="binomial", seed=1234)
#Run GST
results = pygsti.do_stdpractice_gst(ds, mdl_ideal, prep_fiducials, meas_fiducials,
germs, maxLengths, modes="TP,Target", verbosity=1)
mdl_estimate = results.estimates['TP'].models['single']
print("2DeltaLogL(estimate, data): ", pygsti.tools.two_delta_logl(mdl_estimate, ds))
print("2DeltaLogL(true, data): ", pygsti.tools.two_delta_logl(mdl_true, ds))
print("2DeltaLogL(ideal, data): ", pygsti.tools.two_delta_logl(mdl_ideal, ds))
-- Std Practice: [##################################################] 100.0% (Target) -- 2DeltaLogL(estimate, data): 443.55901585984975 2DeltaLogL(true, data): 476.6246206732467 2DeltaLogL(ideal, data): 103821.37412338122
Recall that a lower $2\Delta\log\mathcal{L}$ means a better agreement between model and data. Note that the GST estimate fits the data best (even slightly better than the true model, because it agrees better with the finite sample noise), and the GST estimate is much better than the ideal model.
GST is essentially an automated model-tester that keeps modifying and testing models until it finds one the agrees with the data as well any model of the specified type can. To learn more about how to run GST see the GST overview tutorial and the GST methods tutorial.
The output of GST is an entire Model
(contrasted with the one or several numbers of RB and RPE), there are many ways to assess and understand the performance of a QIP based on GST results. The Results
object in pyGSTi is responsible for holding the GST and other model-based protocol results. The structure and use of a Results
object is explained in the Results tutorial. A common use for Results
objects is to generate "reports". PyGSTi has the ability to generate HTML reports (a directory of files) whose goal is to display relevant model vs. data metrics such as $2\Delta\log\mathcal{L}$ as well as model vs. model metrics like process fidelity and diamond distance. To learn more about generating these "model-explaining" reports see the report generation tutorial.
Here's an example of how to generate a report (it will auto-open in a new tab; if it doesn't display try it in FireFox):
pygsti.report.create_standard_report(results, filename="tutorial_files/myFirstGSTReport",
title="Example GST Report", auto_open=True, verbosity=1)
*** Creating workspace *** *** Generating switchboard *** Found standard clifford compilation from std1Q_XYI Found standard clifford compilation from std1Q_XYI *** Generating tables ***
/Users/enielse/research/pyGSTi/packages/pygsti/report/factory.py:785: UserWarning: Idle tomography failed: Label{layers}
*** Generating plots *** *** Merging into template file *** Output written to tutorial_files/myFirstGSTReport directory Opening tutorial_files/myFirstGSTReport/main.html... *** Report Generation Complete! Total time 21.9227s ***
<pygsti.report.workspace.Workspace at 0x120573630>
pyGSTi contains all the ingredients necessary for multi-qubit tomography, and this functionality is currently in the beta-testing stages. Stay tuned!
#Coming soon
This concludes our overview of what can be done with pyGSTi. There are a few minor topics that haven't been touched on that we've collected within the next tutorial on miscellaneous topics, so you might want to take a quick look there to see if there's anything you're especially interested in.