{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Instability Detection and Characterization\n",
    "\n",
    "This tutorial shows how to implement instability (\"drift\") detection and characterization on time-stamped data. This data can be from *any* quantum circuits, on *any* number of qubits, but we require around 100+ time-stamps per circuit (perhaps fewer if there are multiple measurement outcomes per time-stamp). If you only have data that is binned into a few different time periods then consider instead using the `DataComparator` object demonstrated in the [DataSetComparison](../algorithms/DatasetComparison.ipynb) tutorial.\n",
    "\n",
    "Currently the gap between data collection times for each circuit is required to be approximately constant, both across the data collection times for each circuit, and across circuits. If this is not the case the code should still work, but the analysis it performs may be significantly sub-optimal, and interpretting the results is more complicated. There is beta-level capabilities within the functions used below to properly analyze unequally-spaced data, but it is untested and will not be used with the default options in the analysis code. This limitation will be addressed in a future release of pyGSTi.\n",
    "\n",
    "This notebook is an introduction to these tools, and it will be augmented with further notebooks at a later date."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Importing the drift module is essential\n",
    "from pygsti.extras import drift\n",
    "\n",
    "# Importing all of pyGSTi is optional, but often useful.\n",
    "import pygsti"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quick and Easy Analysis\n",
    "First we import some *time-stamped* data. For more information on the mechanics of using time-stamped `DataSets` see the [TimestampedDataSets](../objects/advanced/TimestampedDataSets.ipynb) tutorial. The data we are importing is from long-sequence GST on $G_i$, $G_x$, and $G_y$ with time-dependent coherent errors on the gates.\n",
    "\n",
    "We load the time-dependent data from the `timestamped_dataset.txt` file included with pyGSTi, and then build a `ProtocolData` object out of it so it can be used as input for `Protocol` objects.  We can pass `None` as the experiment design when constructing `data` because the stability analysis doesn't require any special structure to the circuits - it just requires the data to have timestamps."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialize the circuit structure details of the imported data.\n",
    "\n",
    "# This manually specifies the germ and fiducial structure for the imported data.\n",
    "fiducial_strs = ['{}','Gx','Gy','GxGx','GxGxGx','GyGyGy']\n",
    "germ_strs = ['Gi','Gx','Gy','GxGy','GxGyGi','GxGiGy','GxGiGi','GyGiGi','GxGxGiGy','GxGyGyGi','GxGxGyGxGyGy']\n",
    "log2maxL = 9 # log2 of the maximum germ power\n",
    "\n",
    "# Below we use the maxlength, germ and fuducial lists to create the GST structures needed for box plots.\n",
    "fiducials = [pygsti.objects.Circuit(fs) for fs in fiducial_strs]\n",
    "germs = [pygsti.objects.Circuit(g) for g in germ_strs]\n",
    "max_lengths = [2**i for i in range(0,log2maxL)]\n",
    "gss = pygsti.construction.make_lsgst_structs(('Gx','Gy','Gi'), fiducials, fiducials, germs, max_lengths)[-1] \n",
    "exp_design = pygsti.protocols.CircuitStructuresDesign([gss])\n",
    "\n",
    "ds = pygsti.io.load_dataset(\"../tutorial_files/timestamped_dataset.txt\")  # a DataSet\n",
    "data = pygsti.protocols.ProtocolData(exp_design, ds)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We then simply create a `StabilityAnalysis` protocol, and run it on the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "protocol = pygsti.protocols.StabilityAnalysis()\n",
    "results = protocol.run(data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note** the `StabilityAnalysis` protocol has a variety of optional arguments that can be used to optimize the tool to different circumstances. In this tutorial we won't discuss the full range of analyzes that can be performed using the `drift` module."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Inspecting the results\n",
    "\n",
    "Everything has been calculated, and we can now look at the results. If we print the returned results object (currently, essentially a container `StabilityAnalyzer` object), it will tell us whether instability was detected. If no instability is detected, then there is little else to do: the circuits are, as far as we can tell, stable, and most of the other results contained in a `StabilityAnalyzer` will not be very interesting. However, here instability is detected:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(results.stabilityanalyzer)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a workspace to show plots\n",
    "w = pygsti.report.Workspace()\n",
    "w.init_notebook_mode(connected=False, autodisplay=True) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Instability Detection Results : Power Spectra and the Frequencies of Instabilities\n",
    "\n",
    "The analysis is based on generating power spectra from the data. If the data is sampled from a time-invariant probability distribution, the powers in these spectra have an expected value of 1 and a known distribution (it is $ \\frac{1}{k}\\chi^2_k$ with $k$ depending on exactly what spectrum we're looking at). So we can test for violations of this, by looking for powers that are too high to be consistent with the time-invariant probability distribution hypothesis.\n",
    "\n",
    "A power spectrum is obtained for each circuit, and these can be averaged to obtain a single \"global\" power spectrum. This is plotted below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "w.PowerSpectraPlot(results)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Frequencies with power above the threshold in a spectrum are almost certainly components in the underlying (perhaps) time-varying probability for the corresponding circuit. We can access the frequencies above the threshold in the global power spectrum (which can't be assigned to a particular circuit, but *are* components in the probabilities for one or more circuits) as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(results.get_instability_frequencies())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To get the power spectrum, and detected significant frequencies, for a particular circuit we add an optional argument to the above functions: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "spectrumlabel = {'circuit':pygsti.obj.Circuit('Gx(Gi)^128')}\n",
    "print(\"significant frequencies: \", results.get_instability_frequencies(spectrumlabel))\n",
    "w.PowerSpectraPlot(results, spectrumlabel)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that all frequencies are in 1/units where \"units\" are the units of the time stamps in the `DataSet`. Some of the plotting functions display frequencies as Hertz, which is based on the assumption that these time stamps are in seconds. (in the future we will allow the user to specify the units of the time stamps).\n",
    "\n",
    "We can access a dictionary of all the circuits that we have detected as being unstable, with the values the detected frequencies of the instabilities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "unstablecircuits = results.get_unstable_circuits()\n",
    "# We only display the first 10 circuits and frequencies, as there are a lot of them!\n",
    "for ind, (circuit, freqs) in enumerate(unstablecircuits.items()):\n",
    "    if ind < 10: print(circuit.str, freqs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can jointly plot the power spectra for any set of circuits, by handing a dictionary of circuits to the `PowerSpectraPlot` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "circuits = {L: pygsti.obj.Circuit(None,stringrep='Gx(Gi)^'+str(L)+'Gx') for L in [1,2,4,16,64,128,256]}\n",
    "w.PowerSpectraPlot(results, {'circuit':circuits}, showlegend=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Instability Characterization Results : Probability Trajectories\n",
    "\n",
    "The tools also estimate the probability trajectories for each circuit, i.e., the probabilities to obtain each possible circuit outcome as a function of time. We can plot the estimated probability trajectory for any circuit of interest (or a selection of circuits, if we hand the plotting function a dictionary or list of circuits instead of a single circuit)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "circuit = pygsti.obj.Circuit(None, stringrep= 'Gx(Gi)^256GxGxGx')\n",
    "w.ProbTrajectoriesPlot(results.stabilityanalyzer, circuit, ('1',))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you simply want to access the time-varying distribution, use the `get_probability_trajectory()` method.\n",
    "\n",
    "The size of the instability in a circuit can be summarized by the amplitudes in front of the non-constant basis functions in our estimate of the probability trajectories. By summing these all up (and dividing by 2), we can get an upper bound on the maximum TVD between the instaneous probability distribution (over circuit outcomes) and the mean of this time-varying probability distribution, with this maximization over all times."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "results.get_max_tvd_bound(circuit)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you want to access this quantity for all unstable circuits, you can set `getmaxtvd = True` in the `get_unstable_circuits()` method. We can also access it's maximum over all the circuits as"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "results.get_maxmax_tvd_bound()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. Further plotting for data from structured circuits (e.g., GST circuits)\n",
    "\n",
    "If the data is from GST experiments, or anything with a GST-like structure of germs and fudicials (such as Ramsey circuits), we can create some extra plots, and a drift report.  In our case, our experiment design has `gss` as its circuit structure (see above).\n",
    "\n",
    "We can plot all of the power spectra, and probability trajectories, for amy (preparation fiducial, germ, measurement fiducial) triple. This shows how any instability changes, for this triple, as the germ power is increased."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "w.GermFiducialPowerSpectraPlot(results.stabilityanalyzer, gss, 'Gy', 'Gi', 'Gx', showlegend=True)\n",
    "w.GermFiducialProbTrajectoriesPlot(results.stabilityanalyzer, gss, 'Gy', 'Gi', 'Gx', ('0',), showlegend=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can make a boxplot that shows $\\lambda = -\\log_{10}(p)$ for each circuit, where $p$ is the p-value of the maximum power in the spectrum for that circuit. This statistic is a good summary for the evidence for instability in the circuit. Note that, for technical reasons, $\\lambda$ is truncated above at 16."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "w.ColorBoxPlot('driftdetector', gss, None, None, stabilityanalyzer=results.stabilityanalyzer)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The $\\lambda = -\\log_{10}(p)$ values do not *directly* tell us anything about the size of any detected instabilities. The $\\lambda$ statistics summarizes how certain we are that there is some instability, but if enough data is taken then even tiny instabilities will become obvious (and so we would have $\\lambda \\to \\infty$, except that the code truncates it to 16). \n",
    "\n",
    "The boxplot below summarizes the **size** of any detected instability with the bound on the maximal instantaneous TVD for each circuit, introduced above. Here this is zero for most circuits, as we did not detect any instability in those circuits - so our estimate for the probability trajectory is constant. The gradient of the colored boxes in this plot fairly closely mirror those in the plot above, but this is not always guaranteed to happen (e.g., if there is a lot more data for some of the circuits than others)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a boxplot of the maximum power in the power spectra for each sequence.\n",
    "w.ColorBoxPlot('driftsize', gss, None, None, stabilityanalyzer=results.stabilityanalyzer)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also create a report that contains all of these plots, as well as a few other things. But note that creating this report is currently fairly slow. Moreover, all the plots it contains have been demonstrated above, and everything else it contains can be accessed directly from the `StabilityAnalyzer` object. To explore all the things that are recorded in the `StabilityAnalyzer` object take a look at its `get` methods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "report = pygsti.report.construct_drift_report(results, title='Example Drift Report')\n",
    "report.write_html('../tutorial_files/DriftReport')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can now open the file [../tutorial_files/DriftReport/main.html](../tutorial_files/DriftReport/main.html) in your browser (Firefox works best) to view the report."
   ]
  }
 ],
 "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": 2
}