---
jupytext:
  text_representation:
    extension: .md
    format_name: myst
    format_version: 0.13
    jupytext_version: 1.17.3
kernelspec:
  display_name: Python 3 (ipykernel)
  language: python
  name: python3
---

# Low-level Algorithm API
Once we have data for GST, there are several *algorithms* we can run on it to produce tomographic estimates. Depending on the amount of data you have, and time available for running Gate Set Tomography, one algorithm may be preferable over the others.  **What is typically thought of as "standard GST" is the iterative maximum-likelihood optimization implemented by `run_iterative_gst`** which optimizes one or more objective functions iteratively, meaning it daisy-chains multiple optimizations, changing (typically *adding* to) the circuits under consideration or the objective function itself (e.g. from $\chi^2$ to the log-likelihood).

`pygsti` contains two main low-level algorithms that form the base of the GST protocol(s):

* **Linear gate set tomography (LGST)**: Uses short operation sequences to quickly compute a rough (low accuracy) estimate of a model by linear inversion.  The model estimate produced has a very specific form: it holds the dense process matrices of each "gate" (or circuit layer) under consideration.

* **Long-sequence gate set tomography (LSGST)**: Minimizes an objective function that quantifies a badness-of-fit between a model and data set.  Here "sequence" is synonymous with "circuit".  The optimization is performed in stages, using successively larger sets of longer and longer circuits.

If you're curious, the implementation of these algorithms are found in the `pygsti.algorithms.core` module.  In this tutorial, we'll show how to invoke each of these low-level algorithms directly.  This is probably only useful when building your own protocol that needs to borrow pieces of GST without borrowing the entire protocol.

Additionally, `pygsti` contains **gauge-optimization** algorithms.  Because the outcome data (the input to the GST algorithms above) only determines a model up to some number of un-physical "gauge" degrees of freedom, it is often desirable to optimize the `Model` estimate obtained from a GST algorithm within the space of its gauge freedoms.  This process is called "gauge-optimization" and the final part of this tutorial demonstrates how to gauge-optimize a model using various criteria.

##  Low-level GST Algorithms
The ingredients needed by the LGST and LSGST algorithms are:
- a "target" `Model` which defines the desired gates.  This model is used by LGST to specify the various gate, state preparation, POVM effect, and SPAM labels, as well as to provide an initial guess for the *gauge* degrees of freedom.
- a `DataSet` containing the data that GST attempts to fit using the probabilities generated by a single `Model`.  This data set must at least contain the data for the circuits required by the algorithm that is chosen.
- for LSGST, a list-of-lists of `Circuit` objects that specify which circuits are used during each iteration of the algorithm (the length of the top-level list defines the number of interations).

```{code-cell} ipython3
import pygsti
import json
```

```{code-cell} ipython3
#We'll use the standard I, X(pi/2), Y(pi/2) model that we generated data for in the DataSet tutorial
from pygsti.modelpacks import smq1Q_XYI

target_model = smq1Q_XYI.target_model()
prep_fiducials = smq1Q_XYI.prep_fiducials()
meas_fiducials = smq1Q_XYI.meas_fiducials()
germs = smq1Q_XYI.germs()

maxLengthList = [1,2,4,8,16] #for use in iterative algorithms

ds = pygsti.io.load_dataset("../../tutorial_files/Example_Dataset.txt", cache=True)
dsLowCounts = pygsti.io.load_dataset("../../tutorial_files/Example_Dataset_LowCnts.txt", cache=True)

depol_model = target_model.depolarize(op_noise=0.1)

print("Loaded target model with operation labels: ", target_model.operations.keys())
print("Loaded fiducial lists of lengths: ", (len(prep_fiducials),len(meas_fiducials)))
print("Loaded dataset of length: ", len(ds))
```

### Using LGST to get an initial estimate

An important and distinguising property of the LGST algorithm is that it does *not* require an initial-guess `Model` as an input.  It uses linear inversion and short sequences to obtain a rough model estimate.  As such, it is very common to use the LGST estimate as the initial-guess starting point for more advanced forms of GST.

```{code-cell} ipython3
#Run LGST to get an initial estimate for the gates in target_model based on the data in ds

#run LGST
mdl_lgst = pygsti.run_lgst(ds, prep_fiducials, meas_fiducials, target_model=target_model, verbosity=1)

#Gauge optimize the result to match the target model
mdl_lgst_after_gauge_opt = pygsti.gaugeopt_to_target(mdl_lgst, target_model, verbosity=1)

#Contract the result to CPTP, guaranteeing that the gates are CPTP
mdl_clgst = pygsti.contract(mdl_lgst_after_gauge_opt, "CPTP", verbosity=1)
```

```{code-cell} ipython3
print(mdl_lgst)
```

### Long-sequence GST
Long-sequence GST operates by repeatedly fitting a model to a data set and using the result to seed a similar optimization that compares the data of more circuits.  As such a list-of-lists of circuits, one per iteration, is required.  The elements of these lists are typically repetitions of short "germ" circuits such that the final circuit does not exceed some maximum length (depth).  We use the `make_lsgst_lists` function to create such lists below.

```{code-cell} ipython3
#Get lists of operation sequences for successive iterations of LSGST to use
lsgstListOfLists = pygsti.circuits.create_lsgst_circuit_lists(target_model, prep_fiducials, meas_fiducials, germs, maxLengthList)
```

A single iteration, whereby the parameters of a `Model` are adjusted to find the best fit to a `DataSet`, is performed by the `run_gst_fit` function.  The "best fit" is defined as the model parameters which minimize some objective function.  `run_gst_fit` takes as arguments an a `ModelDatasetCircuitStore`, which simply bundles together a `Model`, `Dataset`, and set of `Circuits`, as well as an objective function builder (an object that encapsulates an objective function but lacks some needed context, and can *build* an objective function given this context). For convience, the `run_gst_fit_simple` function abstracts away the `ModelDatasetCircuitStore` construction.

For example, the cell below finds the model parameters which minimize the log-likelihood computed over the first list of circuits.

```{code-cell} ipython3
opt_result, mdl_single_optimization = \
    pygsti.algorithms.run_gst_fit_simple(ds, mdl_clgst, lsgstListOfLists[0],
                                         optimizer=None, objective_function_builder="logl",
                                         resource_alloc=None, verbosity=1)
```

The iterative algorithm, which essentially runs `run_gst_fit` multiple times, is implemented by `run_iterative_gst`.  The number of iterations is set by the number of circuit lists provided as the 3rd argument.  It takes two lists of objective function builders: the first, `iteration_objfn_builders` gives the objective functions to (sequentially) optimize on each iteration, while the second, `final_objfn_builders` gives additional objective functions to optimize on the final iteration.

```{code-cell} ipython3
models, optimums, final_objfn, mdc_store = \
    pygsti.algorithms.run_iterative_gst(ds, mdl_clgst, lsgstListOfLists,
                                       optimizer={'tol': 1e-5}, resource_alloc=None, verbosity=2,
                                       iteration_objfn_builders=['chi2'],
                                       final_objfn_builders=['logl'])
mdl_lsgst = models[-1]  # the model from the final iteration, i.e. the "final GST model"
```

## Gauge Optimization
All gauge optimization algorithms perform essentially the same task - to find the model which optimizes some objective function from within the set or space of models that are gauge-equivalent to some starting set.  This is accomplished in `pygsti` using the following mechanism:
- one begins with an initial `ExplicitOpModel`, call it `mdl`, to be gauge-optimized.
- a `pygsti.objects.GaugeGroup` instance defines a parameterized group of allowable gauge transformations.  This gauge group must be compatible with the `mdl`'s parameterization, so that `mdl.transform_inplace` (which calls `LinearOperator.transform_inplace` and `SPAMVec.transform_inplace`) is able to process elements of the `GaugeGroup` (obtained via a call to `GaugeGroup.element(params)`).  That is, the gauge transformation must map between models with the *same* parameterization (that give by `mdl`).  Because of the close interplay between a model's parameterization and its allowed gauge transformations, `Model` objects can contain a `GaugeGroup` instance as their `default_gauge_group` member. In many circumstances, `mdl.default_gauge_group` is set to the correct gauge group to use for a given `Model`.
- `pygsti.gaugeopt_custom(...)` takes an intial `ExplicitOpModel`, an objective function, and a `GaugeGroup` (along with other optimization parameters) and returns a gauge-optimized `ExplicitOpModel`.  Note that if its `gauge_group` argument is left as `None`, then the model's default gauge group is used.  And objective function which takes a single `ExplicitOpModel` argument and returns a float can be supplied, giving the user a fair amount of flexiblity.
- since usually the objective function is one which compares the model being optimized to a fixed "target" model, `pygsti.gaugeopt_to_target(...)` is a routine able to perform these common types of gauge optimization.  Instead of an objective function, `gaugeopt_to_target` takes a target `ExplicitOpModel` and additional arguments (see below) from which it constructs a objective function and then calls `gaugeopt_custom`.  It is essetially a convenience routine for constructing common gauge optimization objective functions.  Relevant arguments which affect what objective function is used are:
  - `target_model` : the `ExplicitOpModel` to compare against - i.e., the one you want to gauge optimize toward. **Note that this doesn't have to be a set of ideal gates **- it can be any (imperfect) model that reflects your expectations about what the estimates should look like.
  - `item_weights` : a dictionary of weights allowing different gates and/or SPAM operations to be weighted differently when computing the objective function's value.
  - `CPpenalty` : a prefactor multiplying the sum of all the negative Choi-matrix eigenvalues corresponding to each of the gates.
  - `TPpenalty` : a prefactor multiplying the sum of absoulte-value differences between the first row of each operation matrix and `[1 0 ... 0 ]` and the discrpance between the first element of each state preparation vector and its expected value.
  - `validSpamPenalty` : a prefactor multiplying penalty terms enforcing the non-negativity of state preparation eigenavlues and that POVM effect eigenvalues lie between 0 and 1.
  - `gates_metric` : how to compare corresponding gates in the gauge-optimized and target sets. `"frobenius`" uses the frobenius norm (weighted before taking the final sqrt), `"fidelity"` uses the *squared process infidelity* (squared to avoid negative-infidelity issues in non-TP models), and `"tracedist"` uses the trace distance (weighted after computing the trace distance between corresponding gates).
  - `spamMetric` : how to compare corresponding SPAM vectors. `"frobenius"` (the default) should be used here, as `"fidelity"` and `"tracedist"` compare the "SPAM gates" -- the outer product of state prep and POVM effect vectors -- which isn't a meaningful metric.
  
The cell below demonstrates some of common usages of `gaugeopt_to_target`.

```{code-cell} ipython3
mdl = mdl_lsgst.copy() #we'll use the MLGST result from above as an example
mdl_go1 = pygsti.gaugeopt_to_target(mdl, target_model) # optimization to the perfect target gates
mdl_go2 = pygsti.gaugeopt_to_target(mdl, depol_model) # optimization to a "guess" at what the estimate should be
mdl_go3 = pygsti.gaugeopt_to_target(mdl, target_model, {'gates': 1.0, 'spam': 0.01}) 
  # weight the gates differently from the SPAM operations
mdl_go4 = pygsti.gaugeopt_to_target(mdl, target_model, {'gates': 1.0, 'spam': 0.01, 'Gx': 10.0, 'E0': 0.001}) 
  # weight an individual gate/SPAM separately (note the specific gate/SPAM labels always override
  # the more general 'gates' and 'spam' weight values). 
mdl_go5 = pygsti.gaugeopt_to_target(mdl, target_model, gates_metric="tracedist") #use trace distance instead of frobenius

from pygsti.models.gaugegroup import UnitaryGaugeGroup
print("Default gauge group = ",type(mdl.default_gauge_group)) # default is FullGaugeGroup
mdl_go6 = pygsti.gaugeopt_to_target(mdl, target_model, gauge_group=UnitaryGaugeGroup(1, 'pp'))
  #gauge optimize only over unitary gauge transformations

print("\ngaugeopt_to_target output:")
mdl_go7 = pygsti.gaugeopt_to_target(mdl, target_model, verbosity=3) # show output
print("Final frobenius distance between mdl_go7 and target_model = ", mdl_go7.frobeniusdist(target_model))
```