{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Different ways to run Gate Set Tomography"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `pygsti` package provides multiple ways to use its core Gate Set Tomography (GST) algorithms.  This  tutorial will show you how to work with pyGSTi's GST protocol objects to perform GST in different ways with a minimial amount of effort.  In order to run the GST protocol there are 3 essential ingredients: 1) an \"experiment design\" specifying the structure of the GST circuits and how the data should be collected, 2) the outcome counts for the circuits specified by the experiment design, and 3) a desired, or \"target\", `Model`.  The [GST overview tutorial](GST-Overview.ipynb), gave an end-to-end example of how to construct a GST experiment design, run GST, and generate a report.  This tutorial focuses on the first and second steps in more detail; related information about circuit construction and report generation can be found in the [GST circuits tutorial](../objects/advanced/GSTCircuitConstruction.ipynb) and [report generation tutorial](../reporting/ReportGeneration.ipynb).\n",
    "\n",
    "There are two different `Protocol` objects within pyGSTi for running GST:\n",
    "\n",
    "- [`GateSetTomography`](#gatesettomography) - runs a single model optimization based on a *given* initial model that can have any parameterization you like.  This protocol can be run on any `GateSetTomographyDesign` experiment design, which only needs a target model (to describe what gates occur in the circuits) and a list of circuit lists to specify the circuits used for each iteration of the model optimization.\n",
    "\n",
    "- [`StandardGST`](#standardgst) - runs multiple model optimizations based on an `ExplicitOpModel` target model by parameterizing this model in different ways.  The target model is expected to be a part of the experiment design, and only `StandardGSTDesign`-type experiment designs are allowed since the usual germs-and-fiducials structure of the GST circuits is expected.\n",
    "\n",
    "Overall, the `GateSetTomography` protocol is more flexible than the `StandardGST` protocol, but requires a little more work to get going because its inputs are more complicated.  Both protocols return a `ModelEstimateResults` object when they are run."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import print_function\n",
    "import pygsti"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Setup\n",
    "In the [DataSet tutorial](../objects/DataSet.ipynb) we simulate the circuits required by a GST experiment design and save the results.  In this tutorial, we'll be analyzing that data.  This illustrates a typical workflow where at some earlier time you setup an experiment (a \"GST experiment in this case) and save the experiment design to disk and at some later time (after the data has been collected) you want to analyze it.  Now *is* that later time, and we start by reading the the data we've collected."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = pygsti.io.load_data_from_dir(\"../tutorial_files/Example_GST_Data\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"gatesettomography\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## `GateSetTomography`\n",
    "This protocol performs a single model optimization, and so computes a **single GST estimate** given a `DataSet`, a target `Model`, and other parameters.  (The returned `ModelEstimateResults` object may sometimes contain multiple related estimates in certain cases, but in these cases all the estimates are closely related.)  The experiment design provides all of the information about the GST circuits, in this case a *standard*  (*prep_fiducial + germ^power + meas_fiducial*) set, so the only thing needed by the protocol is an initial `Model` to optimize.  Thus, the `GateSetTomography` protocol is essentially just a model optimizer that you give an initial point.  Importantly, this initial point (a `Model`) also specifies the *parameterization*, i.e. the space of parameters that are optimized over.\n",
    "\n",
    "Minimally, when using `GateSetTomography` you should set the parameterization of the initial model.  This can be viewed as setting the constraints on the optimization.  For instance, when the gates in the model are parameterized as trace-preserving (TP) maps, the optimization will be constrained to trying gate sets with TP gates (because every set of parameters corresponds to a set of TP gates).  In the cell below, we constrain the optimization to TP gate sets by using `.target_model(\"TP\")`, which returns a version of the target model where all the gates are TP-parameterized, the state preparation has trace = 1, and the POVM effects always add to the identity.  This could also be done by calling `set_all_parameterizations(\"TP\")` on the fully-parameterized target model returned by `.target_model()`.  See the [tutorial on explicit models](../objects/ExplicitModel.ipynb) for more information on setting a model's parameterization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from pygsti.modelpacks import smq1Q_XYI\n",
    "target_model_TP = smq1Q_XYI.target_model(\"TP\")\n",
    "proto = pygsti.protocols.GateSetTomography(target_model_TP)\n",
    "results_TP = proto.run(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A summary of what's inside a Results object is obtained by printing it\n",
    "(for more examples of how to use a Results object, see the [Results tutorial](../objects/advanced/Results.ipynb))."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "print(results_TP)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Gauge optimization parameters\n",
    "The `gaugeopt_suite` argument specifies a set of gauge optimizations to be performed on the final GST estimate.  It is a dictionary whose keys are gauge-optimization names (these can be whatever you want) and whose values are dictionaries of arguments ultimately to be passed to the `gaugeopt_to_target` function (which provides full documentation).  (For example, by specifying `itemWeights` we can set the ratio of the state preparation and measurement (SPAM) weighting to the gate weighting when performing a gauge optimization.)  In lieu of a dictionary of `gaugeopt_to_target` arguments, the elements of `gaugeopt_suite` may also be strings which name a built-in set of gauge optimizations (e.g. `\"stdgaugeopt\"` is the name of the standard gauge optimization).\n",
    "\n",
    "If `gaugeopt_suite` is set to a string, this is the same as passing a dictionary with a single key-value pair where both key and value are equal to the string.  Thus, the default `\"stdgaugeopt\"` is equivalent to specifying the dictionary `{\"stdgaugeopt\": \"stdgagueopt\"}`.\n",
    "\n",
    "The example below performs a customized gauge-optimization where the gate parameters are weighted 1000 times more relative to the SPAM parameters.  Mathematically this corresponds to a multiplicative factor of 0.001 preceding the sum-of-squared-difference terms corresponding to SPAM elements in the model.   Typically it is good to weight the gates parameters more heavily since GST amplifies gate parameter errors via long operation sequences but cannot amplify SPAM parameter errors.  For more details on the arguments of `gaugeopt_to_target`, see the previous tutorial on low-level algorithms.  For more infomation, see the [gauge optimization tutorial](advanced/GaugeOpt.ipynb).\n",
    "\n",
    "The cell below also illustrates how you can create a TP target model by calling `set_all_parameterizations` explicitly instead of using the equivalent and more condensed `.target_model(\"TP\")`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "target_model_TP2 = smq1Q_XYI.target_model() # a \"fully parameterized\" (unconstrained) model\n",
    "target_model_TP2.set_all_parameterizations(\"TP\") # change parameterization to TP gates\n",
    "\n",
    "proto = pygsti.protocols.GateSetTomography(\n",
    "    target_model_TP2, name=\"GSTwithMyGO\",\n",
    "    gaugeopt_suite={'my_gauge_opt': {'itemWeights': {'gates': 1.0, 'spam': 0.001}}}\n",
    "    )\n",
    "results_TP2 = proto.run(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(results_TP2.estimates['default'].goparameters.keys())  # names of all the gauge opts that were done\n",
    "custom_gauge_opt_model = results_TP2.estimates['default'].models['my_gauge_opt']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### running GST using a custom set of circuits\n",
    "So far we've giving the `GateSetTomography.run` method an \"standard\" experiment design containing circuits chosen to amplify all of a standard TP (or CPTP) model's parameters (see the `StandardGSTExpermentDesign` used in the [DataSet tutorial](../objects/DataSet.ipynb)).  A `GateSetTomography` protocol can be run on more general experiment designs, namely those that specify the circuits to use as either a list of lists of `Circuit` objects or a list of or single `CircuitStructure` object(s).  A `CircuitStructure` is preferable as it allows the structured plotting of the sequences in report figures.  In this example, we'll just generate a standard set of circuit structures, but with some of the sequences randomly dropped (see the [tutorial on GST circuit reduction](advanced/GST-FiducialPairReduction.ipynb)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Create the same sequences but drop 50% of them randomly for each repeated-germ block.\n",
    "# and only go out to a max-length of 8\n",
    "orig_design = data.edesign  # the original StandardGSTDesign\n",
    "custom_maxlengths = [1, 2, 4, 8]  # a subset of orig_design.maxlengths\n",
    "circuit_structs = pygsti.construction.make_lsgst_structs(\n",
    "    orig_design.target_model, orig_design.prep_fiducials, orig_design.meas_fiducials,\n",
    "    orig_design.germs, custom_maxlengths, keepFraction=0.5, keepSeed=2020)\n",
    "reduced_exp_design = pygsti.protocols.StructuredGSTDesign(orig_design.target_model, circuit_structs)\n",
    "reduced_data = pygsti.protocols.ProtocolData(reduced_exp_design, data.dataset)\n",
    "\n",
    "\n",
    "proto = pygsti.protocols.GateSetTomography(target_model_TP2, name=\"GSTwithReducedData\")\n",
    "results_reduced = proto.run(reduced_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"standardgst\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## `StandardGST`\n",
    "The protocol embodies a standard *set* of GST protocols to be run on a set of data.  It essentially runs multiple `GateSetTomography` protocols on the given data which use different parameterizations of an `ExplicitOpModel`  (the `StandardGST` protocol doesn't work with other types of `Model` objects, e.g. *implicit* models, which don't implement `set_all_parameterizations`).  The `modes` argument is a comma-separated list of the parameterization types that should be run (e.g. `\"TP,CPTP\"` will compute a Trace-Preserving estimate *and* a Completely-Positive & Trace-Preserving estimate). The currently available modes are:\n",
    " - \"full\" : unconstrained gates (fully parameterized)                                                                 \n",
    " - \"TP\"   : TP-constrained gates and state preparations\n",
    " - \"CPTP\" : CPTP-constrained gates and TP-constrained state preparations               \n",
    " - \"H+S\"  : Only Hamiltonian and Pauli stochastic errors allowed (CPTP)                                             \n",
    " - \"S\"    : Only Pauli-stochastic errors allowed (CPTP)                                                           \n",
    " - \"Target\" : use the target (ideal) gates as the estimate     \n",
    "\n",
    "Gauge optimization(s) are controlled by the `gaugeopt_suite` and `gaugeopt_target` arguments, jsut as in `GateSetTomography`.  The `gaugeopt_target` argument may be set to a `Model` that is used as the target for gauge optimization, overriding the (typically ideal) target gates given by within the experiment design."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "results_stdprac = pygsti.protocols.StandardGST().run(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Estimates: \", \", \".join(results_stdprac.estimates.keys()))\n",
    "print(\"TP Estimate's gauge optimized models: \", \", \".join(results_stdprac.estimates[\"TP\"].goparameters.keys()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we'll perform the same analysis but with a **non-default standard suite of gauge optimizations** - this one toggles the SPAM penalty in addition to varying the spam weight (the default suite just varies the spam weight without any SPAM penalty).  See the [gauge optimization tutorial](advanced/GaugeOpt.ipynb) for more details on gauge optmization \"suites\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "proto = pygsti.protocols.StandardGST(gaugeopt_suite=\"varySpam\", name=\"StdGST_varySpam\")\n",
    "results_stdprac_nondefaultgo = proto.run(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Estimates: \", \", \".join(results_stdprac_nondefaultgo.estimates.keys()))\n",
    "print(\"TP Estimate's gauge optimized models: \", \", \".join(results_stdprac_nondefaultgo.estimates[\"TP\"].goparameters.keys()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we'll demonstrate how to specify a fully custom set of gauge optimization parameters and how to use a **separately-specified target model for gauge optimization**.  You can get a more intuitive gauge-optimized `Model` when by placing as much expected noise as possible into the gauge-optimization target, as this essentially tells the algorithm \"this is what I think the estimated model should look like\".  If you just use the perfect or ideal model for this (the default), then the gauge optimizer may make tradeoffs which don't reflect the expected physics (remember, all gauge-equivalent models product the same observables!).  For example, it may spread error across all your gate operations when you expect just the 2-qubit operations are noisy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_goparams = { 'itemWeights': {'gates': 1.0, 'spam': 0.001} }\n",
    "my_gaugeOptTarget = smq1Q_XYI.target_model().depolarize(\n",
    "    op_noise=0.005, spam_noise=0.01) # a guess at what estimate should be\n",
    "\n",
    "proto = pygsti.protocols.StandardGST(gaugeopt_suite={ 'myGO': my_goparams },\n",
    "                                     gaugeopt_target=my_gaugeOptTarget,\n",
    "                                     name=\"StdGST_myGO\")\n",
    "results_stdprac_nondefaultgo = proto.run(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Estimates: \", \", \".join(results_stdprac_nondefaultgo.estimates.keys()))\n",
    "print(\"TP Estimate's gauge optimized models: \", \", \".join(results_stdprac_nondefaultgo.estimates[\"TP\"].goparameters.keys()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To finish up, we'll write the results for processing in other tutorials.  We do this by calling `.write` on the results objects, optionally specifying the root diretory under which the results should be written.  This is the *same* root directory that the experiment design and data are written to, as subdirectories beneath this directory separate these quantities.\n",
    "\n",
    "Two remarks are in order:\n",
    "1. When results are from running a protocol on data that was loaded with the `load_data_from_dir` method (see the beginning of this notebook), then knowledge of this directory is remembered and you don't need to give a directory to `write` (this is the case for all except `results_reduced`, which created a new experiment design containing less experiments).\n",
    "\n",
    "2. Notice how the `name=` arguments given to protocols above are used as sub-directory names, e.g. under the \"tutorial_files/Example_GST_Data/results\" parent directory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "results_TP.write()  # uses \"../tutorial_files/Example_GST_Data\" (where data was loaded from)\n",
    "results_TP2.write() # ditto\n",
    "results_stdprac.write() # ditto\n",
    "results_reduced.write(\"../tutorial_files/Example_Reduced_GST_Data\") # choose a different dir"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "While it is possible to **pickle** a results object, this method of serialization is **not recommended** for long-term storage since pickles are relatively fragile to changes in pyGSTi or other python libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#Not recommended:\n",
    "# import pickle\n",
    "# pickle.dump(results_TP, open('../tutorial_files/exampleResults_TP.pkl',\"wb\"))\n",
    "# pickle.dump(results_TP2, open('../tutorial_files/exampleResults_TP2.pkl',\"wb\"))\n",
    "# pickle.dump(results_reduced, open('../tutorial_files/exampleResults_reduced.pkl',\"wb\"))\n",
    "# pickle.dump(results_stdprac, open('../tutorial_files/exampleResults_stdprac.pkl',\"wb\"))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}