{
"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": [
"from __future__ import print_function\n",
"\n",
"# 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."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Loading ../tutorial_files/timeseries_data.txt: 100%\n"
]
}
],
"source": [
"ds = pygsti.io.load_tddataset(\"../tutorial_files/timeseries_data.txt\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Then we simply pass this data to `drift.do_stability_analysis()`. This has a variety of optional arguments that can be used to optimize the tool to different circumstances, but in this tutorial we won't discuss the full range of analyzes that can be performed using the `drift` module."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" - Formatting the data...done!\n",
" - Calculating power spectra...done!\n",
" - Running instability detection...done!\n",
" - Running instability characterization...done!\n"
]
}
],
"source": [
"# This'll take 5 - 10 minutes. This will be sped up in the future.\n",
"results = drift.do_stability_analysis(ds)"
]
},
{
"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 (a `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": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Instability *has* been detected, from tests at a global significance of 5.0%\n"
]
}
],
"source": [
"print(results)"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"\n",
"
Loading...
\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"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": [
{
"data": {
"text/html": [
"
\n",
"
\n",
"
\n",
"
\n",
"\n",
"
\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
""
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"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": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[0.001 0.003 0.005 0.008 0.011 0.016 0.2 ]\n"
]
}
],
"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": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"significant frequencies: [0.003]\n"
]
},
{
"data": {
"text/html": [
"
\n",
"
\n",
"
\n",
"
\n",
"\n",
"
\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
""
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"spectrumlabel = {'circuit':pygsti.obj.Circuit(None,stringrep='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": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Gx(Gx)^64GxGxGx [200]\n",
"Gx(Gi)^64 [3]\n",
"GxGxGx(Gi)^64 [3]\n",
"(Gx)^128 [200]\n",
"(Gx)^128GxGx [200]\n",
"Gx(Gx)^128GxGxGx [200]\n",
"GxGxGx(Gx)^128GxGxGx [200]\n",
"(Gi)^128Gx [3, 5]\n",
"(Gi)^128GxGxGx [3]\n",
"Gx(Gi)^128 [3]\n"
]
}
],
"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": [
{
"data": {
"text/html": [
"
\n",
"
\n",
"
\n",
"
\n",
"\n",
"
\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
""
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"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": [
{
"data": {
"text/html": [
"
\n",
"
\n",
"
\n",
"
\n",
"\n",
"
\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
""
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"circuit = pygsti.obj.Circuit(None, stringrep= 'Gx(Gi)^256GxGxGx')\n",
"w.ProbTrajectoriesPlot(results, 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": [
{
"data": {
"text/plain": [
"0.3453884720802307"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"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": [
{
"data": {
"text/plain": [
"0.49655622988939285"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"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."
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"# Initialize the circuit structure details of the imported data.\n",
"from pygsti.construction import std1Q_XYI\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(None,stringrep=fs) for fs in fiducial_strs]\n",
"germs = [pygsti.objects.Circuit(None,stringrep=mdl) for mdl in germ_strs]\n",
"max_lengths = [2**i for i in range(0,log2maxL)]\n",
"gssList = pygsti.construction.make_lsgst_structs(std1Q_XYI.gates, fiducials, fiducials, germs, max_lengths) \n",
"gss = gssList[-1]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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": 15,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
""
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"w.GermFiducialPowerSpectraPlot(results, gss, 'Gy', 'Gi', 'Gx', showlegend=True)\n",
"w.GermFiducialProbTrajectoriesPlot(results, 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": 16,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"
\n",
"
\n",
"
\n",
"\n",
"
\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
""
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"w.ColorBoxPlot('driftdetector', gss, None, None, stabilityanalyzer=results)"
]
},
{
"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": 17,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"
\n",
"
\n",
"
\n",
"
\n",
"\n",
"
\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
""
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Create a boxplot of the maximum power in the power spectra for each sequence.\n",
"w.ColorBoxPlot('driftsize', gss, None, None, stabilityanalyzer=results)"
]
},
{
"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": 18,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"*** Creating workspace ***\n",
"*** Generating switchboard ***\n",
"*** Generating tables ***\n",
" driftSummaryTable took 0.303495 seconds\n",
" driftDetailsTable took 0.000219 seconds\n",
"*** Generating plots ***\n",
" GlobalPowerSpectraPlot took 0.098495 seconds\n",
" GermFiducialPowerSpectraPlot took 194.147931 seconds\n",
" GermFiducialProbTrajectoriesPlot took 561.129772 seconds\n",
" driftdetectorColorBoxPlot took 108.564209 seconds\n",
" driftsizeColorBoxPlot took 0.933101 seconds\n",
"*** Merging into template file ***\n",
" Rendering drift_switchBd took 9.4e-05 seconds\n",
" Rendering topSwitchboard took 2.5e-05 seconds\n",
" Rendering driftSummaryTable took 0.001453 seconds\n",
" Rendering driftDetailsTable took 0.001144 seconds\n",
" Rendering GlobalPowerSpectraPlot took 0.064204 seconds\n",
" Rendering GermFiducialPowerSpectraPlot took 54.85668 seconds\n",
" Rendering GermFiducialProbTrajectoriesPlot took 208.715839 seconds\n",
" Rendering driftdetectorColorBoxPlot took 0.252762 seconds\n",
" Rendering driftsizeColorBoxPlot took 0.237954 seconds\n",
"Output written to ../tutorial_files/DriftReport directory\n",
"*** Report Generation Complete! Total time 1129.54s ***\n"
]
},
{
"data": {
"text/plain": [
""
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"drift.driftreport.create_drift_report(results, gss, \"../tutorial_files/DriftReport\", \n",
" title=\"Example Drift Report\", verbosity=10)"
]
},
{
"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
}