UniB_Aer.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# UniB Workshop 26/27 - Qiskit Aer\n",
    "## IBM Q\n",
    "### François Varchon (francois.varchon@fr.ibm.com)\n",
    "##### Thanks to Christopher J. Wood\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## How do we simulate NISQ devices using classical computers?\n",
    "\n",
    "### Outline\n",
    "\n",
    "1. Ideal Quantum Circuits\n",
    "1. Ideal simulations with Qiskit Aer\n",
    "1. Noisy quantum circuits\n",
    "1. Noisy simulations with Qiskit Aer\n",
    "1. Noise model examples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Housekeeping: uncomment this to suppress deprecation warnings\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "# Import the basics\n",
    "import numpy as np\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import qiskit\n",
    "from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister\n",
    "from qiskit.tools import visualization as qplt\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Ideal Quantum Circuits\n",
    "\n",
    "The ideal circuit based model of quantum computation consists of three pieces:\n",
    "1. **States:** A *complex statevector* $|\\psi\\rangle$ satisfying $|\\langle\\psi|\\psi\\rangle|^2=1$\n",
    "1. **Gates:** A *unitary* matrix $U$ satisfying $U^\\dagger U = U U^\\dagger = \\mathbb{1}$.\n",
    "1. **Measurements:** A *projective-valued measure (PVM)* $\\{\\Pi_j\\}$ s.t. $p_j = |\\langle \\psi |\\Pi_j|\\psi\\rangle|^2$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ideal Quantum States\n",
    "\n",
    "Quantum states my be expressed as a complex vector with respect to the computation basis:\n",
    "\n",
    "Eg for a single qubit: $|\\psi\\rangle = \\alpha |0\\rangle + \\beta |1\\rangle = [\\alpha, \\beta]$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "psi = np.array([1, 1j]) / np.sqrt(2)\n",
    "\n",
    "qplt.plot_bloch_multivector(psi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ideal Quantum Gates\n",
    "\n",
    "Ideal quantum gates are operations that map statevectors to statevectors:\n",
    "\n",
    "$ |\\psi_{i+1}\\rangle = U_i |\\psi_i\\rangle $\n",
    "\n",
    "These operations may be represented by *unitary matrices*.\n",
    "\n",
    "* A matrix $U$ is unitary if and only if $U^\\dagger U = U^\\dagger U = \\mathbb{1}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ideal Measurement\n",
    "\n",
    "* In simulation we can direct access to the wavefunction, however in experiment we can only perform measurements to extract information.\n",
    "\n",
    "* An ideal measurement is a non-unitary operation which causes the wavefunction to *collapse*\n",
    "\n",
    "* We will consider measurement in the computational ($Z$) basis:\n",
    "\n",
    "If we measure a state $|\\psi\\rangle = \\alpha_0 |0\\rangle + \\alpha_1 |1\\rangle$\n",
    "\n",
    "For a single measurement:\n",
    "* We will observe outcome \"0\" with probability $p_0 = |\\alpha_0|^2$, \n",
    "* After observing \"0\" the system will have collapsed to the state $|0\\rangle$.\n",
    "* Similarly outcome \"1\" will be observed with probability $p_1 = |\\alpha_1|^2$, leaving the system in state $|1\\rangle$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Ideal Simulations with Qiskit Aer\n",
    "\n",
    "* Qiskit Aer is a local simulation provider for Qiskit.\n",
    "* We will briefly show how to do ideal simulations with Qiskit-Aer.\n",
    "* These are taken from the Aer Provider tutorial in `qiskit-tutorials`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Qiskit Aer simulator backends\n",
    "\n",
    " Qiskit Aer currently includes three simulator backends\n",
    " \n",
    "* `QasmSimulator`: Simulates experiments (returns counts)\n",
    "* `StatevectorSimulator`: Simulates an ideal shot of an experiment (returns statevector)\n",
    "* `UnitarySimulator`: Simulates an ideal circuit (return circuit unitary)\n",
    "\n",
    "These backends are found in the `Aer` provider with the names `qasm_simulstor`, `statevector_simulator` and `unitary_simulator` respectively"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiskit import Aer\n",
    "import qiskit.providers.aer as aer\n",
    "from qiskit.providers.aer import QasmSimulator, StatevectorSimulator, UnitarySimulator"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## QasmSimulator\n",
    "\n",
    "* Designed to mimic an actual device. \n",
    "* Executes `QuantumCircuit` for $N$ shots\n",
    "* Returns count dictionary containing $N$ final vlaues of `ClassicalRegisters` in the circuit.\n",
    "* The circuit may contain *gates*, *measure*, *reset*, *conditionals*\n",
    "\n",
    "### Example: prepare and measure a Bell-state\n",
    "\n",
    "* Circuit prepares a 2-qubit Bell-state $|\\psi\\rangle = \\frac{1}{2}(|0,0\\rangle + |1,1 \\rangle)$\n",
    "* Measures both qubits"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Construct quantum circuit\n",
    "qr = QuantumRegister(2, 'qr')\n",
    "cr = ClassicalRegister(2, 'cr')\n",
    "circ = QuantumCircuit(qr, cr)\n",
    "circ.h(qr[0])\n",
    "circ.cx(qr[0], qr[1])\n",
    "circ.measure(qr, cr)\n",
    "\n",
    "# Plot the circuit\n",
    "print(circ.draw())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Select the QasmSimulator from the Aer provider\n",
    "simulator = Aer.get_backend('qasm_simulator')\n",
    "\n",
    "# Execute and get counts\n",
    "result = qiskit.execute(circ, simulator).result()\n",
    "counts = result.get_counts(circ)\n",
    "qplt.plot_histogram(counts, title='Bell-State counts')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Returning measurements outcomes for each shot\n",
    "\n",
    "* `QasmSimulator` also supports returning a list of outcomes for each shot. \n",
    "* Enabled by setting the keyword argument `memory=True` in `compile` or `execute`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Execute and get memory\n",
    "result = qiskit.execute(circ, simulator, shots=10, memory=True).result()\n",
    "memory = result.get_memory(circ)\n",
    "print(memory)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Starting simulation with a custom initial state\n",
    "\n",
    "* The default initial state for `QasmSimulator` is $|\\psi_0\\rangle = |0^{\\otimes n}\\rangle$.\n",
    "* `QasmSimulator` allows setting a custom initial statevector for the simulation. \n",
    "* Done using `backend_options={\"initial_statevector\": psi_0}` in `QasmSimulator.run` or `qiskit.execute`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set the initial state\n",
    "opts = {\"initial_statevector\": np.array([1, 0, 0, 1] / np.sqrt(2))}\n",
    "\n",
    "# Select the QasmSimulator from the Aer provider\n",
    "simulator = Aer.get_backend('qasm_simulator')\n",
    "\n",
    "# Execute and get counts\n",
    "result = qiskit.execute(circ, simulator, backend_options=opts).result()\n",
    "counts = result.get_counts(circ)\n",
    "qplt.plot_histogram(counts, title=\"Bell initial statevector\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Requirements for initial statevector\n",
    "\n",
    "* $|\\psi_0\\rangle$ must be valid quantum state: $|\\langle\\psi_0|\\psi_0\\rangle|=1$.\n",
    "* If not an error will be raised. \n",
    "* The same $\\psi_0\\rangle$ will be used for *all* circuits in Qobj.\n",
    "* If it is wrong size for number of qubits will raise an error\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Statevector Simulator\n",
    "\n",
    "* Simulates a single-shot of a `QuantumCircuit` and returns the final quantum statevector.\n",
    "* The circuit may contain *gates*, and also *measure*, *reset*, and *conditional* operations.\n",
    "* If measure/reset operations are in the circuit:\n",
    "  * The final state will be a *conditional state* for the outcome of those operations.\n",
    "  * The measurement outcome (counts) will also be returned (with $N=1$ shots).\n",
    "  \n",
    "* A custom initial statevector $|\\psi_0\\rangle$ can be used (same as for `QasmSimulator`)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example: Bell state preparation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Construct quantum circuit without measure\n",
    "qr = QuantumRegister(2, 'qr')\n",
    "circ = QuantumCircuit(qr)\n",
    "circ.h(qr[0])\n",
    "circ.cx(qr[0], qr[1])\n",
    "\n",
    "# Plot the circuit\n",
    "print(circ.draw())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Select the StatevectorSimulator from the Aer provider\n",
    "simulator = Aer.get_backend('statevector_simulator')\n",
    "\n",
    "# Execute and get counts\n",
    "result = qiskit.execute(circ, simulator).result()\n",
    "statevector = result.get_statevector(circ)\n",
    "qplt.plot_state_city(statevector, title='Bell state')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example: Bell state preparation *with measurement*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Construct quantum circuit with measure\n",
    "qr = QuantumRegister(2, 'qr')\n",
    "cr = ClassicalRegister(2, 'cr')\n",
    "circ = QuantumCircuit(qr, cr)\n",
    "circ.h(qr[0])\n",
    "circ.cx(qr[0], qr[1])\n",
    "circ.measure(qr, cr)\n",
    "\n",
    "# Select the StatevectorSimulator from the Aer provider\n",
    "simulator = Aer.get_backend('statevector_simulator')\n",
    "\n",
    "# Execute and get counts\n",
    "result = qiskit.execute(circ, simulator).result()\n",
    "statevector = result.get_statevector(circ)\n",
    "qplt.plot_state_city(statevector, title='Bell state post-measurement')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Unitary Simulator\n",
    "\n",
    "* Constructs the full $n$-qubit unitary matrix for a `QuantumCircuit`.\n",
    "* Done by applying each gate matrix to an initial identity matrix. \n",
    "* Circuit may only contain *gates*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example: Bell-state preparation circuit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Construct an empty quantum circuit\n",
    "qr = QuantumRegister(2)\n",
    "circ = QuantumCircuit(qr)\n",
    "circ.h(qr[0])\n",
    "circ.cx(qr[0], qr[1])\n",
    "\n",
    "# Select the UnitarySimulator from the Aer provider\n",
    "simulator = Aer.get_backend('unitary_simulator')\n",
    "\n",
    "# Execute and get counts\n",
    "result = qiskit.execute(circ, simulator).result()\n",
    "unitary = result.get_unitary(circ)\n",
    "print(\"Final Unitary:\")\n",
    "print(unitary)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Custom initial unitary\n",
    "\n",
    "* Instead of the identity matrix, a custom initial unitary can be used\n",
    "* `backend_options={\"initial_unitary\": U0}`\n",
    "\n",
    "### Requirements for initial unitary\n",
    "\n",
    "* $U_0$ must be unitary ($U_0^\\dagger U_0 = \\mathbb{1}$).\n",
    "* The same $U_0$ will be used for *all* circuits in Qobj.\n",
    "* If these are violated an error will be raised"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Noisy Quantum Circuits\n",
    "\n",
    "A noisy quantum circuit is described using formalism of *Open Quantum Systems*. For the circuit model noisy quantum computation consists of three pieces:\n",
    "\n",
    "1. **States:** A density matrix $\\rho$, satisfying $\\rho \\ge 0$, $Tr[\\rho]=1$.\n",
    "1. **Gates:** A *completely-positive trace-preserving* (CPTP) map $\\cal{E}$.\n",
    "1. **Measurement:** A *positive operator valued measure* (POVM) $\\{E_j\\}$ s.t. $p_j = Tr[E_j\\rho]$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Noisy Quantum States\n",
    "\n",
    "* The density matrix $\\rho$ can be thought of our *average knowledge* about a state averaged over many experiments.\n",
    "\n",
    "### Example: Noisy initialization\n",
    "\n",
    "* We have a device that is supposed to generate a single-qubit state $|0\\rangle$.\n",
    "* 10% of the time it generates the state $|1\\rangle$ instead.\n",
    "* We describe the *average output* of this noisy process as a state:\n",
    "\n",
    "$$\\rho = \\frac{9}{10} * |0\\rangle\\langle0| + \\frac{1}{10} * |1\\rangle\\langle1|$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rho = 0.9 * np.array([[1, 0], [0, 0]]) + 0.1 * np.array([[0, 0,], [0, 1]])\n",
    "qplt.plot_bloch_multivector(rho)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qplt.plot_state_city(rho)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Noisy Quantum Gates\n",
    "\n",
    "Noisy quantum gates (CPTP-maps) a *roughly speaking* operations which map density matrices to density matrices. Thus, like density matrices themselves they represent our *average knowledge* about evolution of a system (under certain assumptions).\n",
    "\n",
    "They may always be written in the **Kraus representation** as\n",
    "\n",
    "$$\\mathcal{E}(\\rho) = \\sum_j K_j \\rho K_j^\\dagger, \\quad\\mbox{where}\\quad \\sum_j K_j^\\dagger K_j = \\mathbb{1}$$\n",
    "\n",
    "The set of matrices $\\{K_j\\}$ *non-uniquely* defines a noisy quantum gate."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A *unique* representation of a CPTP map is the linear *superoperator* $S$, which is the matrix that acts on *vectorized* density matrices:\n",
    "\n",
    "$$|\\mathcal{E}\\rangle\\!\\rangle = S|\\rho\\rangle\\!\\rangle$$\n",
    "\n",
    "In terms of the *Kraus* representation the superoprator is given by $S = \\sum_j K_j^*\\otimes K_j$ and is a $4^n * 4^n$ matrix."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example: Noisy X gate\n",
    "\n",
    "* Suppose we have a device that is supposed to apply an $X$ gate to a qubit.\n",
    "* 10% of the time it does nothing (an identity gate)\n",
    "* We describe the *average output* of this noisy process on an *arbitrary* state as the CPTP map:\n",
    "\n",
    "$$\\mathcal{E}(\\rho) = \\frac{9}{10}\\, X\\rho X + \\frac{1}{10}\\, \\rho$$\n",
    "\n",
    "* In terms of a Kraus representation this is $\\{\\sqrt{9/10}\\, X, \\sqrt{1/10}\\, I\\}$\n",
    "* **NOTE:** This is a special case of a CPTP map called a *mixed unitary channel*, where all the Kraus operators are a unitary matrix applied with some probability"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### How do we apply CPTP maps in simulation?\n",
    "\n",
    "* If we want to do a density matrix simulation we require $16 * 4^n$ complex floating point numbers to store the state itself\n",
    "* We require an additional $16 * 16^m$ complex floats for each $m$-qubit superoperator to apply.\n",
    "* This greatly limits the size of systems we can simulate\n",
    "\n",
    "### Alternative approach: quantum trajectories\n",
    "\n",
    "* We can mimic CPTP evolution by *unravelling* noisy gates into (possibly non-unitary) operations that can be applied to a pure quantum state\n",
    "* We must rescale the statevector after such a gate to preserve normalization\n",
    "* We recover the average dynamics of the map by *averaging* over the output of many shots\n",
    "* **NOTE:** This trades of memory requirement for simulation time (less RAM more shots)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Quantum Trajectory state update\n",
    "\n",
    "* Let the system be in the pure state $|psi\\rangle$\n",
    "* Consider a Kraus gate error $\\{K_j\\}$\n",
    "* Compute probability of applying Kraus operator $j$ to the state from:\n",
    "$$ p_j = |\\langle\\psi|K_j^\\dagger K_j\\psi\\rangle|^2 $$\n",
    "* Sample a single Kraus operator to apply from the discrete distribution give by $\\{p_j\\}$\n",
    "* Update the state to:\n",
    "$$|\\psi^\\prime\\rangle = \\frac{1}{\\sqrt{p_j}} K_j|\\psi\\rangle$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Important: Mixed unitary errors\n",
    "\n",
    "* This can be done more efficiently for mixed unitary error channels\n",
    "* *The probabilities are known ahead of time and are independent of the current state*\n",
    "* In this case for each noisy operation we can sample the noisy unitary to apply in advance and insert it into the circuit. Then apply the noisy circuits to the initial state and average over the output"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Noisy Measurement\n",
    "\n",
    "* Noisy measurement is represented by a POVM $\\{E_j\\}$, where $E_j \\ge 0$, $\\sum_j E_j = \\mathbb{1}$.\n",
    "\n",
    "* Measurement outcome probability is given by $p_j = Tr[E_j\\rho]$.\n",
    "\n",
    "* For simulation purposes we consider *two types* of noisy measurement.\n",
    " \n",
    "  1. **Quantum measurement error:** The error can be expressed as a CPTP map *before* an ideal measurement.\n",
    "  2. **Classical readout error:** The error is expressed as a stochastic map that applies to the *classical register* recording the measurement outcome *after measurment*."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Classical Readout Error\n",
    "\n",
    "* Readout errors are specified by a list of assignment probabilities vectors $P(A|B)$.\n",
    "  * $A$ is the *recorded* classical bit value\n",
    "  * $B$ is the *true* bit value returned from the measurement \n",
    " \n",
    "Eg for 1 qubits: $ P(A|B) = [P(A|0), P(A|1)]$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Noisy Simulations with Qiskit Aer\n",
    "\n",
    "* We will demonstrate how noisy simulations may be done in Qiskit-Aer\n",
    "* The `QasmSimulator` supports a trajector based noise simulation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Qiskit Aer Noise Module\n",
    "\n",
    "* The `aer.noise` module contains Python classes to build customized noise models for simulation.\n",
    "* Three key classes:\n",
    "\n",
    "  1. The `NoiseModel` class which stores a noise model used for noisy simulation \n",
    "  2. The `QuantumError` class which describes CPTP gate errors. These can be applied\n",
    "    * After *gate* or *reset* instructions\n",
    "    * Before *measure* instructions.\n",
    "  3. The `ReadoutError` class which describes classical readout errors.\n",
    "* Functions to generate common `QuantumError` objects in `aer.noise.errors`.\n",
    "* Functions to automatically generate a noise model in `aer.noise.device`\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Importing necessary functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import the noise model object\n",
    "from qiskit.providers.aer.noise import NoiseModel\n",
    "\n",
    "# Import all standard errors\n",
    "from qiskit.providers.aer.noise.errors import *  \n",
    "\n",
    "# Import noise base classes\n",
    "from qiskit.providers.aer.noise.errors.quantum_error import QuantumError\n",
    "from qiskit.providers.aer.noise.errors.readout_error import ReadoutError"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Standard Quantum Errors\n",
    "\n",
    "Rather than deal with the `QuantumError` object directly, many helper functions exist to automatically generate a specific type of parameterized quantum error. These are contained in the `noise.errors` submodule and include functions for many common errors types used in quantum computing research. The function names and the type of error they return are:\n",
    "\n",
    "| Standard error function | Details |\n",
    "| --- | --- |\n",
    "| `kraus_error` | a general n-qubit CPTP error channel given as a list of Kraus matrices $[K_0, ...]$. |\n",
    "| `mixed_unitary_error` | an n-qubit mixed unitary error given as a list of unitary matrices and probabilities $[(U_0, p_0),...]$. |\n",
    "| `coherent_unitary_error` | an n-qubit coherent unitary error given as a single unitary matrix $U$. |\n",
    "| `pauli_error` | an n-qubit Pauli error channel (mixed unitary) given as a list of Pauli's and probabilities $[(P_0, p_0),...]$ |\n",
    "| `depolarizing_error` | an n-qubit depolarizing error channel parameterized by a depolarizion probability $p$. |\n",
    "| `reset_error` | a single-qubit reset error parameterized by a probabilities $p_0, p_1$ of reseting to the $|0\\rangle$, $|1\\rangle$ state.|\n",
    "| `thermal_relaxation_error` | a single qubit thermal relaxation channel parameterized by relaxation time constants $T_1$, $T_2$, gate time $t$, and excited state thermal population $p_1$. |\n",
    "| `phase_amplitude_damping_error` | A single-qubit generalized combined phase and amplitude damping error channel given by an amplitude damping parameter $\\lambda$, a phase damping parameter $\\gamma$, and an excited state thermal population $p_1$. |\n",
    "| `amplitude_damping_error` |  A single-qubit generalized amplitude damping error channel given by an amplitude damping parameter $\\lambda$, and an excited state thermal population $p_1$. |\n",
    "| `phase_damping_error` |  A single-qubit phase damping error channel given by a phase damping parameter $\\gamma$ |\n",
    "\n",
    "### Combining Errors\n",
    "\n",
    "`QuantumError` instances can be combined to produce new Quantum Errors using:\n",
    " * Composition: $\\cal{E}=\\cal{E_1}\\circ\\cal{E_2}$ as `error1.compose(error2)`\n",
    " * Tensor product: $\\cal{E}=\\cal{E_1}\\otimes\\cal{E_2}$ as `error1.kron(error2)`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example: Depolarizing errors\n",
    "\n",
    "For example to construct a 5% single-qubit depolarizing error:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Construct a 1-qubit depolarizing error\n",
    "depol1 = depolarizing_error(0.05, 1)\n",
    "print(depol1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Composition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "depol1sq = depol1.compose(depol1)\n",
    "print(depol1sq)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Tensor Product"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "depol2 = depol1.kron(depol1)\n",
    "print(depol2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Adding errors to a Noise Model\n",
    "\n",
    "When adding a quantum error to a noise model we must specify the type of *instruction* that it acts on, and what qubits to apply it to. There are three cases for Quantum Errors:\n",
    " \n",
    " 1. All-qubit quantum error\n",
    " 2. Specific qubit quantum error\n",
    " 3. Non-local quantum error\n",
    " \n",
    "### All-qubit quantum error\n",
    "\n",
    "This applies the same error to any occurence of an instruciton, regardless of which qubits it acts on.\n",
    "\n",
    "It is added as `noise_model.add_all_qubit_quantum_error(error, instructions)`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create an empty noise model\n",
    "noise_model = NoiseModel()\n",
    "\n",
    "# Add depolarizing error to all single qubit u1, u2, u3 gates\n",
    "noise_model.add_all_qubit_quantum_error(depol1, ['u1', 'u2', 'u3'])\n",
    "\n",
    "# Print noise model info\n",
    "print(noise_model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Specific qubit quantum error\n",
    "\n",
    "This applies the error to any occurence of an instruction acting on a specified list of qubits. Note that the order of the qubit matters: For a 2-qubit gate an error applied to qubits [0, 1] is different to one applied to qubits [1, 0] for example.\n",
    "\n",
    "It is added as `noise_model.add_quantum_error(error, instructions, qubits)`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create an empty noise model\n",
    "noise_model = NoiseModel()\n",
    "\n",
    "# Add depolarizing error to all single qubit u1, u2, u3 gates on qubit 0 only\n",
    "noise_model.add_quantum_error(depol1, ['u1', 'u2', 'u3'], [0])\n",
    "\n",
    "# Print noise model info\n",
    "print(noise_model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Non-local qubit quantum error\n",
    "\n",
    "This applies an error to a specific set of noise qubits after any occurence of an instruction acting on a specific of gate qubits.\n",
    "\n",
    "It is added as `noise_model.add_quantum_error(error, instructions, instr_qubits, error_qubits)`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create an empty noise model\n",
    "noise_model = NoiseModel()\n",
    "\n",
    "# Add depolarizing error on qubit 2 forall single qubit u1, u2, u3 gates on qubit 0 \n",
    "noise_model.add_nonlocal_quantum_error(depol1, ['u1', 'u2', 'u3'], [0], [2])\n",
    "\n",
    "# Print noise model info\n",
    "print(noise_model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Executing a noisy simulation with a noise model\n",
    "\n",
    "* To execute a noisy simulation we pass the noise model object to `QasmSimulator.run` or `execute` using the `noise_model` kwarg.\n",
    "* Eg: `qiskit.execute(circuits, QasmSimulator(), noise_model=noise)`\n",
    "\n",
    "**Important:** *When running a noisy simulation make sure you compile your qobj to the same basis gates as the noise model!*\n",
    "\n",
    "This can be done using `NoiseModel.basis_gates`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "# Noise Model Examples\n",
    "\n",
    "We will now give some examples of noise models\n",
    "\n",
    "### Test circuit\n",
    "\n",
    "* For our demonstrations we wil use a simple test circuit generating a n-qubit GHZ state:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "# Simulator\n",
    "simulator = QasmSimulator()\n",
    "\n",
    "# System Specification\n",
    "n_qubits = 4\n",
    "qr = QuantumRegister(n_qubits)\n",
    "cr = ClassicalRegister(n_qubits)\n",
    "circ = QuantumCircuit(qr, cr)\n",
    "\n",
    "# Test Circuit\n",
    "circ.h(qr[0])\n",
    "for qubit in range(n_qubits - 1):\n",
    "    circ.cx(qr[qubit], qr[qubit + 1])\n",
    "circ.barrier()\n",
    "circ.measure(qr, cr)\n",
    "print(circ.draw())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Test execution\n",
    "result = qiskit.execute(circ, simulator).result()\n",
    "qplt.plot_histogram(result.get_counts(0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Basic Error Model Examples"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "## Example 1: Basic bit-flip error noise model\n",
    "\n",
    "Lets consider a simple toy noise model example common in quantum information theory research:\n",
    "\n",
    "* When applying a single qubit gate, flip the state of the qubit with probability `p_gate1.\n",
    "* When applying a 2-qubit gate apply single-qubit errors to each qubit.\n",
    "* When reseting a qubit reset to 1 instead of 0 with probability `p_reset`\n",
    "* When measuring a qubit, flip the state of the qubit before with probability `p_meas`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Example error probabilities\n",
    "p_reset = 0.03\n",
    "p_meas = 0.1\n",
    "p_gate1 = 0.05\n",
    "\n",
    "# QuantumError objects\n",
    "error_reset = pauli_error([('X', p_reset), ('I', 1 - p_reset)])\n",
    "error_meas = pauli_error([('X',p_meas), ('I', 1 - p_meas)])\n",
    "error_gate1 = pauli_error([('X',p_gate1), ('I', 1 - p_gate1)])\n",
    "error_gate2 = error_gate1.kron(error_gate1)\n",
    "\n",
    "# Add errors to noise model\n",
    "noise_bit_flip = NoiseModel()\n",
    "noise_bit_flip.add_all_qubit_quantum_error(error_reset, \"reset\")\n",
    "noise_bit_flip.add_all_qubit_quantum_error(error_meas, \"measure\")\n",
    "noise_bit_flip.add_all_qubit_quantum_error(error_gate1, [\"u1\", \"u2\", \"u3\"])\n",
    "noise_bit_flip.add_all_qubit_quantum_error(error_gate2, [\"cx\"])\n",
    "\n",
    "print(noise_bit_flip)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Executing the noisy simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# Compile previous GHZ circuit using noise model basis gates\n",
    "qobj = qiskit.compile(circ, simulator,\n",
    "               basis_gates=noise_bit_flip.basis_gates)\n",
    "\n",
    "# Run the noisy simulation\n",
    "result_bit_flip = simulator.run(qobj, noise_model=noise_bit_flip).result()\n",
    "counts_bit_flip = result_bit_flip.get_counts(0)\n",
    "\n",
    "# Plot noisy output\n",
    "qplt.plot_histogram(counts_bit_flip)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 2: T1/T2 thermal relaxation\n",
    "\n",
    "* Now consider a more realistic error model based on thermal relaxation with the qubit environment\n",
    "* Each qubit parameterized by a thermal relaxation time constant $T_1$ and a dephasing time constant $T_2$.\n",
    "* Note that we must have $T_2 \\le 2 T_1$\n",
    "* Error rates on instructions are determined by gate time and qubit $T_1$, $T_2$ values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# T1 and T2 values for qubits 0-3\n",
    "T1s = np.random.normal(50e3, 10e3, 4) # Sampled from normal distribution mean 50 microsec\n",
    "T2s = np.random.normal(70e3, 10e3, 4)  # Sampled from normal distribution mean 50 microsec\n",
    "\n",
    "# Truncate random T2s <= T1s\n",
    "T2s = np.array([min(T2s[j], 2 * T1s[j]) for j in range(4)])\n",
    "\n",
    "# Instruction times (in nanoseconds)\n",
    "time_u1 = 0   # virtual gate\n",
    "time_u2 = 50  # (single X90 pulse)\n",
    "time_u3 = 100 # (two X90 pulses)\n",
    "time_cx = 300\n",
    "time_reset = 1000  # 1 microsecond\n",
    "time_measure = 1000 # 1 microsecond\n",
    "\n",
    "# QuantumError objects\n",
    "errors_reset = [thermal_relaxation_error(t1, t2, time_reset)\n",
    "                for t1, t2 in zip(T1s, T2s)]\n",
    "errors_measure = [thermal_relaxation_error(t1, t2, time_measure)\n",
    "                  for t1, t2 in zip(T1s, T2s)]\n",
    "errors_u1  = [thermal_relaxation_error(t1, t2, time_u1)\n",
    "              for t1, t2 in zip(T1s, T2s)]\n",
    "errors_u2  = [thermal_relaxation_error(t1, t2, time_u2)\n",
    "              for t1, t2 in zip(T1s, T2s)]\n",
    "errors_u3  = [thermal_relaxation_error(t1, t2, time_u3)\n",
    "              for t1, t2 in zip(T1s, T2s)]\n",
    "errors_cx = [[thermal_relaxation_error(t1a, t2a, time_cx).kron(\n",
    "             thermal_relaxation_error(t1b, t2b, time_cx))\n",
    "              for t1a, t2a in zip(T1s, T2s)]\n",
    "               for t1b, t2b in zip(T1s, T2s)]\n",
    "\n",
    "# Add errors to noise model\n",
    "noise_thermal = NoiseModel()\n",
    "for j in range(4):\n",
    "    noise_thermal.add_quantum_error(errors_reset[j], \"reset\", [j])\n",
    "    noise_thermal.add_quantum_error(errors_measure[j], \"measure\", [j])\n",
    "    noise_thermal.add_quantum_error(errors_u1[j], \"u1\", [j])\n",
    "    noise_thermal.add_quantum_error(errors_u2[j], \"u2\", [j])\n",
    "    noise_thermal.add_quantum_error(errors_u3[j], \"u3\", [j])\n",
    "    for k in range(4):\n",
    "        noise_thermal.add_quantum_error(errors_cx[j][k], \"cx\", [j, k])\n",
    "\n",
    "print(noise_thermal)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Executing the noisy simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compile previous GHZ circuit using noise model basis gates\n",
    "qobj = qiskit.compile(circ, simulator,\n",
    "                      basis_gates=noise_thermal.basis_gates)\n",
    "\n",
    "# Run the noisy simulation\n",
    "result_thermal = simulator.run(qobj, noise_model=noise_thermal).result()\n",
    "counts_thermal = result_thermal.get_counts(0)\n",
    "\n",
    "# Plot noisy output\n",
    "qplt.plot_histogram(counts_thermal)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Intermediate Error Model Examples\n",
    "\n",
    "We will now show how one can simulate more complicated error processes using the Qiskit-Aer noise model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 3: Correlated spectator qubit noise\n",
    "\n",
    "* For our second example we will consider correlated spectator qubit noise on the CX gates in our test circuit.\n",
    "\n",
    "* Let us consider the case where our 4 qubits are arrange on the verticies of a square.\n",
    "* Apply a CNOT gate with control $c$ and target $t$.\n",
    "* With probability $p$ there is a of a correlated $ZZ$ error between $t$ and spectator qubit $s$ coupled to the control:\n",
    "\n",
    "```\n",
    "c -- t                     I -- I              I -- Z\n",
    "|    |  --> (1-p) * CNOT * |    | + p * CNOT * |    |\n",
    "s -- x                     I -- I              Z -- I\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "# Define an example coupling map for the square\n",
    "coupling_map = [[0, 1], [0, 3], [2, 1], [2, 3]]\n",
    "\n",
    "# Example error probabilities\n",
    "p_01 = 0.03\n",
    "p_03 = 0.01\n",
    "p_21 = 0.02\n",
    "p_23 = 0.05\n",
    "\n",
    "# QuantumError objects\n",
    "def zz_error(p):\n",
    "    return pauli_error([('ZZ', p), ('II', 1 - p)])\n",
    "\n",
    "# Add errors to noise model\n",
    "noise_corr_cx = NoiseModel()\n",
    "noise_corr_cx.add_nonlocal_quantum_error(zz_error(p_01), \"cx\", [0, 1], [1, 3])\n",
    "noise_corr_cx.add_nonlocal_quantum_error(zz_error(p_01), \"cx\", [0, 3], [1, 3])\n",
    "noise_corr_cx.add_nonlocal_quantum_error(zz_error(p_01), \"cx\", [2, 1], [1, 3])\n",
    "noise_corr_cx.add_nonlocal_quantum_error(zz_error(p_01), \"cx\", [2, 3], [1, 3])\n",
    "\n",
    "print(noise_corr_cx)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Executing the noisy simulation\n",
    "\n",
    "Now when we execute we have to remember to compile using the coupling map so that the CX gates of the circuit will be mapped to those in the coupling map that are specified in the noise model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compile previous GHZ circuit using noise model basis gates\n",
    "qobj = qiskit.compile(circ, simulator,\n",
    "                      coupling_map=coupling_map,\n",
    "                      basis_gates=noise_corr_cx.basis_gates)\n",
    "\n",
    "# Run the noisy simulation\n",
    "result_corr_cx = simulator.run(qobj, noise_model=noise_corr_cx).result()\n",
    "counts_corr_cx = result_corr_cx.get_counts(0)\n",
    "\n",
    "# Plot noisy output\n",
    "qplt.plot_histogram(counts_corr_cx)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 4: Coherent leakage error\n",
    "\n",
    "* *Leakage* is an error that occurs when a qubit is actually a subspace in a higher-dimensional system\n",
    "* Coherent leakage errors are unitary errors between the qubit subspace and the other energy levels.\n",
    "* Now show how one can model *leakage error* even though the simulation is only a qubit simulation. \n",
    "\n",
    "### To model coherent leakage:\n",
    "\n",
    "* We add a second quantum register of ancilla qubits to represent the leakage levels for each original qubit\n",
    "* We add  second set of classical bits to flag when leakage has occured based on the state of the ancilla qubits\n",
    "* For a single qubit we have \n",
    "$$|\\psi\\rangle \\mapsto  |\\phi\\rangle\\otimes|\\psi\\rangle$$ \n",
    "\n",
    "\n",
    "* where our qubit is now a four level system\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "|0\\rangle \\equiv |0\\rangle\\otimes|0\\rangle \\\\\n",
    "|1\\rangle \\equiv |0\\rangle\\otimes|1\\rangle \\\\\n",
    "|2\\rangle \\equiv |1\\rangle\\otimes|0\\rangle \\\\\n",
    "|3\\rangle \\equiv |1\\rangle\\otimes|1\\rangle\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "* If we measure the ancilla qubit state $|\\phi\\rangle$ and get an outcome of 1 we know that the qubit state is in the leakage subspace."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define additional registers\n",
    "qa = QuantumRegister(n_qubits)\n",
    "ca = ClassicalRegister(n_qubits)\n",
    "circ_leakage = QuantumCircuit(qr, qa, cr, ca)\n",
    "\n",
    "# Add original circuit\n",
    "circ_leakage += circ\n",
    "\n",
    "# Measure ancillas to flag leakage\n",
    "circ_leakage.barrier(qr, qa)\n",
    "circ_leakage.measure(qa, ca)\n",
    "print(circ_leakage)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Representing coherent leakage error\n",
    "\n",
    "* Represent a coherent leakage error as a 2-qubit unitary matrix.\n",
    "* Add it to the noise model as a non-local error for gates on the gate qubit\n",
    "* Non-local error  effects both the gate qubit and corresponding ancilla qubit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.linalg import expm\n",
    "\n",
    "# Ideal X-90 drive Hamiltonian\n",
    "H_ideal =  0.25 * np.pi * np.kron(np.eye(2), np.array([[0, 1], [1, 0]]))\n",
    "U_ideal = U_noisy = expm(-1j * H_ideal)\n",
    "\n",
    "# Coherent unitary error due to leakage transitions in Hamiltonian\n",
    "H_error = 0.1 * np.array([[0, 0, 0, 0], [0, 0, 1, 0], [0, 1, 0, 0], [0, 0, 0, 0]])\n",
    "U_error = coherent_unitary_error(np.dot(U_ideal.conj().T, expm(-1j * (H_ideal + H_error))))\n",
    "\n",
    "# Build noise model and add to H and CX\n",
    "# For CX we will just do a tensor-product of the single qubit errors\n",
    "U_error2 = U_error.kron(U_error)\n",
    "\n",
    "# Add errors to noise model, but be careful about qubit numbers!\n",
    "noise_leakage = NoiseModel()\n",
    "for i in range(n_qubits):\n",
    "    noise_leakage.add_nonlocal_quantum_error(U_error, ['u1', 'u2', 'u3'],\n",
    "                                             [i], [i, n_qubits + i])\n",
    "    for j in range(n_qubits):\n",
    "        if j != i:\n",
    "            noise_leakage.add_nonlocal_quantum_error(U_error2, 'cx',\n",
    "                                                     [i, j],\n",
    "                                                     [i, n_qubits + i, j, n_qubits + j])\n",
    "print(noise_leakage)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Executing the noisy simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compile previous GHZ circuit using noise model basis gates\n",
    "qobj = qiskit.compile(circ_leakage, simulator,\n",
    "                      basis_gates=noise_leakage.basis_gates)\n",
    "\n",
    "result_leakage = simulator.run(qobj, noise_model=noise_leakage).result()\n",
    "counts_leakage = result_leakage.get_counts(0)\n",
    "\n",
    "# Plot noisy output\n",
    "qplt.plot_histogram(counts_leakage)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Advanced Noise Model Examples\n",
    "\n",
    "These advanced noise modules use features of *Qiskit Aer* that aren't fully integrated with Qiskit Terra. They involve modifying a compiled qobj to insert special instructions that will be added to the Terra QuantumCircuit API in a future release."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 5: Adding correlated classical readout error\n",
    "\n",
    "* Now let us consider adding a correlated readout error on a multi-qubit measurement. \n",
    "\n",
    "* For 2-qubits our assignment probabilities are\n",
    "$$P(A|B) = [P(A|00), P(A|01), P(A|10), P(A|11)]$$\n",
    "\n",
    "* We will build the noise model for 2-qubit readout errors on pairs of qubits `[0, 1]` and `[2, 3]` \n",
    "* *Note: We could also do a 4-qubit correlated readout error, but we'll stick to 2-qubit to keep the matrices small*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Assignment probabilities for measure qubits [0, 1]\n",
    "p_00_qubits01 = [0.8, 0.1, 0.1, 0.0]\n",
    "p_01_qubits01 = [0.2, 0.6, 0.1, 0.1]\n",
    "p_10_qubits01 = [0.2, 0.1, 0.6, 0.1]\n",
    "p_11_qubits01 = [0.1, 0.2, 0.2, 0.5]\n",
    "p_qubits01 = [p_00_qubits01, p_01_qubits01, p_10_qubits01, p_11_qubits01]\n",
    "\n",
    "# Assignment probabilities for measure qubits [2, 3]\n",
    "p_00_qubits23 = [0.7, 0.1, 0.1, 0.1]\n",
    "p_01_qubits23 = [0.2, 0.6, 0.1, 0.1]\n",
    "p_10_qubits23 = [0.2, 0.1, 0.6, 0.1]\n",
    "p_11_qubits23 = [0.2, 0.2, 0.2, 0.4]\n",
    "p_qubits23 = [p_00_qubits23, p_01_qubits23, p_10_qubits23, p_11_qubits23]\n",
    "\n",
    "# Build readout error objects\n",
    "ro_error01 = ReadoutError(p_qubits01)\n",
    "print(ro_error01)\n",
    "\n",
    "# Build readout error objects\n",
    "ro_error23 = ReadoutError(p_qubits23)\n",
    "print(ro_error23)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we add these errors to the noise model, note that we dont actually have to construct the `ReadoutError` object, we could just add the readout errors using the lists of probabilities `p_qubitsij` directly "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "noise_corr_ro = NoiseModel()\n",
    "noise_corr_ro.add_readout_error(ro_error01, [0, 1])\n",
    "noise_corr_ro.add_readout_error(ro_error01, [2, 3])\n",
    "print(noise_corr_ro)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Executing the noisy simulation\n",
    "\n",
    "* To actually use this model in a simulation we need to modify our qobj\n",
    "* We need to make the `measure` instruction a multi-qubit measurement\n",
    "* Editing is necessary as the Terra compiler decomposes all `measure` circuit instructions into successive single-qubit `measure` qobj instructions.\n",
    "* We can do this using helper functions in `aer.utils.qobj_utils`\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compile previous GHZ circuit using noise model basis gates\n",
    "qobj = qiskit.compile(circ, simulator,\n",
    "                      basis_gates=noise_corr_ro.basis_gates)\n",
    "\n",
    "# Find measure instruction positions in qobj and drop them\n",
    "meas_pos = aer.utils.qobj_utils.get_instr_pos(qobj, 0, 'measure')\n",
    "\n",
    "# Delete measure from qobj\n",
    "instructions = qobj.experiments[0].instructions\n",
    "[instructions.pop(j) for j in reversed(meas_pos)]\n",
    "\n",
    "# Add new multi-qubit measure instructions\n",
    "meas01 = aer.utils.qobj_utils.measure_instr([0, 1], [0, 1])\n",
    "meas23 = aer.utils.qobj_utils.measure_instr([2, 3], [2, 3])\n",
    "instructions += [meas01, meas23]\n",
    "\n",
    "# Replace experiment instructions with new modified ones\n",
    "qobj.experiments[0].instructions = instructions\n",
    "\n",
    "# Run the noisy simulation\n",
    "result_corr_ro = simulator.run(qobj, noise_model=noise_corr_ro).result()\n",
    "counts_corr_ro = result_corr_ro.get_counts(0)\n",
    "\n",
    "# Plot noisy output\n",
    "qplt.plot_histogram(counts_corr_ro)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 6: Always-on cross-talk noise\n",
    "\n",
    "Finally we will add always on cross-talk to do this we will define a noise model on a *custom unitary gate* (which in our case will be an identity), that we can then assign an error to in the noise model. We will insert this gate in between all slices to simulate an always on multi-qubit cross-talk error.\n",
    "\n",
    "Let us consider an error that consists of correlated two-body $XX$ errors between all the qubits on a square"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p01 = 0.01\n",
    "p03 = 0.03\n",
    "p21 = 0.02\n",
    "p23 = 0.04\n",
    "ct_error = pauli_error([('IIXX', p01),\n",
    "                        ('XIIX', p03),\n",
    "                        ('IXXI', p21),\n",
    "                        ('XXII', p23),\n",
    "                        ('IIII', 1 - p01 - p03 - p21 - p23)],\n",
    "                        standard_gates=True)\n",
    "print(ct_error)\n",
    "\n",
    "# Add error to a noise model for custom instruction \"ct\"\n",
    "noise_ct = NoiseModel()\n",
    "noise_ct.add_nonlocal_quantum_error(ct_error, \"ct\", [0], [0, 1, 2, 3])\n",
    "# the following aren't really necessary if we always apply ct to qubit-0\n",
    "noise_ct.add_nonlocal_quantum_error(ct_error, \"ct\", [1], [0, 1, 2, 3]) \n",
    "noise_ct.add_nonlocal_quantum_error(ct_error, \"ct\", [2], [0, 1, 2, 3])\n",
    "noise_ct.add_nonlocal_quantum_error(ct_error, \"ct\", [3], [0, 1, 2, 3])\n",
    "\n",
    "print(noise_ct)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Executing noisy simulation\n",
    "\n",
    "* We need to modify the Qobj for custom gate errors\n",
    "* First we compile to a qobj\n",
    "* Next we insert our custom identity matrix with name `\"ct\"` after the Hadamard and CX gates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compile previous GHZ circuit using noise model basis gates\n",
    "qobj = qiskit.compile(circ, simulator,\n",
    "                      basis_gates=noise_ct.basis_gates)\n",
    "\n",
    "# Create custom gate\n",
    "ct_gate = aer.utils.qobj_utils.unitary_instr(np.eye(2), [0], label=\"ct\")\n",
    "\n",
    "# Because of our noise basis gates the circuit will be compiled to U and CX gates:\n",
    "U_pos = aer.utils.qobj_utils.get_instr_pos(qobj, 0, 'U')\n",
    "CX_pos = aer.utils.qobj_utils.get_instr_pos(qobj, 0, 'CX')\n",
    "for pos in reversed(sorted(U_pos + CX_pos)):\n",
    "    aer.utils.qobj_utils.insert_instr(qobj, 0, ct_gate, pos + 1)\n",
    "\n",
    "# Run the noisy simulation\n",
    "result_ct = simulator.run(qobj, noise_model=noise_ct).result()\n",
    "counts_ct = result_ct.get_counts(0)\n",
    "\n",
    "# Plot noisy output\n",
    "qplt.plot_histogram(counts_ct)"
   ]
  },
  {
   "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}