{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The Extended Stabilizer Simulator"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Extended Simulator is a new method for classically simulating quantum circuits available in the latest release of [Qiskit-Aer](https://github.com/qiskit/qiskit-aer). \n",
    "\n",
    "This method is an implementation of the ideas published in the paper _Simulation of quantum circuits by low-rank stabilizer decompositions_ by Bravyi, Browne, Calpin, Campbell, Gosset & Howard, 2018, [**arXiv:1808.00128**](https://arxiv.org/abs/1808.00128).\n",
    "\n",
    "It uses a different representation of a quantum circuit, that gives it some unique capabilities. This notebook will give some examples of what the extended stabilizer method can do.\n",
    "\n",
    "For example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-08-19T17:17:32.305227Z",
     "start_time": "2019-08-19T17:17:32.300959Z"
    }
   },
   "outputs": [],
   "source": [
    "from qiskit import QuantumCircuit, transpile\n",
    "from qiskit_aer import AerSimulator\n",
    "from qiskit.visualization import plot_histogram\n",
    "\n",
    "import random"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-08-19T17:21:43.756467Z",
     "start_time": "2019-08-19T17:21:43.742436Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<qiskit.circuit.instructionset.InstructionSet at 0x7fbd18a23310>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "circ = QuantumCircuit(40, 40)\n",
    "\n",
    "# Initialize with a Hadamard layer\n",
    "circ.h(range(40))\n",
    "# Apply some random CNOT and T gates\n",
    "qubit_indices = [i for i in range(40)]\n",
    "for i in range(10):\n",
    "    control, target, t = random.sample(qubit_indices, 3)\n",
    "    circ.cx(control, target)\n",
    "    circ.t(t)\n",
    "circ.measure(range(40), range(40))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We've created a random circuit with just 60 gates, that acts on 40 qubits. But, because of the number of qubits, if we wanted to run this with say the statevector simulator then I hope you have access to terabytes of RAM!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-08-19T17:21:58.295935Z",
     "start_time": "2019-08-19T17:21:58.246239Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Simulation failed and returned the following error message:\n",
      "ERROR:  [Experiment 0] Insufficient memory to run circuit \"circuit-2\" using the statevector simulator.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "This succeeded?: False\n",
      "Why not? ERROR:  [Experiment 0] Insufficient memory to run circuit \"circuit-2\" using the statevector simulator.\n"
     ]
    }
   ],
   "source": [
    "# Create statevector method simulator\n",
    "statevector_simulator = AerSimulator(method='statevector')\n",
    "\n",
    "# Transpile circuit for backend\n",
    "tcirc = transpile(circ, statevector_simulator)\n",
    "\n",
    "# Try and run circuit\n",
    "statevector_result =  statevector_simulator.run(tcirc, shots=1).result()\n",
    "print('This succeeded?: {}'.format(statevector_result.success))\n",
    "print('Why not? {}'.format(statevector_result.status))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Extended Stabilizer method, in contrast, handles this circuit just fine. (Though it needs a couple of minutes!)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-08-19T17:22:26.158636Z",
     "start_time": "2019-08-19T17:22:22.534910Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "This succeeded?: True\n"
     ]
    }
   ],
   "source": [
    "# Create extended stabilizer method simulator\n",
    "extended_stabilizer_simulator = AerSimulator(method='extended_stabilizer')\n",
    "\n",
    "# Transpile circuit for backend\n",
    "tcirc = transpile(circ, extended_stabilizer_simulator)\n",
    "\n",
    "extended_stabilizer_result = extended_stabilizer_simulator.run(tcirc, shots=1).result()\n",
    "print('This succeeded?: {}'.format(extended_stabilizer_result.success))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## How does this work?\n",
    "If you're interested in how exactly we can handle such large circuits, then for a detailed explanation you can [read the paper!](https://arxiv.org/abs/1808.00128)\n",
    "\n",
    "For running circuits, however, it's important to just understand the basics.\n",
    "\n",
    "The Extended Stabilizer method is made up of two parts. The first is a method of decomposing quantum circuits into _stabilizer circuits_, a special class of circuit that can be efficiently simulated classically. The second is then a way of combining these circuits to perform measurements.\n",
    "\n",
    "The number of terms you need scales with the number of what we call _non-Clifford Gates_. At the moment, the method knows how to handle the following methods:\n",
    "```python\n",
    "circ.t(qr[qubit])\n",
    "circ.tdg(qr[qubit])\n",
    "circ.ccx(qr[control_1], qr[control_2], qr[target])\n",
    "circ.u1(rotation_angle, qr[qubit])\n",
    "```\n",
    "\n",
    "The simulator is also able to handle circuits of up to 63 qubits.\n",
    "     \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One thing that's important to note is these decompositions are approximate. This means that the results aren't exactly the same as with the State Vector simulator."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-08-19T17:23:52.873028Z",
     "start_time": "2019-08-19T17:23:52.863305Z"
    }
   },
   "outputs": [],
   "source": [
    "small_circ = QuantumCircuit(2, 2)\n",
    "small_circ.h(0)\n",
    "small_circ.cx(0, 1)\n",
    "small_circ.t(0)\n",
    "small_circ.measure([0, 1], [0, 1])\n",
    "# This circuit should give 00 or 11 with equal probability...\n",
    "expected_results ={'00': 50, '11': 50}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-08-19T17:23:54.879606Z",
     "start_time": "2019-08-19T17:23:53.249255Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "100 shots in 0.4958779811859131s\n"
     ]
    }
   ],
   "source": [
    "tsmall_circ = transpile(small_circ, extended_stabilizer_simulator)\n",
    "result = extended_stabilizer_simulator.run(\n",
    "    tsmall_circ, shots=100).result()\n",
    "counts = result.get_counts(0)\n",
    "print('100 shots in {}s'.format(result.time_taken))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-08-19T17:23:55.222602Z",
     "start_time": "2019-08-19T17:23:55.101958Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 18592 characters]\n",
      "text/plain": [
       "<Figure size 504x360 with 1 Axes>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plot_histogram([expected_results, counts],\n",
    "               legend=['Expected', 'Extended Stabilizer'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can control this approximation error using the `extended_stabilizer_approximation_error` in Qiskit Aer. The default error is 0.05. The smaller the error, the more precise the results, but also the longer your simulation will take and the more memory it will require."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-08-19T17:24:01.638958Z",
     "start_time": "2019-08-19T17:23:57.206364Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "100 shots in 1.404871940612793s\n"
     ]
    },
    {
     "data": {
      "image/png": "[encoded data, 17248 characters]\n",
      "text/plain": [
       "<Figure size 504x360 with 1 Axes>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Add runtime options for extended stabilizer simulator\n",
    "opts = {'extended_stabilizer_approximation_error': 0.03}\n",
    "\n",
    "reduced_error = extended_stabilizer_simulator.run(\n",
    "    tsmall_circ, shots=100, **opts).result()\n",
    "\n",
    "reduced_error_counts = reduced_error.get_counts(0)\n",
    "print('100 shots in {}s'.format(reduced_error.time_taken))\n",
    "plot_histogram([expected_results, reduced_error_counts],\n",
    "               legend=['Expected', 'Extended Stabilizer'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simulator Options"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are several other options you can tweak to control how the extended stabilizer method performs. What these options are and their explanation can all be found in the Qiskit Aer documentation. However, I want to highlight two important ones that can help to optimize your simulations.\n",
    "\n",
    "To perform measurements, the extended stabilizer method uses a Markov chain method to sample outcomes at random. This Markov chain has to be run for some time we call the 'mixing time' before it will start sampling, and has to be re-mixed for every circuit shot.\n",
    "\n",
    "If you expect your circuit output to be concentrated on just a few output states, then you can likely optimize your simulations by reducing the `extended_stabilizer_mixing_time` option.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-08-19T17:24:03.270544Z",
     "start_time": "2019-08-19T17:24:03.132881Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The circuit above, with 100 shots at precision 0.03 and default mixing time, needed 1s\n",
      "Dialing down the mixing time, we completed in just 1.4710919857025146s\n"
     ]
    }
   ],
   "source": [
    "print(\"The circuit above, with 100 shots at precision 0.03 \"\n",
    "      \"and default mixing time, needed {}s\".format(int(reduced_error.time_taken)))\n",
    "\n",
    "opts = {\n",
    "    'extended_stabilizer_approximation_error': 0.03,\n",
    "    'extended_stabilizer_mixing_time': 100\n",
    "}\n",
    "\n",
    "optimized = extended_stabilizer_simulator.run(\n",
    "    tsmall_circ, shots=100, **opts).result()\n",
    "\n",
    "print('Dialing down the mixing time, we completed in just {}s'.format(optimized.time_taken))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Similarly, if your circuit has some non-zero probability on all amplitudes (e.g. if it's a random circuit), then you can avoid this expensive re-mixing step to take multiple shots from the output at once. This can be enabled by setting `extended_stabilizer_measure_sampling=True`.\n",
    "\n",
    "For example, let's look again at the random circuit from the start of the tutorial, running for 100 shots:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-08-19T17:24:20.656799Z",
     "start_time": "2019-08-19T17:24:04.332417Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "100 shots took 29.634929895401 s\n"
     ]
    }
   ],
   "source": [
    "# We set these options here only to make the example run more quickly.\n",
    "opts = {'extended_stabilizer_mixing_time': 100}\n",
    "\n",
    "multishot = extended_stabilizer_simulator.run(\n",
    "    tcirc, shots=100, **opts).result()\n",
    "print(\"100 shots took {} s\".format(multishot.time_taken))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-08-19T17:24:21.070233Z",
     "start_time": "2019-08-19T17:24:20.918300Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "With the optimization, 100 shots took 0.4958779811859131 s\n"
     ]
    }
   ],
   "source": [
    "opts = {\n",
    "    'extended_stabilizer_measure_sampling': True,\n",
    "    'extended_stabilizer_mixing_time': 100\n",
    "}\n",
    "\n",
    "measure_sampling = extended_stabilizer_simulator.run(\n",
    "    circ, shots=100, **opts).result()\n",
    "print(\"With the optimization, 100 shots took {} s\".format(result.time_taken))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## When shall I use it?\n",
    "\n",
    "If you have smaller circuits with lots of non-Clifford gates, then the statevector method will likely perform better than the extended stabilizer. If however you want to look at circuits on many qubits, without needing access to high performance computation, then give this method a try!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-08-19T17:24:45.930560Z",
     "start_time": "2019-08-19T17:24:45.920568Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<h3>Version Information</h3><table><tr><th>Qiskit Software</th><th>Version</th></tr><tr><td>Qiskit</td><td>0.25.0</td></tr><tr><td>Terra</td><td>0.17.0</td></tr><tr><td>Aer</td><td>0.8.0</td></tr><tr><td>Ignis</td><td>0.6.0</td></tr><tr><td>Aqua</td><td>0.9.0</td></tr><tr><td>IBM Q Provider</td><td>0.12.2</td></tr><tr><th>System information</th></tr><tr><td>Python</td><td>3.7.7 (default, May  6 2020, 04:59:01) \n",
       "[Clang 4.0.1 (tags/RELEASE_401/final)]</td></tr><tr><td>OS</td><td>Darwin</td></tr><tr><td>CPUs</td><td>6</td></tr><tr><td>Memory (Gb)</td><td>32.0</td></tr><tr><td colspan='2'>Fri Apr 02 12:28:14 2021 EDT</td></tr></table>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<div style='width: 100%; background-color:#d5d9e0;padding-left: 10px; padding-bottom: 10px; padding-right: 10px; padding-top: 5px'><h3>This code is a part of Qiskit</h3><p>&copy; Copyright IBM 2017, 2021.</p><p>This code is licensed under the Apache License, Version 2.0. You may<br>obtain a copy of this license in the LICENSE.txt file in the root directory<br> of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.<p>Any modifications or derivative works of this code must retain this<br>copyright notice, and modified files need to carry a notice indicating<br>that they have been altered from the originals.</p></div>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import qiskit\n",
    "qiskit.__version__\n"
   ]
  },
  {
   "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.7.5"
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
