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.

- computing circuit outcome probabilities
- simulating observed data based on a model
- testing how well a model describes a set of data
- running Randomized Benchmarking (RB)
- running Robust Phase Estimation (RPE)
- performing data set comparison tests
- running Gate Set Tomography (GST)
- idle tomography (IDT)
- drift characterization
- time-dependent GST
- multi-qubit tomography

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 [1]:

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

In [2]:

```
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`
```

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.

In [3]:

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

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 [4]:

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

```
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
```

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.

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 [6]:

```
#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)
```

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

```
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,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])
```

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.

In [8]:

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

```
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, 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['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 [10]:

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

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.

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.

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.

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

In [11]:

```
#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.