Things you can do with pyGSTi's "essential objects"

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.

Contents

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).

In [ ]:
import pygsti
import numpy as np
ws = pygsti.report.Workspace()
ws.init_notebook_mode(autodisplay=True)

Computing circuit outcome probabilities

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. For more information on circuit simulation, see the circuit simulation tutorial.

In [ ]:
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`

Simulating observed data based on a model

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 sample_error='none', as demonstrated below.

In [ ]:
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, n_samples=100,
                                                 sample_error='multinomial', seed=1234)
print("Normal:")
print(ds_fake)

ds_nosampleerr = pygsti.construction.generate_fake_data(mdl, circuit_list, n_samples=100,
                                                 sample_error='none')
print("Without any sample error:")
print(ds_nosampleerr)

Testing how well a model describes a set of data

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:

  • the larger $\log\mathcal{L}$ is, and the smaller $\chi^2$ is, the better the model agrees with the data.
  • the value of $\log\mathcal{L}$ doesn't mean anything in isolation - only when compared to other $\log\mathcal{L}$ values.
  • one can compute the $\log\mathcal{L}$ of a "maximal model" that agrees with the data exactly. We call this value $\log\mathcal{L}_{max}$.
  • in the limit of many samples, $\chi^2 \approx 2(\log\mathcal{L}_{max}-\log\mathcal{L})$, and we denote the latter $2\Delta\log\mathcal{L}$.
  • let $k$ be the the number of independent degrees of freedom in the data (e.g. each row of 4 (00, 01, 10, and 11) counts in the example below contributes $3$ degrees of freedom because the counts are constrained to add to 100, so $k=5*3=15$). If the model is "valid" (i.e. it could have generated the data) then $2\Delta\log\mathcal{L}$ should have come from a $\chi^2_k$ distribution, i.e. it has expectation value $k$ and standard deviation $\sqrt{2k}$.

Here's how we compute the $\chi^2$ and $2\Delta\log\mathcal{L}$ between some data and a model:

In [ ]:
import numpy as np
dataset_txt = \
"""## Columns = 00 count, 01 count, 10 count, 11 count
{}@(0,1)            100   0   0   0
Gx:[email protected](0,1)           55   5  40   0
Gx:0Gy:[email protected](0,1)       20  27  23  30
Gx:0^[email protected](0,1)         85   3  10   2
Gx:0Gcnot:0:[email protected](0,1)  45   1   4  50
"""
with open("tutorial_files/Example_Short_Dataset.txt","w") as f:
    f.write(dataset_txt)
ds = pygsti.io.load_dataset("tutorial_files/Example_Short_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)

We can also look at these values on a per-circuit basis:

In [ ]:
logl_percircuit = pygsti.tools.logl_per_circuit(mdl, ds)
max_logl_percircuit = pygsti.tools.logl_max_per_circuit(mdl, ds)
print("2DeltaLogL per circuit = ", 2*(max_logl_percircuit - logl_percircuit))

#ws.ColorBoxPlot('logl', pygsti.obj.LsGermsSerialStructure([0],)) TODO

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(using protocol objects) and the functions for model testing.

Randomized Benchmarking (RB)

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. If you're interested using pyGSTi noise models to simulate RB results, see the tutorials for Clifford RB simulation using explicit-op models and Clifford RB simulation using implicit models.

In [ ]:
#TODO REMOVE
# 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)

Robust Phase Estimation (RPE)

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.

In [ ]:
from pygsti.extras import rpe
from pygsti.modelpacks.legacy 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
print(stringListsRPE['totalStrList'][0])
mdl_real = std1Q_XY.target_model().randomize_with_unitary(.01,seed=0)
ds_rpe = pygsti.construction.generate_fake_data(mdl_real,stringListsRPE['totalStrList'],
                                            1000,sample_error='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])

Data set comparison tests

The DataComparator object is designed to check multiple DataSet objects for consistency. This procedure essentially answers the question: "Is it better to describe DataSets $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.

In [ ]:
from pygsti.modelpacks import smq1Q_XYI

#Generate data from two different models
mdlA = smq1Q_XYI.target_model().randomize_with_unitary(.01,seed=0)
mdlB = smq1Q_XYI.target_model().randomize_with_unitary(.01,seed=1)

circuits = pygsti.construction.make_lsgst_experiment_list(
    smq1Q_XYI.target_model(),smq1Q_XYI.prep_fiducials(),
    smq1Q_XYI.meas_fiducials(),smq1Q_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)

Gate Set Tomography (GST)

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.

In [ ]:
from pygsti.modelpacks import smq1Q_XYI

# 1) get the target Model
mdl_ideal = smq1Q_XYI.target_model()

# 2) get the building blocks needed to specify which circuits are needed
prep_fiducials, meas_fiducials = smq1Q_XYI.prep_fiducials(), smq1Q_XYI.meas_fiducials()
germs = smq1Q_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, n_samples=1000,
                                            sample_error="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['stdgaugeopt']
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))

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 using functions that act on essential objects see the function-based GST overview tutorial and the GST driver functions tutorial. GST can also be run in a more object-oriented way, using Protocol objects as described in the GST overview 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 ModelEstimateResults object in pyGSTi is responsible for holding the GST and other model-based protocol results. The structure and use of a ModelEstimateResults 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):

In [ ]:
pygsti.report.construct_standard_report(
    results, title="Example GST Report", verbosity=1
).write_html("tutorial_files/myFirstGSTReport", auto_open=True, verbosity=1)

Idle tomography

Idle tomography estimates the error rates of an $n$-qubit idle operation using relatively few sequences. To learn more about how to use it, see the idle tomography tutorial.

Drift Characterization

Time-series data can be analyzed for significant indications of drift (time variance in circuit outcome probabilities). See the tutorial on drift characterization for more details.

Time-dependent gate set tomography

pyGSTi has recently added support for time-dependent models and data sets, allowing the GST to be performed in a time-dependent fashion. See the time-dependent GST tutorial for more details.

Multi-qubit tomography

pyGSTi contains all the ingredients necessary for multi-qubit tomography, and this functionality is currently in the beta-testing stages. Stay tuned!

In [ ]:
#Coming soon

What's next?

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.