{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Mirror circuit benchmarks\n",
    "This notebook shows how to generate and analyze experiments that are similar to those performed in \"Measuring the Capabilities of Quantum Computers\" arXiv:2008.11294. For more details, see arXiv:2008.11294"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pygsti\n",
    "from pygsti.processors import QubitProcessorSpec, CliffordCompilationRules\n",
    "from pygsti.protocols import MirrorRBDesign as RMCDesign\n",
    "from pygsti.protocols import PeriodicMirrorCircuitDesign as PMCDesign\n",
    "from pygsti.protocols import ByDepthSummaryStatistics as SummaryStats\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_qubits = 4\n",
    "qubit_labels = ['Q'+str(i) for i in range(n_qubits)] \n",
    "gate_names = ['Gcnot'] + ['Gc{}'.format(i) for i in range(24)] \n",
    "availability = {'Gcnot':[('Q'+str(i),'Q'+str((i+1) % n_qubits)) for i in range(n_qubits)]}\n",
    "pspec = QubitProcessorSpec(n_qubits, gate_names, availability=availability, qubit_labels=qubit_labels)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clifford_compilations = {'absolute': pygsti.processors.CliffordCompilationRules.create_standard(pspec, verbosity=0)}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This cell is utility code for deciding how long the longest circuits should be. You\n",
    "# may want to choose this by hand instead.\n",
    "\n",
    "# A guess at the rough per-qubit error rate for picking the maximum depth\n",
    "# to go to at each width. Change as appropriate. Putting this too low will\n",
    "# just mean you run longer and more circuits than necessary.\n",
    "estimated_qubit_error_rate = 0.005\n",
    "\n",
    "# Heuristic for removing depths that are so long that you'll get no\n",
    "# useful data for w-qubit circuits. You could do this another way.\n",
    "def trim_depths(depths, w):\n",
    "    target_polarization = 0.01 \n",
    "    maxdepth = np.log(target_polarization)/(w * np.log(1 - estimated_qubit_error_rate))\n",
    "    trimmed_depths = [d for d in depths if d < maxdepth]\n",
    "    numdepths = len(trimmed_depths)\n",
    "    if numdepths < len(depths) and trimmed_depths[-1] < maxdepth:\n",
    "        trimmed_depths.append(depths[numdepths])\n",
    "    return trimmed_depths"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# This cell sets the circuit sampling parameters. These parameters are chosen\n",
    "# so as to replicate the experiments shown in Figs. 2 and 3 of arXiv:2008.11294.\n",
    "\n",
    "# The number of circuits per circuit shape (width and depth). Set this to 40\n",
    "# to replicate the experiments in arXiv:2008.11294. A smaller number might\n",
    "# be a good idea for a fast test run.\n",
    "circuits_per_shape = 10\n",
    "\n",
    "# The circuit widths to include in the benchmark.\n",
    "widths = [i for i in range(1, n_qubits + 1)]\n",
    "\n",
    "# The circuit depths to include, as a function of width.\n",
    "base_depths = [0,] +  [int(d) for d in 2**np.arange(1, 15)]\n",
    "depths = {w:trim_depths(base_depths,w) for w in widths}\n",
    "\n",
    "# The one-or-more qubit subsets to test for each width. You might want\n",
    "# to choose them more intelligently than here (in arXiv:2008.11294 we used the\n",
    "# qubits that were \"best\" according to their RB calibration data, although\n",
    "# there's nothing special about that strategy).\n",
    "qubit_lists = {w:[tuple([q for q in qubit_labels[:w]])] for w in widths}\n",
    "\n",
    "for w in widths:\n",
    "    print(w, depths[w], qubit_lists[w])\n",
    "    \n",
    "# Sets the two-qubit gate density in the circuits (where \"density\" refers\n",
    "# to the number of circuit locations occupied by a CNOT, with each CNOT\n",
    "# occupying two locations).\n",
    "xi = 1/8"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Samples randomized mirror circuits, using the `edgegrab` sampler with the two-qubit gate\n",
    "# density specified above. The `edgegrab` sampler is what is used in the experiments of\n",
    "# Figs. 2 and 3 in arXiv:2008.11294 (and not what was used for the experiments of Fig. 1d).\n",
    "edesigns = {}\n",
    "for w in widths:\n",
    "    for qs in qubit_lists[w]:\n",
    "        print(w, qs)\n",
    "        key = str(w) + '-' + '-'.join(qs) + '-' + 'RMCs'\n",
    "        edesigns[key] = RMCDesign(pspec, depths[w], circuits_per_shape, \n",
    "                                  clifford_compilations=clifford_compilations,\n",
    "                                  qubit_labels=qs, sampler='edgegrab', \n",
    "                                  samplerargs=[2 * xi,])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Samples periodic mirror circuits using the random germ selection algorithm specified\n",
    "# in arXiv:2008.11294, designed to match the RMCs sampled above (except that they're\n",
    "# periodic, not disordered).\n",
    "for w in widths:\n",
    "    for qs in qubit_lists[w]:\n",
    "        print(w, qs)\n",
    "        key = str(w) + '-' + '-'.join(qs) + '-' + 'PMCs'\n",
    "        edesigns[key] = PMCDesign(pspec, depths[w], circuits_per_shape,\n",
    "                                  clifford_compilations=clifford_compilations, \n",
    "                                  qubit_labels=qs, sampler='edgegrab', \n",
    "                                  samplerargs=[xi,])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "edesign = pygsti.protocols.CombinedExperimentDesign(edesigns)\n",
    "\n",
    "pygsti.io.write_empty_protocol_data('test_mirror_benchmark', edesign, clobber_ok=True)\n",
    "\n",
    "# All the circuits that need to be run\n",
    "circuits = edesign.all_circuits_needing_data\n",
    "\n",
    "# Shuffle the circuits: best to run them in a random order.\n",
    "np.random.shuffle(circuits)\n",
    "\n",
    "# Write this circuit list to file (the non-random order list is in the edesign folder).\n",
    "pygsti.io.write_circuit_list('test_mirror_benchmark/randomized_circuits.txt', circuits)\n",
    "\n",
    "# Convert to a list of OpenQASM format circuits. You may or may not want to use this \n",
    "# to get the pyGSTi circuits into a format that you can run on your device.\n",
    "qasm = [c.convert_to_openqasm(standard_gates_version='x-sx-rz') for c in circuits]\n",
    "\n",
    "# You'd then run the circuits of `qasm` or `circuits` and put them into a pyGSTi dataset\n",
    "# that replaces the empty dataset `test_mirror_benchmark/data/dataset.txt`. Below we instead\n",
    "# create simulated data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Creates simulated data.\n",
    "lindblad_error_coeffs = {g:{('H','Z'):0.01, ('S','X'):0.001} for g in gate_names if g != 'Gcnot'}\n",
    "lindblad_error_coeffs['Gcnot'] = {('H','XX'):0.05, ('S','XI'):0.005}\n",
    "noisemodel = pygsti.models.create_crosstalk_free_model(pspec, lindblad_error_coeffs=lindblad_error_coeffs)\n",
    "ds = pygsti.io.fill_in_empty_dataset_with_fake_data('test_mirror_benchmark/data/dataset.txt', noisemodel, num_samples=1000, seed=1234)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Loads the data\n",
    "data = pygsti.io.read_data_from_dir('test_mirror_benchmark')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The statistics to compute for each circuit.\n",
    "statistics = ['polarization', 'success_probabilities', 'success_counts', 'total_counts', 'two_q_gate_count']\n",
    "stats_generator = pygsti.protocols.SimpleRunner(SummaryStats(statistics_to_compute=statistics))\n",
    "\n",
    "# Computes the stats\n",
    "summary_data = stats_generator.run(data)\n",
    "\n",
    "# Turns this \"summary\" data into a DataFrame\n",
    "df = summary_data.to_dataframe('ValueName', drop_columns=['ProtocolName','ProtocolType'])\n",
    "\n",
    "# Adds a row that tells us which type of circuit the row is for. Will not work if the `keys` in the\n",
    "# edesign are changed to not include `RMCs` or `PMCs`. \n",
    "df['CircuitType'] = ['RMC' if 'RMCs' in p[0] else 'PMC' for p in df['Path']]\n",
    "\n",
    "# Redefines \"depth\" as twice what is in the Depth column, because the circuit generation code currently\n",
    "# uses a different convention to that used in arXiv:2008.11294.\n",
    "df['Depth'] = 2*df['Depth']\n",
    "\n",
    "# Puts the DataFrame into VBDataFrame object that can be used to create VB plots\n",
    "vbdf = pygsti.protocols.VBDataFrame(df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Creates a capability region plot like those in Fig. 3 of arXiv:2008.11294.\n",
    "fig, ax = pygsti.report.capability_region_plot(vbdf, figsize=(6, 8), scale=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extracts the data for a plot like Fig. 2a of arXiv:2008.11294.\n",
    "vb_min = {}\n",
    "for circuit_type in ('RMC', 'PMC'):   \n",
    "    vbdf1 = vbdf.select_column_value('CircuitType', circuit_type)\n",
    "    vb_min[circuit_type] = vbdf1.vb_data(metric='polarization', statistic='monotonic_min', no_data_action='min')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Creates the plot like those in Fig. 2a of arXiv:2008.11294. The inner squares\n",
    "# are the randomized mirror circuits, and the outer squares are the periodic\n",
    "# mirror circuits.\n",
    "spectral = pygsti.report.spectral\n",
    "fig, ax = pygsti.report.volumetric_plot(vb_min['PMC'], scale=1.9, cmap=spectral, figsize=(5.5,8))\n",
    "fig, ax = pygsti.report.volumetric_plot(vb_min['RMC'], scale=0.4, cmap=spectral, fig=fig, ax=ax, linescale=0.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Creates a plot like those in Fig. 1d of arXiv:2008.11294. But note\n",
    "# that these RMCs don't have the same sampling as those in Fig. 1d:\n",
    "# this is just the same type of plot from RMC data, not the same\n",
    "# type of RMCs. To get the same color map as in Fig. 1d, set cmap=None\n",
    "vbdf1 = vbdf.select_column_value('CircuitType', 'RMC')\n",
    "fig, ax = pygsti.report.volumetric_distribution_plot(vbdf1, figsize=(5.5,8), cmap=spectral)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Randomized mirror circuit data can also be used for \"RB\", i.e., estimating an average\n",
    "# gate error rate by fitting data to an exponential. Below shows how to do this, and\n",
    "# shows the RB error rates versus the number of qubits.\n",
    "rb = pygsti.protocols.RB(datatype='adjusted_success_probabilities', defaultfit='A-fixed')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rb_results = {}\n",
    "r = {}\n",
    "for key, subdata in data.items():\n",
    "    if 'RMCs' in key:\n",
    "        rb_results[key] = rb.run(subdata)\n",
    "        n_qubits = int(key.split('-')[0])\n",
    "        r[n_qubits] = rb_results[key].fits['A-fixed'].estimates['r']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "plt.plot(widths, [r[w] for w in widths], 'o')\n",
    "plt.xlabel('Number of Qubits')\n",
    "plt.ylabel(\"RB Error Rate\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}