EntanglementPurification.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Entanglement Purification: The BBPSSW Protocol\n",
    "***\n",
    "\n",
    "\n",
    "## Summary and Outline\n",
    "***\n",
    "In this tutorial we demonstrate the use of IBM Qiskit and the IBM Q Experience by implementing a protocol for entanglement purification, the BBPSSW protocol. This technique from quantum information processing, which is especially required for quantum communication, allows to generate highly-entangled qubit pairs from noisy sources. In the protocol presented here, we will use two entangled qubit pairs to increase the fidelity of one of the pairs with respect to a maximally-entangled Bell state.\n",
    "\n",
    "In **Part I** we will first give a very brief introduction to entanglement purification itself. We present the necessary details to this topic in a depth, so that we understand the following implementation in Qiskit with this knowledge.\n",
    "\n",
    "In **Part II** we give a detailed description of how to implement the entanglement purification protocol in Qiskit. We focus on explaining the necessary steps in detail and providing a comprehensive description of the use of Qiskit. The presented contents include the creation of quantum circuits from single- and two-qubit gates, the implementation of quantum measurements and the visualization of quantum circuits. Furthermore, we present how to use the tomography functions integrated in Qiskit to perform quantum state tomography, which is necessary to demonstrate the functioning of entanglement purification.\n",
    "\n",
    "After we have generated the quantum circuit for entanglement purification, we will perform a simulation of this circuit with the built-in quantum simulator of Qiskit in **Part III**. We show how the density matrix of the qubit pair for which we want to increase the entanglement can be reconstructed step by step from the simulation results. In order to verify the correctness of our Qiskit implementation of the purification protocol, we compare the simulation results with exact calculations.\n",
    "\n",
    "The simulation in Part III is based on a perfect quantum computer. In **Part IV** we extend this simulation and account for imperfections in the setup. Through the IBM Q Experience, we will read the properties of a real IBM Q device, such as gate and measurement accuracy, and relaxation and decoherence rates. From these properties we will show how to create a simple noise model with the help of integrated Qiskit functions, which can then be used in the simulation. From the simulation results, we will again reconstruct the density matrix of our entangled qubit pair and we can gain a first impression of how well purification could work on a real device.\n",
    "\n",
    "Finally, in **Part V**, we will run the purification protocol on a real IBM Q device. We show how to access such a device through the IBM Q Experience and execute the quantum circuit for entanglement purification that we have previously constructed. As in the previous parts, we will reconstruct the density matrix of the real qubit pair and check whether the entanglement can be increased.\n",
    "\n",
    "In the last part, **Part VI**, we will summarize our results and propose further steps that could logically extend this tutorial.\n",
    "\n",
    "In between, we will provide a few short exercises that can be solved with the gained knowledge. (So far the solutions of the exercises are given directly in the notebook due to the participation in the IBM Q Awards. However, they could also be stored in a separate file, e.g. if the tutorial is officially made available.)\n",
    "\n",
    "## Learning Objectives\n",
    "***\n",
    "\n",
    "We give a comprehensive explanation of the following basic features of Qiskit and the IBM Q Experience:\n",
    "\n",
    "* Composition of quantum circuits from single- and two-qubit gates\n",
    "* Generation of parametrized quantum circuits\n",
    "* Implementation of qubit measurements\n",
    "* Execution of quantum circuits using the simulator integrated in Qiskit\n",
    "* Execution of quantum circuits on real IBM Q devices using the IBM Q Experience\n",
    "* Extraction of measurement results\n",
    "\n",
    "In addition, we explain more advanced uses of Qiskit and the topics are:\n",
    "\n",
    "* Performance of quantum state tomography\n",
    "* Postprocession of measurement results for conditional quantum state tomography\n",
    "* Extraction of noise parameters from real devices using the IBM Q Experience\n",
    "* Creation of simple noise models\n",
    "* Simulation of noisy systems using Qiskit's integrated simulator\n",
    "* Calculation of an ensemble average over many quantum systems\n",
    "\n",
    "We also give an introduction to entanglement purification and explain the associated quantum protocol, which we will implement step by step in Qiskit. However, our main focus is on the use of Qiskit and we refer readers interested in in-depth explanations of entanglement purification to the references we provide.\n",
    "\n",
    "\n",
    "## Prior Knowledge\n",
    "***\n",
    "\n",
    "To understand the tutorial, the following topics from the field of quantum information theory are required:\n",
    "\n",
    "* Common single- and two-qubit gates\n",
    "* Qubit entanglement or the notion of maximal entanglement\n",
    "* Bell states\n",
    "* Mixed quantum states and density matrices\n",
    "* Quantum state fidelity\n",
    "\n",
    "The following topics are helpful, but not absolutely necessary, to work through this tutorial:\n",
    "\n",
    "* Quantum state tomography\n",
    "\n",
    "Since this tutorial is provided as Jupyter notebook, basic knowledge of the following Python libraries is recommended:\n",
    "\n",
    "* NumPy\n",
    "\n",
    "\n",
    "## Contents\n",
    "***\n",
    "* **Part I: Introduction to Entanglement Purification and the BBPSSW Protocol**\n",
    "  * 1.1 Concept\n",
    "  * 1.2 Quantum Protocol\n",
    "  \n",
    "  \n",
    "* **Part II: Implementation of Quantum Circuit for BBPSSW Protocol in Qiskit**\n",
    "  * 2.1 Basic Code Imports\n",
    "  * 2.2 Definitions\n",
    "  * 2.3 Quantum Register and Quantum Circuit\n",
    "    * 2.3.1 Quantum Register\n",
    "    * 2.3.2 Circuit for Initialization of Source State\n",
    "    * 2.3.3 Implementation of Twirl Operation by Random Bilateral Rotation\n",
    "    * 2.3.4 Exchange of $|\\Psi^-\\rangle$ and $|\\Phi^+\\rangle$ Components\n",
    "    * 2.3.5 Application of Bilateral CNOT Gate\n",
    "    * 2.3.6 Concatenation of Quantum Circuits\n",
    "  * 2.4 (Conditional) Quantum State Tomography\n",
    "  \n",
    "  \n",
    "* **Part III: Simulation of the BBPSSW Protocol using Qiskit**\n",
    "   * 3.1 Single Quantum System\n",
    "   * 3.2 Ensemble Average\n",
    "     * 3.2.1 Definition of Auxiliary Functions\n",
    "     * 3.2.2 Simulation\n",
    "\n",
    "\n",
    "* **Part IV: Including Noise in the Simulation Model**\n",
    "  * 4.1 IBM Q Account\n",
    "  * 4.2 Creation of Noise Model in Qiskit\n",
    "  * 4.3 Simulation Including Noise\n",
    "    * 4.3.1 Initial State and Initial Fidelity\n",
    "    * 4.3.2 Entanglement Purification in the Presence of Noise\n",
    "\n",
    "\n",
    "* **Part V: Entanglement Purification on a Real IBM Q Device**\n",
    "  * 5.1 IBM Q Experience\n",
    "    * 5.1.1 Device Queue, Job Status, and Job Monitor\n",
    "    * 5.1.2 Personal Credits\n",
    "  * 5.2 Initial State and Initial Fidelity\n",
    "  * 5.3 Execution of the BBPSSW Protocol on a Real Device\n",
    "\n",
    "\n",
    "* **Part VI: Conclusion and Further Ideas**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Part I: Introduction to Entanglement Purification and the BBPSSW Protocol\n",
    "***\n",
    "\n",
    "## 1.1 Concept\n",
    "\n",
    "The concept of **entanglement purification** was introduced by C. H. Bennett *et al.* [1]. Following the names of the authors, the original protocol is also called the **BBPSSW protocol**. The basic idea of entanglement purification is to create an (almost) maximally entangled Bell state by distilling entanglement from a larger number of qubit pairs that are not maximally entangled. For this reason, this concept is also referred to as **entanglement distillation**. Purification usually has to be performed by two remote parties, so only local operations, measurements and classical two-way communication are allowed for entanglement purification protocols. \n",
    "\n",
    "Many applications of quantum information processing, such as quantum teleportation [2] or entanglement-based quantum cryptography [3], use the maximally entangled *Bell states* $\\{|\\Phi^{\\pm}\\rangle, |\\Psi^{\\pm}\\rangle\\}$,\n",
    "\n",
    "$$ |\\Phi^{\\pm}\\rangle = \\frac{1}{\\sqrt{2}}(|00\\rangle \\pm |11\\rangle), \\\\\n",
    "|\\Psi^{\\pm}\\rangle = \\frac{1}{\\sqrt{2}}(|01\\rangle \\pm |10\\rangle). $$\n",
    "\n",
    "The feasibility or security of these protocols depends on the \"quality\" of the available Bell states. However, real sources of entangled qubit pairs or real lossy quantum channels are very likely to produce imperfect Bell states due to preparation errors or imperfect qubit transmission. The aim of entanglement purification is to compensate for this loss in entanglement and restore a maximally entangled Bell state. That is why entanglement purification is also a fundamental component of a quantum repeater [4].\n",
    "\n",
    "At the beginning we assume a source $S$ of entangled qubit pairs, which generates two-qubit states $\\rho_\\textrm{S}$ with some **state fidelity** $F$ towards the maximally-entangled Bell state $|\\Psi^-\\rangle$,\n",
    "\n",
    "$$\n",
    "F = \\langle \\Psi^- | \\rho_\\textrm{S} | \\Psi^- \\rangle.\n",
    "$$\n",
    "\n",
    "As already mentioned, imperfection can be caused by the production process itself or by the transmission of the qubits via lossy quantum channels.\n",
    "\n",
    "The two parties that each have one of the qubits of the source state $\\rho_\\textrm{S}$ are usually called Alice and Bob. The whole scenario is shown in the following figure:\n",
    "\n",
    "<img vspace=\"30\">\n",
    "<img src=\"images/fig01.png\" width=\"400\" vspace=\"30\">\n",
    "\n",
    "The basic idea of the BBPSSW protocol is that Alice and Bob take two copies of the state $\\rho_\\textrm{S}$ and perform certain quantum operations locally on their qubits to finally generate a single two-qubit state with larger fidelity $F' > F$. When this is achieved, the process can be iterated to produce even larger fidelities until a maximally entangled Bell state $| \\Psi^- \\rangle$ is obtained. \n",
    "\n",
    "In the following section we will explain the BBPSSW protocol step by step and show how Alice and Bob can increase the entanglement of the state $\\rho_\\textrm{S}$.\n",
    "\n",
    "## 1.2 Quantum Protocol\n",
    "\n",
    "### Step 1: Creation of a Werner State through Random Bilateral Rotation\n",
    "\n",
    "From the two-qubit source state $\\rho_\\textrm{S}$ the two parties Alice and Bob must first create a so-called **Werner state** $W_F$,\n",
    "\n",
    "$$\n",
    "W_F = F |\\Psi^-\\rangle \\langle \\Psi^-| + \\frac{1-F}{3}\\left(|\\Psi^+\\rangle \\langle \\Psi^+|+|\\Phi^+\\rangle \\langle \\Phi^+|+|\\Phi^-\\rangle \\langle \\Phi^-|\\right).  \n",
    "$$\n",
    "\n",
    "Such a state is diagonal in the Bell basis $\\{ |\\Phi^{+}\\rangle, |\\Phi^{-}\\rangle, |\\Psi^{+}\\rangle, |\\Psi^{-}\\rangle \\}$, has a fidelity of $F$ towards the Bell state $|\\Psi^-\\rangle$, and overlaps equally with the other three Bell states by $(1-F)/3$. \n",
    "\n",
    "Alice and Bob can achieve this by applying a so-called **twirl** operation [5] that keeps the $|\\Psi^-\\rangle$ component of the state $\\rho_\\textrm{S}$, balances the components of the other three Bell states, and removes all off-diagonal elements. A twirl can be implemented by bilateral operations $B_k (k \\in \\{x, y, z\\})$ that rotate each qubit of the pair by an angle of $\\pi/2$ around the $k$ axis. A random selection from the following set $\\{ U_i \\}$ of 12 combinations of bilateral rotations $B_k$,\n",
    "\n",
    "$$\n",
    "\\{U_i\\}\n",
    "=\n",
    "\\{\n",
    "I, \\\\ B_x B_x, \\\\ B_y B_y, \\\\ B_z B_z, \\\\\n",
    "B_x B_y, \\\\ B_y B_z, \\\\ B_z B_x, \\\\ B_y B_x, \n",
    "\\\\ B_x B_y B_x B_y, \\\\ B_y B_z B_y B_z, \\\\ B_z B_x B_z B_x, \\\\ B_y B_x B_y B_x\n",
    "\\},  \n",
    "$$  \n",
    "\n",
    "transforms the state $\\rho_\\textrm{S}$ into the arithmetic average\n",
    "\n",
    "$$\n",
    "W_F = \\frac{1}{12} \\sum_{i=1}^{12} U_i \\rho_\\textrm{S} U_i^\\dagger,\n",
    "$$\n",
    "\n",
    "which will be exactly the desired Werner state $W_F$. Here $I$ denotes the identity operation. An intuitive explanation of the twirl, which comes from molecular physics, can be found in Ref. [2] and is left to the interested reader.\n",
    "\n",
    "Since Alice and Bob randomly select an operation from the set $\\{ U_i \\}$ in a realistic quantum communication scenario, we will also randomly select such an operation in the later implementation in Qiskit. In order to obtain the arithmetic average, we will have to perform several runs of the purification protocol.\n",
    "\n",
    "We summarize the quantum operations performed locally by Alice and Bob in a graphical illustration for the sake of clarity:\n",
    "\n",
    "<img vspace=\"30\">\n",
    "<img src=\"images/fig02.png\" width=\"1000\" vspace=\"30\">\n",
    "\n",
    "Here Alice's and Bob's qubit are referred to as $q_A$ and $q_B$, and $U_{i, A(B)}$ denotes Alice's (Bob's) single-qubit part of the randomly selected operation $U_i$.\n",
    "\n",
    "### Step 2: Exchange of $|\\Psi^-\\rangle$ and $|\\Phi^+\\rangle$ Components\n",
    "\n",
    "After Alice and Bob have prepared the state $W_F$, one of the parties rotates its qubit by an angle of $\\pi$ around the $y$ axis, which corresponds to using a $Y$ gate,\n",
    "\n",
    "$$\n",
    "Y = \n",
    "\\begin{pmatrix}\n",
    "0 & -i \\\\\n",
    "i & 0 \\\\\n",
    "\\end{pmatrix}.\n",
    "$$\n",
    "\n",
    "The result of such a unilateral rotation is that the components of $|\\Psi^-\\rangle$ and $|\\Phi^+\\rangle$ are exchanged in the Werner state $W_F$. Therefore, the operation creates a state $\\rho_F$,\n",
    "\n",
    "$$\n",
    "\\rho_F = F |\\Phi^+\\rangle \\langle \\Phi^+| + \\frac{1-F}{3}\\left(|\\Psi^+\\rangle \\langle \\Psi^+|+|\\Psi^-\\rangle \\langle \\Psi^-|+|\\Phi^-\\rangle \\langle \\Phi^-|\\right), \n",
    "$$  \n",
    "\n",
    "where total phases were omitted in this representation. Without loss of generality, we assume that Alice performs this unilateral rotation, so that the quantum circuit of the protocol up to this point looks like this:\n",
    "\n",
    "<img vspace=\"30\">\n",
    "<img src=\"images/fig03.png\" width=\"1000\" vspace=\"30\">\n",
    "\n",
    "\n",
    "### Step 3: Application of Bilateral CNOT Gate\n",
    "\n",
    "Now that Alice and Bob can prepare the states $\\rho_F$ from the source states $\\rho_\\textrm{S}$, they will use two copies of $\\rho_F$ to increase the entanglement of one of these copies. In the first step, the two parties must apply a so-called bilateral CNOT gate to the two qubit pairs $\\rho_F$.\n",
    "\n",
    "A bilateral CNOT gate means that Alice and Bob each apply a CNOT gate to the two qubits they have locally available. The qubit from the first pair serves as the control qubit in the local CNOT operation, while the corresponding qubit from the second pair serves as the target qubit. Remember: In a CNOT gate, the state of the target qubit is flipped if and only if the control qubit is in the state $|1\\rangle$.\n",
    "\n",
    "If we now consider all four qubits from the two source qubit pairs $\\rho_\\textrm{S}$, the local quantum circuits of Alice and Bob look like this:\n",
    "\n",
    "<img vspace=\"30\">\n",
    "<img src=\"images/fig04.png\" width=\"1000\" vspace=\"30\">\n",
    "\n",
    "Here we have denoted Alice's (Bob's) control qubit $q_{A(B),c}$ and the corresponding target qubit $q_{A(B),t}$.\n",
    "\n",
    "### Step 4: Measurement of Target Qubits\n",
    "\n",
    "In the second step, Alice and Bob measure the target qubits $q_{A,t}$ and $q_{B,t}$ in the computational basis $\\{|0\\rangle, |1\\rangle \\}$:\n",
    "\n",
    "<img vspace=\"30\">\n",
    "<img src=\"images/fig05.png\" width=\"1000\" vspace=\"0\">\n",
    "\n",
    "After performing the measurements, Alice and Bob compare their measurement results. If both measurement results are the same, i.e. the target qubits were projected into state $|00\\rangle$ or $|11\\rangle$, the purification was successful and the control qubits are retained. Otherwise they will be discarded. \n",
    "\n",
    "To see that the purification process is successful exactly in this case, one must calculate the fidelity $F'$ of the resulting state of the control qubits, which we will call $\\rho'_{F'}$,\n",
    "\n",
    "$$\n",
    "F' = \\langle \\Phi^+ | \\rho'_{F'}| \\Phi^+ \\rangle.\n",
    "$$\n",
    "\n",
    "Note that the state $\\rho'_{F'}$ does not have the same form as $\\rho_{F}$, which is why we chose this notation. We will not elaborate on the calculation of $\\rho'_{F'}$ in this introduction, as it is not required for the implementation in Qiskit, but only provide the result for $F'$[1],\n",
    "\n",
    "$$\n",
    "F' \\equiv F'(F) = \\langle \\Phi^+ | \\rho'_{F'}| \\Phi^+ \\rangle = \\frac{F^2 + \\frac{1}{9}(1-F)^2}{F^2 + \\frac{2}{3}F(1-F) + \\frac{5}{9}(1-F)^2}.\n",
    "$$\n",
    "\n",
    "Since the function $F'(F)$ is required at several points in this tutorial, we define the corresponding function `output_fidelity`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def output_fidelity(F_in):\n",
    "    \"\"\" Returns output fidelity of BBPSSW protocol.\n",
    "    \n",
    "    Keyword arguments:\n",
    "    F_in -- input fidelity\n",
    "    \"\"\"\n",
    "    F_out = (np.square(F_in) + np.square(1-F_in)/9) \\\n",
    "            / (np.square(F_in) + 2*F_in*(1-F_in)/3 + 5*np.square(1-F_in)/9)\n",
    " \n",
    "    return F_out"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The output fidelity $F'$ is larger than the input fidelity for $1/2 < F < 1$, which you can see very well when plotting the function $F'(F)$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "F_in = np.arange(0.5, 1.0, 0.01)\n",
    "plt.plot(F_in, F_in, 'k', label = \"$F' = F$\")\n",
    "plt.plot(F_in, output_fidelity(F_in), label = \"$F'(F)$\")\n",
    "plt.legend()\n",
    "plt.xlabel('Input Fidelity $F$')\n",
    "plt.ylabel(\"Output Fidelity $F'$\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, if the input fidelity $F$ of the source state $\\rho_\\textrm{S}$ is larger than 1/2, Alice and Bob can distill entanglement from two copies of that state and increase the fidelity of one of those copies. The two parties can then use two qubit pairs whose fidelity has already been increased to create another qubit pair with even higher fidelity, and so on. By iterating the procedure of the BBPSSW protocol described above, the fidelity of a two-qubit state can thus be brought arbitrarily close to 1, which in turn corresponds to a maximally-entangled Bell state.\n",
    "\n",
    "### Step 5: Restoration of Werner State\n",
    "\n",
    "Note that the state $\\rho'_{F'}$ does not usually have the form of $\\rho_{F}$, i.e. it can contain off-diagonal elements in the Bell basis and the components of $|\\Phi^{-}\\rangle$, $|\\Psi^{+}\\rangle$ and $|\\Psi^{-}\\rangle$ differ in size. But in order to perform a subsequent purification round, we must bring $\\rho'_{F'}$ into the form of $\\rho_{F}$, i.e. we must create something like a state $\\rho_{F'}$,\n",
    "\n",
    "$$\n",
    "\\rho_{F'} = F' |\\Phi^+\\rangle \\langle \\Phi^+| + \\frac{1-F'}{3}\\left(|\\Psi^+\\rangle \\langle \\Psi^+|+|\\Psi^-\\rangle \\langle \\Psi^-|+|\\Phi^-\\rangle \\langle \\Phi^-|\\right).  \n",
    "$$  \n",
    "\n",
    "This can be achieved by performing the above steps backwards again. In the state $\\rho'_{F'}$, the components of $|\\Psi^-\\rangle$ and $|\\Phi^+\\rangle$ are first swapped again by a unilateral $Y$ gate (see step 2). Then another twirl operation is performed to create a Werner state $W_{F'}$,\n",
    "\n",
    "$$\n",
    "W_{F'} = F' |\\Psi^-\\rangle \\langle \\Psi^-| + \\frac{1-F'}{3}\\left(|\\Psi^+\\rangle \\langle \\Psi^+|+|\\Phi^+\\rangle \\langle \\Phi^+|+|\\Phi^-\\rangle \\langle \\Phi^-|\\right).  \n",
    "$$\n",
    "\n",
    "This state $W_{F'}$ can now serve as the initial state for a new round of entanglement purification, , since we can easily create the state $\\rho_{F'}$ from it.\n",
    "\n",
    "In the following implementation in Qiskit we will shorten the purification protocol and omit this fifth step. We only want to do a single purification round here and therefore do not convert the output state $\\rho'_{F'}$ of the control qubits back to the Werner form $W_{F'}$ after the target qubits have been measured. Instead, we will directly evaluate the fidelity of the control qubits to the Bell state $| \\Phi^+ \\rangle$,\n",
    "\n",
    "$$\n",
    "\\langle \\Phi^+ | \\rho'_{F'}| \\Phi^+ \\rangle = F' = \\langle \\Psi^- | W_{F'}| \\Psi^- \\rangle,\n",
    "$$\n",
    "\n",
    "because we know the twirl operation leaves the $|\\Psi^- \\rangle$ component unchanged.\n",
    "\n",
    "Herewith we finish the introductory part and turn to the implementation in Qiskit. We will explain in detail how steps 1 - 4 described above can be implemented in Qiskit to determine the state $\\rho'_{F'}$ and calcutate the final $F'$ fidelity."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Part II: Implementation of Quantum Circuit for BBPSSW Protocol in Qiskit\n",
    "***\n",
    "\n",
    "## 2.1 Basic Code Imports\n",
    "\n",
    "To use Qiskit in this initial part of the tutorial in which we create a quantum circuit for the BBPSSW protocol, we must first import the necessary base code from the Qiskit module:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import Qiskit\n",
    "from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister\n",
    "\n",
    "# Import Parameter class\n",
    "from qiskit.circuit import Parameter"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also need the external library `deepcopy` and use the constant $\\pi$ from the `math` library:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from copy import deepcopy\n",
    "from math import pi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.2 Definitions\n",
    "\n",
    "As described in **Part I**, to demonstrate entanglement purification, the fidelity of quantum states towards a maximally-entangled Bell state must be calculated. This calculation is later performed with an integrated Qiskit function that returns the fidelity between two states when specified as numpy arrays. \n",
    "\n",
    "We therefore define the four Bell states in the basis $\\{|00\\rangle, |01\\rangle, |10\\rangle, |11\\rangle\\}$ as numpy arrays:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "phi_plus = np.array([1, 0, 0, 1])/np.sqrt(2)   # | Phi^+ >\n",
    "phi_minus = np.array([1, 0, 0, -1])/np.sqrt(2) # | Phi^- >\n",
    "psi_plus = np.array([0, 1, 1, 0])/np.sqrt(2)   # | Psi^+ >\n",
    "psi_minus = np.array([0, 1, -1, 0])/np.sqrt(2) # | Psi^- >"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will also have to check later how well we can initialize a Werner state $W_F$. For this we define a function `Werner_state` which returns a numpy array corresponding to a Werner state with a given fidelity $F$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Werner_state(F):\n",
    "    \"\"\"Returns Werner state of fidelity F.\n",
    "    \n",
    "    Keyword arguments:\n",
    "    F -- fidelity of Werner state, range [0, 1]\n",
    "    \"\"\"\n",
    "    if F < 0 or F > 1:\n",
    "        raise Exception('Fidelity must be between 0 and 1.')\n",
    "    \n",
    "    state = F*np.outer(psi_minus, psi_minus.conj().T) \\\n",
    "            + (1 - F)/3*(np.outer(phi_plus, phi_plus.conj().T) \\\n",
    "                        + np.outer(phi_minus, phi_minus.conj().T) \\\n",
    "                        + np.outer(psi_plus, psi_plus.conj().T))\n",
    "    \n",
    "    return state"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.3 Quantum Register and Quantum Circuit\n",
    "\n",
    "This section describes how to implement the quantum circuit for entanglement purififcation in Qiskit. For a better understanding, we divide the quantum circuit into its individual parts according to the individual steps in **Part I**. This means that we first simulate the source of the entangled qubits to create two qubits in the source state $\\rho_\\textrm{S}$. We then integrate steps 1 (random bilateral rotation) to 3 (bilateral CNOT gate) individually. Step 4 (measurement of target qubits) will then be dealt with separately later, as the implementation is a little trickier here.\n",
    "\n",
    "\n",
    "### 2.3.1 Quantum Register\n",
    "\n",
    "The BBPSSW protocol works with two entangled qubit pairs. In order to implement this protocol, we therefore need a quantum register consisting of four qubits on which our quantum circuit to be defined runs.\n",
    "\n",
    "In Qiskit a quantum register is an instance of the class `QuantumRegister` and we can create a register with four qubits in the following way:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qr = QuantumRegister(4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.3.2 Circuit for Initialization of Source State\n",
    "\n",
    "First of all we have to create qubit pairs in the state $\\rho_\\textrm{S}$. Our goal at the beginning is therefore to simulate the source $S$, which produces entangled qubits in the state $\\rho_\\textrm{S}$, which have an initial fidelity $F > 1/2$. Remember that such a high fidelity is required for the purification protocol to work.\n",
    "\n",
    "To achieve this, we will implement a quantum circuit that creates the following state $|\\psi_\\textrm{S} \\rangle$ between two qubits,\n",
    "\n",
    "$$\n",
    "|\\psi_\\textrm{S}\\rangle = \\cos\\left(\\frac{\\theta}{2}\\right) |01\\rangle - \\sin\\left(\\frac{\\theta}{2}\\right) |10\\rangle ,\n",
    "$$\n",
    "\n",
    "where the two amplitudes are parametrized by an angle $\\theta$ ($- \\pi/2 < \\theta < \\pi/2$). From the state vector $|\\psi_\\textrm{S} \\rangle$ we can then calculate the density matrix $\\rho_\\textrm{S} = |\\psi_\\textrm{S}\\rangle\\langle \\psi_\\textrm{S}|$,\n",
    "\n",
    "$$\n",
    "\\rho_\\textrm{S} = |\\psi_\\textrm{S}\\rangle\\langle \\psi_\\textrm{S} |\n",
    "=\n",
    "\\cos^2\\left(\\frac{\\theta}{2}\\right) |01\\rangle \\langle 01|\n",
    "-\n",
    "\\cos\\left(\\frac{\\theta}{2}\\right)\\sin\\left(\\frac{\\theta}{2}\\right) \n",
    "\\left( |01\\rangle\\langle 10| + |10\\rangle\\langle 01| \\right)\n",
    "+\n",
    "\\sin^2\\left(\\frac{\\theta}{2}\\right) |10\\rangle \\langle 10|.\n",
    "$$\n",
    "\n",
    "The fidelity of $\\rho_\\textrm{S}$ towards the Bell state $| \\Psi^- \\rangle$, which corresponds to the initial fidelity $F$, is then given by\n",
    "\n",
    "$$\n",
    "F \\equiv F(\\theta) = \\langle \\Psi^- | \\rho_\\textrm{S} | \\Psi^- \\rangle  = \\frac{1}{2} \\left( \\sin(\\theta) + 1 \\right).\n",
    "$$\n",
    "\n",
    "We define a function `F_to_theta`, which returns the angle $\\theta$ for a given fidelity $F$, since we will need this calculation several times in the tutorial:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def F_to_theta(F):\n",
    "    \"\"\"Returns the angle theta corresponding to an input fidelity F.\n",
    "    \n",
    "    Keyword arguments:\n",
    "    F -- fidelity between 0 and 1\n",
    "    \"\"\"\n",
    "    if F < 0 or F > 1:\n",
    "        raise Exception('Fidelity must be between 0 and 1.')\n",
    "    else:\n",
    "        theta = np.arcsin(2*F - 1)\n",
    "    \n",
    "    return theta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In Qiskit quantum circuits are objects of the `QuantumCircuit` class. We can create a quantum circuit on our quantum register `qr` by creating an instance of the `QuantumCircuit` class in the following way:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "source_circ = QuantumCircuit(qr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The quantum circuit `source_circ` shall finally generate two pairs of qubits in the state $|\\psi_\\textrm{S}\\rangle$ in the quantum register `qr`, but is so far an empty quantum circuit, because we have not introduced quantum gates yet. With the method `draw` quantum circuits can be visualized nicely and we can get a better picture of what our circuit does:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "source_circ.draw(output='mpl')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see how the quantum register `qr` is constructed from the four qubits $q0_0$ to $q0_3$. (q0 is actually the name of the register itself.) Note that the qubit label starts with index 0 according to the indexing in Python.\n",
    "\n",
    "We will create a source state $|\\psi_\\textrm{S}\\rangle$ between qubits $q0_0$ (qubit 0) and $q0_1$ (qubit 1) and between qubits $q0_2$ (qubit 2) and $q0_3$ (qubit 3). Qubit 0 and qubit 1 should serve as control qubits, while the target qubits are qubit 2 and qubit 3. If we compare this with the representation in the introductory **Part I**, we find the following correspondences:\n",
    "\n",
    "| Qubit  | Notation  | Notation in Qiskit  |\n",
    "|:-:|:-:|:-:|\n",
    "| Control qubit Alice  |  $q_{A,c}$ | $q0_0$  | \n",
    "| Control qubit Bob  |  $q_{B,c}$ | $q0_1$  |\n",
    "| Target qubit Alice  |  $q_{A,t}$ | $q0_2$  |\n",
    "| Target qubit Bob  |  $q_{B,t}$ | $q0_3$  |\n",
    "\n",
    "The following shows how to set the control pair to $|\\psi_\\textrm{S}\\rangle$. The implementation for the target pair is left as an exercise.\n",
    "\n",
    "By default, all qubits are initialized to the state $|0\\rangle$, i.e. the state $|\\psi_\\textrm{control}\\rangle$ of the control pair is \n",
    "\n",
    "$$\n",
    "|\\psi_\\textrm{control}\\rangle = |0\\rangle_{0} |0\\rangle_{1},\n",
    "$$\n",
    "\n",
    "where the index is to denote the corresponding qubit. As first step we flip qubit 1 using a Pauli-$X$ gate,\n",
    "\n",
    "$$\n",
    "X = \n",
    "\\begin{pmatrix}\n",
    "0 & 1 \\\\\n",
    "1 & 0 \\\\\n",
    "\\end{pmatrix}.\n",
    "$$\n",
    "\n",
    "In Qiskit a Pauli-$X$ gate can be implemented in a quantum circuit with the method `x` of the class`QuantumCircuit`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "source_circ.x(qr[1]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Every qubit $i$ in the register can be called with `qr[i]`. So `x(qr[1])` implements an $X$ gate on qubit 1. If we draw the quantum circuit again, we see that a single-qubit gate $X$ is applied to qubit $q0_1$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "source_circ.draw(output='mpl')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After applying the $X$ gate, the control pair is in the state\n",
    "\n",
    "$$\n",
    "|\\psi_\\textrm{control}\\rangle = |0\\rangle_{0} |1\\rangle_{1}.\n",
    "$$\n",
    "\n",
    "We then apply a specific single-qubit gate from the general single-qubit rotations $u(\\theta, \\phi, \\lambda)$ to qubit 0,\n",
    "\n",
    "$$\n",
    "u(\\theta, \\phi , \\lambda) = \\begin{pmatrix}\n",
    "\\cos(\\theta/2) & -e^{i\\lambda}\\sin(\\theta/2) \\\\\n",
    "e^{i\\phi}\\sin(\\theta/2) & e^{i\\lambda+i\\phi}\\cos(\\theta/2) \n",
    "\\end{pmatrix}.\n",
    "$$\n",
    "\n",
    "We set the parameters $\\phi = \\pi$ and $\\lambda = 0$,\n",
    "\n",
    "$$\n",
    "u(\\theta, \\phi = \\pi , \\lambda = 0) = \\begin{pmatrix}\n",
    "\\cos(\\theta/2) & -\\sin(\\theta/2) \\\\\n",
    "-\\sin(\\theta/2) & -\\cos(\\theta/2) \\\\\n",
    "\\end{pmatrix},\n",
    "$$\n",
    "\n",
    "to bring the control qubits into the following state\n",
    "\n",
    "$$\n",
    "|\\psi_\\textrm{control}\\rangle = \\cos\\left(\\frac{\\theta}{2}\\right) |0\\rangle_{0} |1\\rangle_{1} - \\sin\\left(\\frac{\\theta}{2}\\right) |1\\rangle_{0} |1\\rangle_{1}.\n",
    "$$\n",
    "\n",
    "Such a one-qubit rotation is realized with the method `u`, to which we pass the rotation parameters as input arguments:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a parameter theta\n",
    "theta = Parameter('theta')\n",
    "\n",
    "# Implement a parametrized single-qubit gate u on qubit 0\n",
    "source_circ.u(theta, pi, 0, qr[0]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Later in this tutorial, we will perform entanglement purification for different input fidelities $F$, i.e. for different values of $\\theta$. Therefore, we have defined a parameter called `theta`, which is an instance of the class `Parameter` and to which a value can be assigned when the quantum circuit is executed. The input parameters of the method `u` are the angles $\\theta$, $\\phi$ and $\\lambda$, and the qubit to which the gate is applied. If we draw the circuit again, we now see another single-qubit operation $U$ applied to qubit $q0_0$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "source_circ.draw(output='mpl')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we have to use a CNOT gate between qubit 0 and qubit 1, where qubit 1 is flipped if qubit 0 is in the state $|1\\rangle$ to create the desired source state $|\\psi_\\textrm{S}\\rangle$,\n",
    "\n",
    "$$\n",
    "|\\psi_\\textrm{control}\\rangle =\n",
    "|\\psi_\\textrm{S}\\rangle =\n",
    "\\cos\\left(\\frac{\\theta}{2}\\right) |0\\rangle_{0} |1\\rangle_{1} - \\sin\\left(\\frac{\\theta}{2}\\right) |1\\rangle_{0} |0\\rangle_{1}.\n",
    "$$\n",
    "\n",
    "In Qiskit a CNOT gate can be added to a quantum circuit using the method `cx`, which uses the control and the target qubit as input arguments:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Implement a CNOT between qubit 0 as control qubit and qubit 1 as target qubit\n",
    "source_circ.cx(qr[0], qr[1]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(To avoid any misunderstandings here due to the term control and target qubit: The CNOT used here serves to generate the source state $|\\psi_\\textrm{S}\\rangle$ and has nothing to do with the bilateral CNOT of the BBPSSW protocol.)\n",
    "\n",
    "The complete quantum circuit to simulate a source of entangled qubit pairs and generate the state $|\\psi_\\textrm{S}\\rangle$ between qubits $q0_0$ and $q0_1$ looks like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "source_circ.draw(output='mpl')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since we need two copies of the state $\\rho_\\textrm{S}$ to perform entanglement purification, we also need to set qubits 2 and 3 to the source state $|\\psi_\\textrm{S}\\rangle$,\n",
    "\n",
    "$$\n",
    "|\\psi_\\textrm{target}\\rangle =\n",
    "|\\psi_\\textrm{S}\\rangle =\n",
    "\\cos\\left(\\frac{\\theta}{2}\\right) |0\\rangle_{2} |1\\rangle_{3} - \\sin\\left(\\frac{\\theta}{2}\\right) |1\\rangle_{2} |0\\rangle_{3}.\n",
    "$$\n",
    "\n",
    "We call this state $|\\psi_\\textrm{target}\\rangle$ because the two qubits are the target pair of the BBPSSW protocol.\n",
    "\n",
    "> **Exercise 1**  \n",
    "Add the necessary quantum gates to the circuit `source_circ` so that the state $|\\psi_\\textrm{S}\\rangle$ is generated between qubits 2 and 3. The implementation can be done in analogy to the above steps in which this state was generated between the control qubits."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "###\n",
    "### Enter your code below\n",
    "###\n",
    "\n",
    "# Implement a Pauli-X gate on qubit 3\n",
    "source_circ.x(qr[3])\n",
    "\n",
    "# Implement a single-qubit gate u on qubit 2\n",
    "source_circ.u(theta, pi , 0, qr[2])\n",
    "\n",
    "# Implement a CNOT between qubit 2 as control qubit and qubit 3 as target qubit\n",
    "source_circ.cx(qr[2], qr[3]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can check your solution for Exercise 1 by drawing the circuit `source_circ` and comparing the two parts that create the state $|\\psi_\\textrm{S}\\rangle$ between the control and the target qubits:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "source_circ.draw(output='mpl')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A single-qubit rotation $U(\\textrm{theta}, \\pi, 0)$ should now also be applied to qubit $q0_2$, an $X$ gate should also be applied to qubit $q0_3$ and a CNOT should be applied between the two qubits.\n",
    "\n",
    "To clearly define which qubits are the control pair and which qubits are the target pair of the BBPSSW protocol, we define the two arrays `control_qubits` and `target_qubits`, which contain the corresponding qubits:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "control_qubits = [qr[0], qr[1]]\n",
    "target_qubits = [qr[2], qr[3]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.3.3 Implementation of Twirl Operation by Random Bilateral Rotation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To implement the twirl operation, we create another instance of a quantum circuit called `twirl_circ` that performs the random bilateral rotation explained in step 1 of **Sec. 1.2**:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "twirl_circ = QuantumCircuit(qr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remember: To implement the twirl, Alice and Bob randomly select an element from the set of 12 bilateral operations $\\{U_i\\}$ and apply the corresponding rotations locally to their qubit.\n",
    "\n",
    "To implement an operation randomly selected from $\\{U_i\\}$, we define an auxiliary function `bilateral_rotation`. This function has the following input parameters:\n",
    "* `circuit`: The instance of a quantum circuit to which the rotations are to be added.\n",
    "* `qubits`: A list of qubits, each of which is rotated.\n",
    "* `rot_idx`: An integer number that numbers the elements of the set $\\{U_i\\}$. For example, `rot_idx = 0` corresponds to the identity operation, `rot_idx = 1` corresponds $B_x B_x$, and so on."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def bilateral_rotation(circuit, qubits, rot_idx):\n",
    "    \"\"\"Adds a bilateral rotation to a quantum circuit.\n",
    "    \n",
    "    Keyword arguments:\n",
    "    circuit --  quantum circuit to which the rotations are added\n",
    "    qubits -- alist of qubits, each of which is rotated\n",
    "    rot_idx -- integer that numbers the rotations from the set of all rotations, ranges from 0 to 11\n",
    "    \"\"\"\n",
    "    for qubit in qubits:\n",
    "        if rot_idx == 0:\n",
    "            pass\n",
    "        \n",
    "        elif rot_idx == 1:\n",
    "            circuit.rx(pi/2, qubit)\n",
    "            circuit.rx(pi/2, qubit)\n",
    "            \n",
    "        elif rot_idx == 2:\n",
    "            circuit.ry(pi/2, qubit)\n",
    "            circuit.ry(pi/2, qubit)\n",
    "            \n",
    "        elif rot_idx == 3:\n",
    "            circuit.rz(pi/2, qubit)\n",
    "            circuit.rz(pi/2, qubit)\n",
    "            \n",
    "        elif rot_idx == 4:\n",
    "            circuit.rx(pi/2, qubit)\n",
    "            circuit.ry(pi/2, qubit)\n",
    "            \n",
    "        elif rot_idx == 5:\n",
    "            circuit.ry(pi/2, qubit)\n",
    "            circuit.rz(pi/2, qubit)\n",
    "            \n",
    "        elif rot_idx == 6:\n",
    "            circuit.rz(pi/2, qubit)\n",
    "            circuit.rx(pi/2, qubit)\n",
    "            \n",
    "        elif rot_idx == 7:\n",
    "            circuit.ry(pi/2, qubit)\n",
    "            circuit.rx(pi/2, qubit)\n",
    "            \n",
    "        elif rot_idx == 8:\n",
    "            circuit.rx(pi/2, qubit)\n",
    "            circuit.ry(pi/2, qubit)\n",
    "            circuit.rx(pi/2, qubit)\n",
    "            circuit.ry(pi/2, qubit)\n",
    "            \n",
    "        elif rot_idx == 9:\n",
    "            circuit.ry(pi/2, qubit)\n",
    "            circuit.rz(pi/2, qubit)\n",
    "            circuit.ry(pi/2, qubit)\n",
    "            circuit.rz(pi/2, qubit)\n",
    "            \n",
    "        elif rot_idx == 10:\n",
    "            circuit.rz(pi/2, qubit)\n",
    "            circuit.rx(pi/2, qubit)\n",
    "            circuit.rz(pi/2, qubit)\n",
    "            circuit.rx(pi/2, qubit)\n",
    "            \n",
    "        elif rot_idx == 11:\n",
    "            circuit.ry(pi/2, qubit)\n",
    "            circuit.rx(pi/2, qubit)\n",
    "            circuit.ry(pi/2, qubit)\n",
    "            circuit.rx(pi/2, qubit)\n",
    "            \n",
    "        else:\n",
    "            raise Exception('Input argument rot_idx has to be an integer between 0 and 11.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The single-qubit rotations around the $x$, $y$ and $z$ axes required for the bilateral rotation can be implemented in Qiskit with the following methods:\n",
    "\n",
    "* `rx`: Rotation about the $x$ axis.\n",
    "* `ry`: Rotation about the $y$ axis.\n",
    "* `rz`: Rotation about the $z$ axis.\n",
    "\n",
    "The angle of rotation and the respective qubit must be specified as input arguments.\n",
    "\n",
    "Remember that the bilateral operations $B_k$ from which the elements of $\\{ U_i\\}$ are constructed are rotations around the $k$ axis by an angle of $\\pi/2$ and that the rotations are applied to both qubits of the entangled pair. With this knowledge it is understandable how the operations in the function `bilateral_rotation` are built up from single-qubit rotations `rx`, `ry` and `rz` by angles of $\\pi/2$.\n",
    "\n",
    "Now we add a random bilateral rotation to the circuit `twirl_circ` with the function `bilateral_rotation`, once for the control qubit pair and once for the target qubit pair:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a random integer between 0 and 11\n",
    "random_idx = np.random.randint(12)\n",
    "\n",
    "# Add a random bilateral rotation to the circuit, applied to qubits 0 and 1 (control qubits)\n",
    "bilateral_rotation(twirl_circ, control_qubits, random_idx)\n",
    "\n",
    "# Add a random bilateral rotation to the circuit, applied to qubits 2 and 3 (target qubits)\n",
    "bilateral_rotation(twirl_circ, target_qubits, random_idx)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The circuit `twirl_circ` thus has the following structure:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "twirl_circ.draw(output='mpl')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If this operation is applied to a large number of qubit pairs in the state $|\\psi_\\textrm{S}\\rangle$, the average output state is a Werner state $W_F$ (see step 1, **Sec. 1.2**),\n",
    "\n",
    "$$\n",
    "W_F = F |\\Psi^-\\rangle \\langle \\Psi^-| + \\frac{1-F}{3}\\left(|\\Psi^+\\rangle \\langle \\Psi^+|+|\\Phi^+\\rangle \\langle \\Phi^+|+|\\Phi^-\\rangle \\langle \\Phi^-|\\right).  \n",
    "$$\n",
    "\n",
    "This means that on average the control and target pairs are brought into a Werner state after the circuit `twirl_circ` has been applied.\n",
    "\n",
    "### 2.3.4 Exchange of $|\\Psi^-\\rangle$ and $|\\Phi^+\\rangle$ Components\n",
    "\n",
    "The next step is that Alice applies a unilateral $\\pi$ rotation around the $y$ axis to her qubits, so that the components of $|\\Psi^-\\rangle$ and $|\\Phi^+\\rangle$ are exchanged in the control and the target pair, which leads to the transformation of a Werner state $W_F$ into the state $\\rho_F$ (see step 2, **Sec 1.2**),\n",
    "\n",
    "$$\n",
    "\\rho_F = F |\\Phi^+\\rangle \\langle \\Phi^+| + \\frac{1-F}{3}\\left(|\\Psi^+\\rangle \\langle \\Psi^+|+|\\Psi^-\\rangle \\langle \\Psi^-|+|\\Phi^-\\rangle \\langle \\Phi^-|\\right). \n",
    "$$  \n",
    "\n",
    "This means that the rotations are applied to one qubit of the control pair and one qubit of the target pair. Such a rotation corresponds to a Pauli-$Y$ gate, i.e. a bit flip and a phase flip gate,\n",
    "\n",
    "$$\n",
    "Y = \n",
    "\\begin{pmatrix}\n",
    "0 & -i \\\\\n",
    "i & 0 \\\\\n",
    "\\end{pmatrix}.\n",
    "$$\n",
    "\n",
    "With Qiskit, such a rotation can be applied to a qubit with the method `y` on a circuit. We create a new quantum circuit called `purification_circ` and add the $Y$ gates to the corresponding qubits:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "purification_circ = QuantumCircuit(qr)\n",
    "\n",
    "# Implement a Pauli-Y gate on Alice's control qubit\n",
    "purification_circ.y(control_qubits[0])\n",
    "\n",
    "# Implement a Pauli-Y gate on Alice's target qubit\n",
    "purification_circ.y(target_qubits[0]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We call this quantum circuit `purification_circ`, because this operation is actually the first part of the BBPSSW protocol, which starts from qubit pairs in the Werner state $W_F$. (In the next section we therefore add the bilateral CNOT gate to this circuit.)\n",
    "\n",
    "In this quantum circuit, the single-qubit operations $Y$ are applied to the qubits $q0_0$ and $q0_2$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "purification_circ.draw(output='mpl')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the quantum circuit `purification_circ` is applied to two qubit pairs in the Werner state $W_F$, this results in two qubit pairs in the state $\\rho_F$.\n",
    "\n",
    "### 2.3.5 Application of Bilateral CNOT Gate\n",
    "\n",
    "The next step of the BBPSSW protocol is that Alice and Bob apply a bilateral CNOT gate. This means that Alice and Bob each perform a CNOT gate between the two qubits they hold locally (step 3, **Sec 1.2**).\n",
    "\n",
    "> **Exercise 2**  \n",
    "We have already defined the arrays `control_qubits` and `target_qubits` above, which contain the qubits of the control and the target pair. The first element of these arrays is Alice's qubit and the second is Bob's.\n",
    ">  \n",
    "Add the bilateral CNOT operation to the quantum circuit `purification_circ` by applying CNOT gates between the corresponding elements of the `control_qubits` and `target_qubits` arrays. I.e. a CNOT between Alice's qubits and a CNOT between Bob's qubits. Remember that in Qiskit a CNOT gate is implemented using the method `cx`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "###\n",
    "### Enter your code below\n",
    "###\n",
    "\n",
    "for i in range(2):\n",
    "    purification_circ.cx(control_qubits[i], target_qubits[i])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you now draw the circuit `purification_circ`, you should see two CNOT gates in addition to the two $Y$ gates, one between qubits $q0_0$ and $q0_2$, and one between qubits $q0_1$ and $q0_3$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "purification_circ.draw(output='mpl')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the subsequent step of the BBPSSW protocol (step 4, **Sec 1.2**) the target qubits (qubits 2 and 3) should be measured in the computational basis. The measurement results must then be compared a successful round of entanglement purification.\n",
    "\n",
    "We will implement this step during the process of quantum state tomography, but this is described and explained below (**Sec 2.4**). Briefly summarized: It will work in such a way that we will only perform state tomography of the control qubits if the results of the measurements of the target qubits are the same. \n",
    "\n",
    "We would like to point out once again that after the bilateral CNOT gate and the measurement of the target qubits, the control qubit pair must be returned to Werner form in order to perform a subsequent purification round (step 5, **Sec. 1.2**). As already explained, this step is not necessary to prove an increase in fidelity and can therefore be omitted for our purposes.\n",
    "\n",
    "### 2.3.6 Concatenation of Quantum Circuits\n",
    "\n",
    "We have divided the entire quantum circuit for the BBPSSW protocol into three parts for the initialization of the entangled qubit pairs (`source_circ`), the twirl operation (`twirl_circ`) and the purification step itself (`purification_circ`). \n",
    "\n",
    "(Note that this was only done for educational reasons. We could have built all gates into a single quantum circuit from the beginning.)\n",
    "\n",
    "In Qiskit it is possible to concatenate quantum circuits simply by adding. Therefore, the entire quantum circuit for the purification experiment can be constructed from the three circuits `source_circ`, `twirl_circ`, and `purification_circ` as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "experiment_circ = source_circ + twirl_circ + purification_circ"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Overall, the complete circuit for the BBPSSW protocol has the following structure:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "experiment_circ.draw(output='mpl')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now compare this with the quantum circuits from the introductory part (**Part I**). The gates up to the first CNOTs are not part of Alice's and Bob's circuits presented in the introduction. This is because these gates are needed to first simulate the source that creates the source states $\\rho_\\textrm{S}$. Then we recognize the single-qubit parts $U_{i,A(B)}$ of the bilateral rotation $U_i$, which of course depend here on the random index `random_idx` chosen above. (Only in the case of `random_idx = 0` do you actually see no parts of the bilateral rotation, as this corresponds to the identity operation.) Then the unilateral rotations Alice performs are implemented on the qubits $q0_0$ and $q0_2$, followed by the bilateral CNOT.\n",
    "\n",
    "## 2.4 (Conditional) Quantum State Tomography\n",
    "\n",
    "Finally, we want to find out the state $\\rho'_{F'}$ of the control qubit pair, i.e. the state of the qubits 0 and 1, so that we can calculate its fidelity $F'$ (see step 4, **Sec 1.2**),\n",
    "\n",
    "$$\n",
    "F' = \\langle \\Phi^+ | \\rho'_{F'} | \\Phi^+ \\rangle.\n",
    "$$\n",
    "\n",
    "To reconstruct the density matrix $\\rho'_{F'}$, we can use methods of **quantum tomography**. The library for tomography is part of the Qiskit element **Ignis**, which is a framework for handling noise in quantum circuits. Qiskit Ignis offers methods for the implementation of state and process tomography. For the purification protocol we only have to implement **quantum state tomography** and will therefore concentrate on it in this tutorial.\n",
    "\n",
    "The functions and classes necessary for performing quantum state tomography can be imported from the tomography library:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiskit.ignis.verification.tomography import state_tomography_circuits, StateTomographyFitter"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to perform quantum state tomography, we must first generate tomography circuits with the function `state_tomography_circuits`. From the measurement results of these tomography circuits we can then reconstruct the state $\\rho'_{F'}$ with methods of the fitter class `StateTomographyFitter`.\n",
    "\n",
    "The tomography circuits themselves consist of the original circuit to which specific gates and qubit measurements are attached. The function `state_tomography_circuits` therefore has the following input arguments:\n",
    "\n",
    "* `circuit`: Quantum circuit on which the tomography is performed.\n",
    "* `measured_qubits`: List of qubits whose quantum state is being reconstructed.\n",
    "\n",
    "In our case we want to reconstruct the state of the control qubits after executing the circuit `experiment_circ`. The desired state tomography circuits can therefore be generated as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tomography_circs = state_tomography_circuits(experiment_circ, control_qubits)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`tomography_circs` is a list of quantum circuits and from the measurement results of all these circuits the state $\\rho'_{F'}$ can be reconstructed.\n",
    "\n",
    "As an example, the first tomography circuit looks like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tomography_circs[0].draw(output='mpl')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can recognize the circuit `experiment_circ` up to the bilateral CNOT gate and see that two additional Hadamard gates are applied to the qubits $q0_0$ and $q0_1$. The control qubits are then measured in the computational basis and the measurement results are stored in a classical register $c0$. This classical register was also added by the function `state_tomography_circuits` and its size corresponds to the number of qubits whose quantum state is to be reconstructed. \n",
    "\n",
    "From the execution results of `tomography_circs` we would reconstruct the density matrix of the control qubits if no measurement was performed on the target qubits. Thus, one can and should rightly remark that **Step 4 (Measurement of Target Qubits)** of the BBPSSW protocol is missing in the quantum circuit `experiment_circ`. Remember that purification only worked if the target qubits were measured and the measurements gave the same results. This is a crucial step, which must, of course, also be taken into account here! For this purpose, we will perform **conditional quantum state tomography** and adapt a strategy presented in [6] to our needs.\n",
    "\n",
    "By conditional state tomography, we mean here that the tomography is only performed if the target qubits were previously measured either both in the state $|0\\rangle$ or both in the state $|1\\rangle$. This corresponds exactly to the discarding of all failed purifications. Therefore we have to add the measurements of the target qubits before building the tomography circuits with the function `state_tomography_circuits`. \n",
    "\n",
    "Now an instance of the class `StateTomographyFitter`, from which we can finally reconstruct a density matrix, is constructed from a set of tomography circuits and the corresponding execution results. The constructor of the class `StateTomographyFitter` needs as input argument a set of quantum circuits in the form of the tomography circuits `tomography_circs`, which contain exactly one classical register $c0$ consisting of two classical bits. However, if we first add the additional measurements of the target qubits, we introduce another classical register to store the measurement results. And in this case we can no longer construct the fitter directly, since the tomography circuits would now contain two classical registers.\n",
    "\n",
    "What we will do instead is to make a copy of the circuit `experiment_circ`, add the measurements of the target qubits, and create a set of corresponding tomography circuits. These circuits are then actually executed. The results of this execution can then be edited to keep only those cases where the target qubits were measured with the same results. The edited results can then be fed into the fitter along with our original circuits `tomography_circs` to reconstruct the state $\\rho'_{F'}$ of the control qubits.\n",
    "\n",
    "The following shows how to create these ancilla circuits and **Part III** explains in detail how to edit the execution results correctly.\n",
    "\n",
    "The conditional state tomography is thus implemented as follows. First we make a copy of our original circuit `experiment_circ`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ancilla_circ = deepcopy(experiment_circ)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We call this circuit `ancilla_circ` because it helps to generate the information needed to reconstruct $\\rho'_{F'}$.\n",
    "\n",
    "We then add a classical register to `ancilla_circ` which is an instance of the class `ClassicalRegister` to store the results of the measurement of the target qubits. We need a two-bit register here:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a classical register cr consisting of two bits\n",
    "cr = ClassicalRegister(2)\n",
    "\n",
    "# Add the classical register to the circuit\n",
    "ancilla_circ.add_register(cr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we must actually add the measurements of the target qubits. A measurement process can be added to a quantum circuit using the method `measure`. The method requires two input arguments, namely the qubit to be measured and the bit in a classical register for storing the measurement result.\n",
    "\n",
    "The desired measurement of the target qubits can therefore be implemented as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Implement a measurement of Alice's target qubit and store the result in bit 0\n",
    "ancilla_circ.measure(target_qubits[0], cr[0])\n",
    "\n",
    "# Implement a measurement of Bob's target qubit and store the result in bit 1\n",
    "ancilla_circ.measure(target_qubits[1], cr[1]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's have a look at the circuit `ancilla_circ`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ancilla_circ.draw(output='mpl')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the circuit diagram we can see that a classical register $c1$ has been added to the circuit `ancilla_circ` compared to the original circuit `experiment_circ` and that the qubits $q0_2$ and $q0_3$ are measured after the bilateral CNOT.\n",
    "\n",
    "For the circuit `ancilla_circ` we also create tomography circuits `ancilla_tomography_circs` using the function `state_tomography_circuits` to reconstruct the state of the control qubits:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ancilla_tomography_circs = state_tomography_circuits(ancilla_circ, control_qubits)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is now these ancilla circuits `ancilla_tomography_circs` that are actually executed. As an example we draw one of the ancilla circuits:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "ancilla_tomography_circs[0].draw(output='mpl')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have already seen that the target qubits $q0_2$ and $q0_3$ are measured in the classical register $c1$. In addition, the results for the tomography of the control qubits $q0_0$ and $q0_1$ are stored in the automatically added classical register $c2$.\n",
    "\n",
    "At this point, we already import a usefull function from the Python toolkit Mitiq [7], which removes all barriers from a quantum circuit:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from mitiq.mitiq_qiskit.conversions import _remove_barriers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The reason for removing the barriers is that circuits as shown above cannot be executed on real IBMQ hardware because no gates are allowed after measurements. Although in the example circuit `ancilla_tomography_circs[0]` above the two Hadamard gates on qubits 0 and 1 \"follow\" measurements on qubits 2 and 3, this circuit is still considered to be not executable on real IBMQ hardware. Therefore we will simply remove the barrier and by this make it executable:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_remove_barriers(ancilla_tomography_circs[0]).draw(output='mpl')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At the end of this Part II, we are at the point where we have generated all the quantum circuits necessary to implement the full BBPSSW protocol. In the following part we can therefore continue with the simulation of the purification experiment and have our implementation run with the simulator integrated in Qiskit."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Part III: Simulation of the BBPSSW Protocol using Qiskit\n",
    "***\n",
    "\n",
    "In this part we will simulate the quantum circuit for the BBPSSW protocol. For this we use a simulator that is integrated in Qiskit. In the first section we simulate the BBPSSW protocol for a single quantum system. To finally demonstrate a successful purification, we calculate an average of many quantum systems to generate approximately the state $\\rho'_{F'}$ and extract the purified fidelity $F'$.\n",
    "\n",
    "To execute a quantum circuit in Qiskit, we must first select a backend on which on which to run the quantum program. Such a backend can either be a simulator or a real IBM Q device.\n",
    "\n",
    "The QASM simulator is part of the Qiskit element **Aer** and to use it we need to import the following libraries:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiskit import Aer, execute"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To select the simulator as backend for the execution of the quantum circuit, we use the method `get_backend`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "backend_QASM = Aer.get_backend('qasm_simulator')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.1 Single Quantum System\n",
    "\n",
    "In the following we will first execute the circuits `tomography_circs` and reconstruct the density matrix for the control qubits. This is done mainly for educational reasons to introduce the functionality of the class `StateTomographyFitter`. Then we run the ancilla circuits `ancilla_tomography_circs` and explain how to edit the execution results to be able to use the fitter.\n",
    "\n",
    "We can execute a quantum circuit with the function `execute`, which requires the following input arguments:\n",
    "* `circuits`: A single quantum circuit or a list of quantum circuits to be executed.\n",
    "* `backend`: The backend to run the circuits.\n",
    "* `shots`: Number of repeated executions of each circuit.\n",
    "* `parameter_binds`: List of parameter-value pairs (given as python dictionary) that bind the parameters to specific values.\n",
    "\n",
    "To simulate the circuits `tomography_circs` with an input fidelity of e.g. $F=0.75$, we can execute the following code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "job = execute(tomography_circs, backend_QASM, shots=1024, parameter_binds=[{theta: F_to_theta(0.75)}])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here the parameter $\\theta$ is bound to the value corresponding to a fidelity of $F = 0.75$ via the function `F_to_theta` defined at the beginning of the tutorial.\n",
    "\n",
    "What is returned after the simulation is an object of the class `AerJob`, which we call `job` and from which we can extract the results of the simulation with the method `result`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "job_results = job.result()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have performed the tomography circuits, we can reconstruct the density matrix of the control qubits using methods from quantum state tomography. All this is implemented in Qiskit by the class `StateTomographyFitter`.\n",
    "\n",
    "Using the object `job_results` we first create a fitter object of the class `StateTomographyFitter`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fitter = StateTomographyFitter(job_results, tomography_circs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The fitter provides a method `fit` that returns the reconstructed density matrix that we call `rho_fit`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "rho_fit = fitter.fit()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The two-qubit state `rho_fit` is returned as a 4x4 numpay array. Although for all further steps we only need this array, it might be helpful to explain a way to visualize quantum states in Qiskit simply for better illustration. We would like to mention the function `plot_state_city` which draws two bar blots of the real and imaginary parts of `rho_fit`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiskit.tools.visualization import plot_state_city"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_state_city(rho_fit)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have seen how to use Qiskit's tools for quantum state tomography, we can proceed and perform conditional quantum state tomography.\n",
    "\n",
    "We start with the same steps as above, i.e. we first execute the ancillas `ancilla_tomography_circs` and retrieve the job results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "job = execute(ancilla_tomography_circs, backend_QASM, shots=1024, parameter_binds=[{theta: F_to_theta(0.75)}])\n",
    "ancilla_results = job.result()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So far the results object `ancilla_results` also contains results of unsuccessful executions (i.e. where the target qubits were measured with different outcomes) and it has a structure from which we cannot directly reconstruct the density matrix of the control qubits. So what we need to do is:\n",
    "\n",
    "1. Discard the results of unsuccessful executions.\n",
    "2. Remove the parts in the object `ancilla_results` that prevent us from using the fitter.\n",
    "\n",
    "To better see excatly what needs to be done, we compare the property `results` of the `job_results` and the `ancilla_results`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "job_results.results[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ancilla_results.results[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The property `results` is a list containing information about the results of each circuit executed.\n",
    "\n",
    "First, the property `counts` must be adjusted so that all events where the target qubits had different outcomes are removed and the counts for the outcomes of the control qubits are corrected accordingly. In the following, we will explain in more detail what we mean by this.\n",
    "\n",
    "We can get the counts as a dictionary using the method `get_counts`. For the first circuit in `ancilla_tomography_circs` the counts dictionary can be created like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ancilla_results.get_counts(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the bit strings, the left-most bit corresponds to the measurement outcome of qubit $q0_3$ and the right-most bit to the outcome of qubit $q0_0$. Therefore, the first two bits are the outcomes of the target qubits and so we must remove all key-value pairs whose first two bits in the key are either `'01'` or `'10'`. In addition, we need to correct the counts of the control qubit outcomes. In the results `job_result`, the counts are:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "job_results.get_counts(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Of course, here we can only have the events `'00'`, `'01'`, `'10'`, and `'11'` because we only measure two qubits. From the counts of the `ancilla_results`, we must therefore add the counts of the successful events which have the same outcome of the control-qubit measurement. For example, if there were events such as `'00 11'` and `'11 11'`, then the counts of these events must be added for a control qubit outcome of `'11'`.\n",
    "\n",
    "To store the edited counts, we first create an empty dictionary `conditional_counts`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "conditional_counts = {}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we loop through all entries in the dictionary `ancilla_results.get_counts(0)` and process the counts exactly as we described above:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for result, count in ancilla_results.get_counts(0).items():\n",
    "    # Get measurement out comes\n",
    "    target_measurement = result.split(' ')[0]\n",
    "    control_measurement = result.split(' ')[1]\n",
    "    \n",
    "    # Keep counts if target qubits were measured '00' or '11'\n",
    "    if target_measurement in ('00', '11'):\n",
    "            \n",
    "            # Add counts to existing events or create new key-value pair\n",
    "            if control_measurement in conditional_counts:\n",
    "                conditional_counts[control_measurement] = conditional_counts[control_measurement] + count\n",
    "            else:\n",
    "                conditional_counts[control_measurement] = count"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The edited counts dictionary, which only contains successful purification events, will then look like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "conditional_counts"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now set the `conditional_counts` as the field `counts` of some other results object, e.g. following the above steps as `counts` of `ancilla_results.results[0].data`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ancilla_results.results[0].data.counts = conditional_counts"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we compare the results `job_results.results.results[0]` and `ancilla_results.results[0]` further, we have to do the following additional steps:\n",
    "* Remove the bits from the classical register $c1$ in `clbit_labels` and in `creg_sizes`. \n",
    "* Set the value of `memory_slots` to 2.\n",
    "\n",
    "In the following we will implement all steps explained above for all results in `ancilla_results`. We store the processed results in an object called `conditional_results`, which we create from the result object `ancilla_results`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "conditional_results = deepcopy(ancilla_results)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The implementation of the code is given below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for ii,_ in enumerate(ancilla_results.results):\n",
    "    counts = ancilla_results.get_counts(ii)\n",
    "    conditional_counts = {}\n",
    "    \n",
    "    # Keep counts of successful purifications\n",
    "    for result, count in counts.items():\n",
    "        target_measurement = result.split(' ')[0]\n",
    "        control_measurement = result.split(' ')[1]\n",
    "        \n",
    "        if target_measurement in ('00', '11'):\n",
    "            if control_measurement in conditional_counts:\n",
    "                conditional_counts[control_measurement] = conditional_counts[control_measurement] + count\n",
    "            else:\n",
    "                conditional_counts[control_measurement] = count\n",
    "    \n",
    "    # Create new counts property from conditional_counts\n",
    "    conditional_results.results[ii].data.counts = conditional_counts\n",
    "    \n",
    "    # Remove classical bits in ancilla register\n",
    "    old_labels = conditional_results.results[ii].header.clbit_labels\n",
    "    conditional_results.results[ii].header.clbit_labels = old_labels[0:2]\n",
    "    \n",
    "    # Remove ancilla register\n",
    "    old_sizes = conditional_results.results[ii].header.creg_sizes\n",
    "    conditional_results.results[ii].header.creg_sizes = [old_sizes[0]]\n",
    "    \n",
    "    # Set memory_slots to 2\n",
    "    conditional_results.results[ii].header.memory_slots = 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we look at the new results `conditional_results`, we see that they now have the same form as the `job_results`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "conditional_results.results[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(To the best knowledge of the author, the property `shots` does not need to be adjusted and the fitter will take care of a possibly smaller number of remaining events.)\n",
    "\n",
    "At this point it is certainly helpful to summarize again which steps we have performed above and why. Remember: All the above steps were necessary because we cannot directly use the results of the ancilla tomography circuit execution to reconstruct the density matrix of the control qubits using the Qiskit functions for quantum state tomography. Therefore we took a small detour and carried out the following steps to perform conditional quantum state tomography:\n",
    "1. Executed the ancilla tomography circuits `ancilla_tomography_circs`.\n",
    "2. Edited the results of this execution so that all failed purifications were discarded.\n",
    "3. Further edited the results to obtain the object `conditional_results`, which has the same structure as a result coming from the execution of the tomography circuits `tomography_circs`.\n",
    "\n",
    "Now we are ready to use the Qiskit features to reconstruct the density matrix of the control qubits, only considering the cases of successful purification.\n",
    "\n",
    "Again with the class `StateTomographyFitter` we reconstruct the density matrix with the method `fit`, but this time we provide the result object `conditional_results` instead of `job_results`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "conditional_rho_fit = StateTomographyFitter(conditional_results, tomography_circs).fit()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The state `conditional_rho_fit` is now the reconstructed state of the control qubits, if only cases of successful purification are considered:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_state_city(conditional_rho_fit)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The state `conditional_rho_fit` depends on the bilateral rotation contained in the twirl circuit `twirl_circ`, since Alice and Bob randomly choose which bilateral rotation is performed at the beginning.\n",
    "\n",
    "Remember that the reason for this was that the BBPSSW protocol works on Werner states $W_F$ and these are arithmetic averages of bilaterally-rotated density matrices,\n",
    "\n",
    "$$\n",
    "W_F = \\frac{1}{12} \\sum_{i=1}^{12} U_i \\rho_\\textrm{S} U_i^\\dagger .\n",
    "$$\n",
    "\n",
    "In practice, Alice and Bob choose a random element from the set $\\{U_i\\}$ of bilateral operations and apply it to a source state $\\rho_\\textrm{S}$. If Alice and Bob apply this random operation to a number $N$ of source states $\\rho_\\textrm{S}$, then the average state for large $N$ is approximately a Werner state,\n",
    "\n",
    "$$\n",
    "W_F \\approx \\frac{1}{N} \\sum_{n=1}^N \\rho_n,\n",
    "$$\n",
    "\n",
    "where $\\rho_n$ denotes the state after the bilateral rotation.\n",
    "\n",
    "The interesting figure is therefore not a single control-qubit state `conditional_rho_fit` after execution of the BBPSSW protocol, but rather an average of these states over many execution runs $N$ with randomly selected operations $U_i$ in the twirl circuit `twirl_circ`. Therefore, in the following section we examine a large number of sources $\\rho_\\textrm{S}$ and calculate the average result of the BBPSSW protocol to demonstrate its working principle.\n",
    "\n",
    "## 3.2 Ensemble Average\n",
    "\n",
    "### 3.2.1 Definition of Auxiliary Functions\n",
    "\n",
    "To make the following section easier to read, we define some helpful functions to implement the tomography for the case of an ensemble. Essentially, these functions summarize the steps from the previous section in a compact way.\n",
    "\n",
    "First we define a function `experiment_circuit`, which has the index of the bilateral rotation as input parameter and returns a corresponding experiment circuit for the whole purification procedure:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def experiment_circuit(index):\n",
    "    \"\"\" Returns experiment circuit corresponding to an index.\n",
    "    \n",
    "    Keyword arguments:\n",
    "    index -- integer between 0 and 11\n",
    "    \n",
    "    Returns:\n",
    "    experiment_circ -- quantum circuit\n",
    "    \"\"\"\n",
    "    if index < 0 or index > 11:\n",
    "        raise Exception(\"Index for bilateral rotation exceeds range of possible values.\")\n",
    "        \n",
    "    # Create twirl circuit from the bilateral rotations\n",
    "    twirl_circ = QuantumCircuit(qr)\n",
    "    bilateral_rotation(twirl_circ, control_qubits, index)\n",
    "    bilateral_rotation(twirl_circ, target_qubits, index)\n",
    "    \n",
    "    experiment_circ = source_circ + twirl_circ + purification_circ\n",
    "    \n",
    "    return experiment_circ"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we define a function `ancilla_circuit` which takes a quantum circuit and returns a corresponding ancilla circuit which is a copy of the original circuit and in which certain qubits are measured:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ancilla_circuit(quantum_circuit, qubits):\n",
    "    \"\"\" Return an ancilla circuit that includes additional measurements.\n",
    "    \n",
    "    Keyword arguments:\n",
    "    quantum_circuit -- single quantum circuit\n",
    "    qubits -- list of qubits to be measured\n",
    "    \n",
    "    Returns:\n",
    "    ancilla_circ -- ancilla quantum circuit including measurement of qubits\n",
    "    \"\"\"\n",
    "    # Copy original circuit\n",
    "    ancilla_circ = deepcopy(quantum_circuit)\n",
    "    \n",
    "    # Add classical register\n",
    "    cr = ClassicalRegister(len(qubits))\n",
    "    ancilla_circ.add_register(cr)\n",
    "    \n",
    "    # Measure qubits\n",
    "    for i, qubit in enumerate(qubits):\n",
    "        ancilla_circ.measure(qubit, cr[i])\n",
    "    \n",
    "    return ancilla_circ"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To perform conditional quantum state tomography, we have seen that the ancilla tomography circuits must be executed, and we need the tomography circuits for the original experiment circuit along with the execution results to fit the state of the control qubits. For this purpose we define a function `conditional_tomography_circuits`, which returns both sets of tomography circuits:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def conditional_tomography_circuits(quantum_circuit):\n",
    "    \"\"\" Returns all tomography circuits for conditional quantum state tomography. \n",
    "    \n",
    "    Keyword arguments:\n",
    "    quantum_circuit -- single quantum circuit\n",
    "    \n",
    "    Returns:\n",
    "    tomography_circs -- tomography circuits for quantum circuit\n",
    "    ancilla_tomography_circs -- ancilla tomography circuits with additional measurements on target qubits\n",
    "    \"\"\"\n",
    "    # Create ancilla circuit\n",
    "    ancilla_circ = ancilla_circuit(quantum_circuit, target_qubits)\n",
    "    \n",
    "    # Create tomography circuits\n",
    "    tomography_circs = state_tomography_circuits(quantum_circuit, control_qubits)\n",
    "    ancilla_tomography_circs = state_tomography_circuits(ancilla_circ, control_qubits)\n",
    "    \n",
    "    return tomography_circs, ancilla_tomography_circs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have also seen how to edit the results of the execution of ancilla tomography circuits so that all shots that do not meet the condition for successful purification are discarded. We summarize those steps in the function `purification_shots`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def purification_shots(job_results):\n",
    "    \"\"\" Keep counts for successful purifications and return corresponding result.\n",
    "    \n",
    "    Keyword arguments:\n",
    "    job_results -- execution results\n",
    "    \n",
    "    Returns:\n",
    "    conditional_results -- results without unsuccessful shots\n",
    "    \"\"\"\n",
    "    # Copy original results\n",
    "    conditional_results = deepcopy(job_results)\n",
    "    \n",
    "    for ii,_ in enumerate(job_results.results):\n",
    "        counts = job_results.get_counts(ii)\n",
    "        conditional_counts = {}\n",
    "        \n",
    "        # Keep counts of successful purifications\n",
    "        for result, count in counts.items():\n",
    "            target_measurement = result.split(' ')[0]\n",
    "            control_measurement = result.split(' ')[1]\n",
    "\n",
    "            if target_measurement in ('00', '11'):\n",
    "                if control_measurement in conditional_counts:\n",
    "                    conditional_counts[control_measurement] = conditional_counts[control_measurement] + count\n",
    "                else:\n",
    "                    conditional_counts[control_measurement] = count\n",
    "        \n",
    "        # Adapt the properties of conditional_results accordingly\n",
    "        conditional_results.results[ii].data.counts = conditional_counts\n",
    "        conditional_results.results[ii].header.clbit_labels = conditional_results.results[ii].header.clbit_labels[0:2]\n",
    "        conditional_results.results[ii].header.creg_sizes = [conditional_results.results[ii].header.creg_sizes[0]]\n",
    "        conditional_results.results[ii].header.memory_slots = 2\n",
    "            \n",
    "    return conditional_results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to reconstruct a density matrix from tomography results, we have seen that an instance of the class `StateTomographyFitter` must be initiated from the tomography circuits and the corresponding execution results. In the following steps of the tutorial it is therefore helpful to save the tomography circuits and the corresponding execution job in some way and we choose a Python dictionary for it. To easily create a dictionary of the desired structure, we define the function `job_dictionary`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def job_dictionary(circuits, job):\n",
    "    \"\"\" Returns dictionary.\n",
    "    \n",
    "    Keyword arguments:\n",
    "    circuits -- list of quantum circuits\n",
    "    job -- execution job, should correspond to execution of keyword argument 'circuits'\n",
    "    \"\"\"\n",
    "    job_dict = {}\n",
    "    job_dict[\"circuits\"] = circuits\n",
    "    job_dict[\"job\"] = job\n",
    "    \n",
    "    return job_dict"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The reason why we do not save the execution results directly becomes apparent at a later point. In short, running a job on a real device can take time because other users' jobs can also be in the queue. By saving the job object itself, we can regularly check its state until its execution is complete.\n",
    "\n",
    "When a job has been executed, a density matrix can be reconstructed from a dictionary returned by the above function using the function `rho_from_dict` defined below. An additional input argument `cond_tomo` also allows to specify whether conditional tomography should be performed (i.e. all unsuccessful purification shots are discarded from the tomography results) or not:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def rho_from_dict(job_dict, cond_tomo):\n",
    "    \"\"\" Returns a reconstructed density matrix using quantum state tomography.\n",
    "    \n",
    "    Keyword arguments:\n",
    "    job_dict -- dictionary containing execution job and corresponding quantum circuits\n",
    "    cond_tomo -- boolean that indicates whether conditional tomography is being performed (cond_tomo=True)\n",
    "    \n",
    "    Returns:\n",
    "    rho_fit -- reconstructed density matrix\n",
    "    \"\"\"\n",
    "    if cond_tomo:\n",
    "        # Conditional quantum state tomography\n",
    "        conditional_results = purification_shots(job_dict[\"job\"].result())\n",
    "        rho_fit = StateTomographyFitter(conditional_results, job_dict[\"circuits\"]).fit()\n",
    "    else:\n",
    "        # Regular quantum state tomography\n",
    "        rho_fit = StateTomographyFitter(job_dict[\"job\"].result(), job_dict[\"circuits\"]).fit()\n",
    "    \n",
    "    return rho_fit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We summarize some of the previous functions and the execution of the tomography circuits in the function `execute_purification`, which returns exactly the dictionary containing the execution job and the corresponding tomography circuits:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def execute_purification(F_in, index, backend, shots, noise_model=None, coupling_map=None, basis_gates=None,\n",
    "                         remove_barriers=False):\n",
    "    \"\"\" Execute entire BBPSSW protocol on a target backend.\n",
    "    \n",
    "    Keyword arguments:\n",
    "    F_in            -- input fidelity\n",
    "    index           -- integer between 0 and 11, defines bilateral rotation\n",
    "    backend         -- backend to execute circuits\n",
    "    shots           -- number of shots per execution\n",
    "    noise_model     -- noise model for simulation of real device (required for noisy simulation)\n",
    "    coupling_map    -- list that specifies the qubit couplings (required for noisy simulation)\n",
    "    basis_gates     -- set of basis gates on real device (required for noisy simulation)\n",
    "    remove_barriers -- bool; if True, removes all barriers from the circuits\n",
    "    \n",
    "    Returns:\n",
    "    job_dict -- dictionary containing execution job and tomography circuits\n",
    "    \"\"\"\n",
    "    # Create experiment circuit and corresponding tomography circuits\n",
    "    experiment_circ = experiment_circuit(index)\n",
    "    tomography_circs, ancilla_tomography_circs = conditional_tomography_circuits(experiment_circ)\n",
    "    \n",
    "    # Remove barriers\n",
    "    if remove_barriers:\n",
    "        ancilla_tomography_circs[:] = [_remove_barriers(circ) for circ in ancilla_tomography_circs]\n",
    "    \n",
    "    # Execute tomography circuits on simulator (w/ noise parameters) or real device (w/o noise parameters)\n",
    "    if backend.configuration().simulator == True:\n",
    "        job = execute(ancilla_tomography_circs, backend, shots=shots,\n",
    "                      noise_model=noise_model, coupling_map=coupling_map, basis_gates=basis_gates,\n",
    "                      parameter_binds=[{theta: F_to_theta(F_in)}])\n",
    "    else:   \n",
    "        job = execute(ancilla_tomography_circs, backend, shots=shots,\n",
    "                      parameter_binds=[{theta: F_to_theta(F_in)}])\n",
    "    \n",
    "    # Construct job dictionary\n",
    "    job_dict = job_dictionary(tomography_circs, job)\n",
    "    \n",
    "    return job_dict"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the function `execute_purification` we already consider the case of a simulation with noise, which will be dealt with in the next part (**Part IV: Noise in the simulation model**). To perform such a simulation, we need additional objects (`noise_model`, `coupling_map` and `basis_gates`), which are already included in the above function. We will explain what these objects are when required in the tutorial.\n",
    "\n",
    "In addition, we have already prepared the case that the circuits will run on a real IBM Q device. If the desired backend is a real device, i.e. the property `backend.configuration().simulator` is `False`, we do not pass the parameters for a noisy simulation (`noise_model`, `coupling_map` and `basis_gates`) to the Qiskit function `execute`.\n",
    "\n",
    "With the auxiliary functions defined above, we essentially perform the following two steps to purify entanglement (whether as a unitary or noisy simulation, or on a real device):\n",
    "1. Create and execute quantum circuits using `execute_purification`.\n",
    "2. Reconstruct a density matrix using `rho_from_dict`.\n",
    "\n",
    "### 3.2.2 Simulation\n",
    "\n",
    "In Section **3.1 Single Quantum System** we have executed the BBPSSW protocol for a single system consisting of two qubit pairs in the state $\\rho_\\textrm{S}$. In this case, the result of the protocol, i.e. the state of the control qubits, depends on the randomly selected bilateral rotation.\n",
    "\n",
    "However, we are interested in the state $\\rho'_{F'}$, which is the average state of the control qubits after applying the BBPSSW protocol. To get this state $\\rho'_{F'}$ in our simulation, we calculate the average\n",
    "\n",
    "$$\n",
    "\\rho'_{F'} \\approx \\frac{1}{N} \\sum_{n=1}^{N} \\rho_n',\n",
    "$$\n",
    "\n",
    "where $\\rho_n'$ denotes the control-qubit state after the application of the BBPSSW protocol to a single system in execution run $n$ and $N$ is the total number of execution runs. (This is equivalent to applying the protocol to $N$ quantum systems and averaging across all systems.) Here, each execution run varies in the random selection of the bilateral rotation. The number of $N$ runs must therefore be large enough to make the result independent of the selected rotations. (We would get the exact state $\\rho'_{F'}$ in the limit $N \\rightarrow \\infty$.)\n",
    "\n",
    "With the functions defined above, the implementation in Qiskit for $N=100$ runs and an initial fidlity of $F = 0.75$ looks like this to approximately get the average state $\\rho'_{F'}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "N = 100\n",
    "F = 0.75\n",
    "\n",
    "# Empty list for storing the reconstructed density matrices of all runs\n",
    "fitted_states = []\n",
    "\n",
    "# Perform runs\n",
    "for ii in range(N):\n",
    "    # Create experiment circuit and corresponding tomography circuits\n",
    "    random_idx = np.random.randint(12)\n",
    "    experiment_circ = experiment_circuit(random_idx)\n",
    "    tomography_circs, ancilla_tomography_circs = conditional_tomography_circuits(experiment_circ)\n",
    "    \n",
    "    # Execute ancilla tomography circuits on simulator\n",
    "    job = execute(ancilla_tomography_circs, backend_QASM, shots=1024, parameter_binds=[{theta: F_to_theta(F)}])\n",
    "    job_results = job.result()\n",
    "    \n",
    "    # Keep counts of successful purifications\n",
    "    conditional_results = purification_shots(job_results)\n",
    "\n",
    "    # Reconstruct density matrix of control qubits\n",
    "    rho_fit = StateTomographyFitter(conditional_results, tomography_circs).fit()\n",
    "    \n",
    "    # Add reconstructed density matrix to list\n",
    "    fitted_states.append(rho_fit)\n",
    "    \n",
    "# Calculate average density matrix\n",
    "rho_average = sum(fitted_states)/N"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Also here we can depict the average state `rho_average`, e.g. by representing its real and imaginary parts:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "plot_state_city(rho_average)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we want to see how well our simulation worked and calculate the fidelity $F'$,\n",
    "\n",
    "$$\n",
    "F' = \\langle \\Phi^+ | \\rho'_{F'} | \\Phi^+ \\rangle.\n",
    "$$\n",
    "\n",
    "We can use the Qiskit function `state_fidelity` which returns the fidelity between two quantum states given as input arguments."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiskit.quantum_info import state_fidelity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The output fidelity $F'$ can thus be obtained from the state `rho_average` by"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "simulated_output_fidelity = state_fidelity(rho_average, phi_plus)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We compare the fidelity `simulated_output_fidelity` coming from our simulation with the theoretical output value of the BBPSSW protocol:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "print('Input fidelity =', F)\n",
    "print('Theoretical output fidelity = {:.4f}'.format(output_fidelity(F)))\n",
    "print('Simulated output fidelity = {:.4f}'.format(simulated_output_fidelity))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The output fidelity of the simulation is close to the theoretically expected one, which can also be seen from plotting the result:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "F_in = np.arange(0.5, 1.0, 0.01)\n",
    "plt.plot(F_in, F_in, 'k', label = \"$F' = F$\")\n",
    "plt.plot(F_in, output_fidelity(F_in), label = \"$F'(F)$\")\n",
    "plt.plot(F, simulated_output_fidelity, 'ro', label='QASM simulator')\n",
    "plt.legend()\n",
    "plt.xlabel('Input Fidelity $F$')\n",
    "plt.ylabel(\"Output Fidelity $F'$\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After seeing that our Qiskit implementation of the BBPSSW protocol works, we will perform entanglement purification for different input fidelities $F$ in the range between 0.5 and 1. We define these fidelities in a numpy array `input_fidelities`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "input_fidelities = np.arange(0.5, 1, 0.05)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We simulate the BBPSSW protocol for all fidelities using the previously defined function `execute_purification`, which further shortens the implementation. We also reconstruct the density matrix in each run with the previously defined function `rho_from_dict`, in which we set the input argument to `cond_tomo=True`, since we want to perform conditional quantum state tomography here. All calculated output fidelities are then stored in a list called `output_fidelities`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 100\n",
    "\n",
    "# Empty list for storing all output fidelities\n",
    "output_fidelities = []\n",
    "\n",
    "for jj, fidelity in enumerate(input_fidelities):\n",
    "    fitted_states = []\n",
    "    \n",
    "    for ii in range(N):\n",
    "        # Create circuits, execute circuits on simulator, reconstruct density matrix of control qubits\n",
    "        random_idx = np.random.randint(12)\n",
    "        job_dict = execute_purification(fidelity, random_idx, backend_QASM, shots=1024, remove_barriers=True)\n",
    "        rho_fit = rho_from_dict(job_dict, cond_tomo=True)\n",
    "        fitted_states.append(rho_fit)\n",
    "        \n",
    "    # Calculate average density matrix\n",
    "    rho_average = sum(fitted_states)/N\n",
    "    \n",
    "    # Calculate fidelity towards Phi^+ and add fidelity to list\n",
    "    output_fidelities.append(state_fidelity(rho_average, phi_plus))\n",
    "    \n",
    "    print(str(jj+1) + \"/\" + str(len(input_fidelities)) + \" simulation runs done.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we now compare the purified fidelities with the input fidelities, we can see a good agreement between the simulated results and the theoretical expectation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "F_in = np.arange(0.5, 1.0, 0.01)\n",
    "plt.plot(F_in, F_in, 'k', label = \"$F' = F$\")\n",
    "plt.plot(F_in, output_fidelity(F_in), label = \"$F'(F)$\")\n",
    "plt.plot(input_fidelities, output_fidelities, 'ro', label='QASM simulator')\n",
    "plt.legend()\n",
    "plt.xlabel('Input Fidelity $F$')\n",
    "plt.ylabel(\"Output Fidelity $F'$\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the simulation, in which we have not included any error sources, we see that the implementation of the BBPSSW protocol for entanglement purification works.\n",
    "\n",
    "The ultimate goal is to run the BBPSSW protocol on a real IBM Q device. But before we do this, we will include noise in our simulation model to consider imperfections in the quantum computation and see how much the effects of noise affect output fidelity."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Part IV: Including Noise in the Simulation Model\n",
    "***\n",
    "\n",
    "In this part we will simulate the quantum circuit for the BBPSSW protocol using a model that includes noise. Noise in this context means errors that occur on real devices such as single-qubit gate errors, two-qubit gate errors and single-qubit readout errors. A simple noise model can e.g. be created using error parameters directly provided through the calibration data of the real devices. This is going to be the way we construct the noise model here.\n",
    "\n",
    "In the following, we show how to choose a real IBM Q device as backend and read out the necessary data using an IBM Q account. Afterwards we explain how to generate a noise model using integrated Qiskit functions and run a noisy simulation using the QASM simulator.\n",
    "\n",
    "## 4.1 IBM Q Account\n",
    "\n",
    "To use the IBM Q Experience, we first have to import the necessary library:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiskit import IBMQ"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you have previously saved your IBM Q credentials locally, you can load your IBM Q account using the `load_account` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "IBMQ.load_account()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For more information on managing your IBM Q account, please refer to the corresponding Qiskit notebook [8].\n",
    "\n",
    "The real backends are provided by so-called providers. In the case of the authors, we only have access to one provider, namely the one that grants access to the public IBM Q devices:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "IBMQ.providers()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can access this provider using the method `get_provider` by specifying `hub='ibm-q'` as input argument:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_provider = IBMQ.get_provider(hub='ibm-q')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can use the provider instance `my_provider` to access a list of available or accessible devices using the method `backends`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "my_provider.backends()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A real IBM Q device can now be selected as backend using the method `get_backend` and for this tutorial we use the device `ibmqx2`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "backend_IBMQ = my_provider.get_backend('ibmqx2')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4.2 Creation of Noise Model in Qiskit\n",
    "\n",
    "It has already been mentioned that Qiskit provides integrated functions to construct simple noise models from device parameters, such as\n",
    "\n",
    "* single-qubit gate errors,\n",
    "* two-qubit gate errors, and\n",
    "* single-qubit readout errors.\n",
    "\n",
    "A noise model can be created by providing values for the above parameters or by reading them directly from the calibration data of a real IBM Q device. In order to use this Qiskit feature, we must first import the class `NoiseModel` from the Qiskit element Aer:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiskit.providers.aer.noise import NoiseModel"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The noise parameters of a real device listed above are stored in the backend properties, and a noise model can be directly created from the real backend using the method `from_backend`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "noise_model = NoiseModel.from_backend(backend_IBMQ)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The object `noise_model` is an instance of the `NoiseModel` class and can be used in the QASM simulator to perform noisy simulations.\n",
    "\n",
    "Since not all qubits on a real device are directly coupled, it is very likely that additional gates will have to be introduced if a circuit is run on a real device. These gates swap quantum states between qubits so that multi-qubit gates can be executed. To include such swap processes in the simulation, the function `execute` takes an additional input parameter `coupling_map`, in which the qubit couplings can be provided.\n",
    "\n",
    "We can also read the qubit couplings directly from the device itself using the property `coupling_map` of a backend configuration:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "coupling_map = backend_IBMQ.configuration().coupling_map"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In addition, a real device can only execute a finite set of quantum operations directly, and all other gates are constructed from that set. The construction of generic quantum gates from these basis gates and the generation of necessary swap operations based on the qubit coupling structure are part of the compilation process of a quantum circuit.\n",
    "\n",
    "When we run a noisy simulation, we want to run exactly the same circuits that would run on the real device. We can pass the set of basis gates to the function `execute` to ensure that the simulated quantum circuits are assembled from the same basis gates available on the real device. The noise model `noise_model` is also defined with respect to the basis gates that are available on the device `ibmqx2`. Therefore the set of basis gates can e.g. be obtained from the property `basis_gates` of a noise model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "basis_gates = noise_model.basis_gates"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the `noise_model`, the `basis_gates` and the `coupling_map` we have all information at hand that is necessary for a (simple) noisy simulation of the quantum circuit for the BBPSSW protocol. \n",
    "\n",
    "## 4.3 Simulation Including Noise\n",
    "\n",
    "To perform a noisy simulation with the QASM simulator, we can again use the function `execute` and provide three additional input parameters, which are\n",
    "* `noise_model`: An instance of the `NoiseModel` class.\n",
    "* `coupling_map`: A list that specifies the qubit couplings.\n",
    "* `basis_gates`: A set of basis gates.\n",
    "\n",
    "Our function `execute_purification` has already been defined in such a way that we simply specify the additional parameters `noise_model`, `basis_gates` and `coupling_map` as input arguments, which are then fed into the function `execute` to perform the noisy simulation from which we finally obtain the reconstructed density matrix of the control qubits $\\rho'_{F'}$ .\n",
    "\n",
    "**However, to determine the output fidelity as a function of the input fidelity taking noise into account, we need to perform an additional simulation.** The reason for this is explained below.\n",
    "\n",
    "Remember that the BBPSSW protocol works for Werner states $W_F$. The input fidelity $F$ into the BBPSSW protocol is thus the fidelity of the Werner state $W_F$, which we created after the application of the twirl circuit `twirl_circ`. In the case of a noisy simulation, this fidelity $F$ is most likely less than the fidelity of the source state $\\rho_\\textrm{S}$, which we will now call $F_\\textrm{S}$,\n",
    "\n",
    "$$\n",
    "F_\\textrm{S} = \\langle \\Psi^- | \\rho_\\textrm{S} | \\Psi^- \\rangle,\n",
    "$$\n",
    "\n",
    "due to errors in the preparation process of $W_F$. I.e. in general we expect\n",
    "\n",
    "$$\n",
    "F_\\textrm{S} \\leq F,\n",
    "$$\n",
    "\n",
    "with equality only in the case of unitary time evolution without any sources of noise. The actual input fidelity $F$ is thus the fidelity of the two-qubit state after the initialization circuit `source_circ` and the application of the twirl procedure in `twirl_circ`.\n",
    "\n",
    "In order to compare the input fidelity $F$ with the output fidelity $F'$, we therefore have to simulate the circuits `source_circ` and `twirl_circ` first to find the actual input state of the BBPSSW protocol. We must also check whether it is actually a Werner state of the form $W_F$ and then calculate the actual input fidelity $F$ from this state.\n",
    "\n",
    "### 4.3.1 Initial State and Initial Fidelity\n",
    "\n",
    "As we described above, the input fidelity $F$ to the BBPSSW protocol is the fidelity of the state $W_F$ towards the Bell state $|\\Psi^-\\rangle$ obtained after the initialization circuit and twirl operation. (We still have to show that it is actually a Werner state, what we can do after the noisy simulation.)\n",
    "\n",
    "In the following we will perform a noisy simulation of a quantum circuit consisting of the circuits `source_circ` and `twirl_circ`, and to reconstruct the state $W_F$ we will perform quantum state tomography of the control qubits.\n",
    "\n",
    "Here, too, we define some helpful functions that make the implementation easier to read and understand. First we define a function `initialization_circuit` (analogous to `experiment_circuit`, see **Sec. 3.2.1**), which takes the index of the bilateral rotation as input argument and returns a corresponding experimental circuit for the initialization process including the twirl:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def initialization_circuit(index):\n",
    "    \"\"\" Returns initialization circuit (source + twirl) corresponding to an index.\n",
    "    \n",
    "    Keyword arguments:\n",
    "    index -- integer between 0 and 11\n",
    "    \n",
    "    Returns:\n",
    "    experiment_circ -- quantum circuit\n",
    "    \"\"\"\n",
    "    if index < 0 or index > 11:\n",
    "        raise Exception(\"Index for bilateral rotation exceeds range of possible values.\")\n",
    "    \n",
    "    # Create twirl circuit from the bilateral rotations\n",
    "    twirl_circ = QuantumCircuit(qr)\n",
    "    bilateral_rotation(twirl_circ, control_qubits, index)\n",
    "    bilateral_rotation(twirl_circ, target_qubits, index)\n",
    "    \n",
    "    experiment_circ = source_circ + twirl_circ\n",
    "    \n",
    "    return experiment_circ"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All steps to execute the initialization process are summarized in the function `execute_initialization` (analogous to `execute_purification`, see **Sec. 3.2.1**). This function creates the initialization circuit via the function `initialization circuit` from above, creates the corresponding tomography circuits with the Qiskit function `state_tomography_circuits` and executes the tomography circuits on the backend, which is specified as input parameter:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def execute_initialization(F_S, index, backend, shots, noise_model=None, coupling_map=None, basis_gates=None):\n",
    "    \"\"\" Execute initialization procedure of BBPSSW protocol (simulation of source and twirl) on a target backend.\n",
    "    \n",
    "    Keyword arguments:\n",
    "    F_S -- source fidelity\n",
    "    index -- integer between 0 and 11, defines bilateral rotation\n",
    "    backend -- backend to execute circuits\n",
    "    shots -- number of shots per execution\n",
    "    noise_model -- noise model for simulation of real device\n",
    "    coupling_map -- list that specifies the qubit couplings\n",
    "    basis_gates -- set of basis gates on real device\n",
    "    \n",
    "    Returns:\n",
    "    job_dict -- dictionary containing execution job and corrresponding quantum circuits\n",
    "    \"\"\"\n",
    "    # Create experiment circuit and corresponding tomography circuits\n",
    "    experiment_circ = initialization_circuit(index)\n",
    "    tomography_circs = state_tomography_circuits(experiment_circ, control_qubits)\n",
    "\n",
    "    # Execute tomography circuits on simulator (w/ noise parameters) or real device (w/o noise parameters)\n",
    "    if backend.configuration().simulator == True:\n",
    "        job = execute(tomography_circs, backend, shots=shots,\n",
    "                      noise_model=noise_model, coupling_map=coupling_map, basis_gates=basis_gates,\n",
    "                      parameter_binds=[{theta: F_to_theta(F_S)}])\n",
    "    else:\n",
    "        job = execute(tomography_circs, backend, shots=shots,\n",
    "                      parameter_binds=[{theta: F_to_theta(F_S)}])\n",
    "    \n",
    "    # Construct job dictionary\n",
    "    job_dict = job_dictionary(tomography_circs, job)\n",
    "    \n",
    "    return job_dict"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function `execute_initialization` also returns a job dictionary containing the execution job and the corresponding tomography circuits.\n",
    "\n",
    "We define a list of source fidelities $F_\\textrm{S}$ as numpy array, for which we want to extract the resulting input fidelities $F$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "source_fidelities = np.arange(0.5, 1, 0.05)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In analogy to our simulation of the entire BBPSSW protocol in the absence of noise in **Part III** we calculate an average over $N=100$ execution runs, \n",
    "\n",
    "$$\n",
    "W_F \\approx \\frac{1}{N} \\sum_{n=1}^N \\rho_n,\n",
    "$$\n",
    "\n",
    "where $\\rho_n$ denotes the state after the random bilateral rotation in run $n$ which can be obtained by executing the function `execute_initialization`. We store the reconstructed initial states and the corresponding input fidelities in the lists `simulated_initial_state` and `input_fidelities_noise`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "N = 100\n",
    "\n",
    "# Empty lists for storing all input states and input fidelities\n",
    "simulated_initial_states = []\n",
    "input_fidelities_noise = []\n",
    "\n",
    "for jj, fidelity in enumerate(source_fidelities):\n",
    "    fitted_states = []\n",
    "    \n",
    "    for ii in range(N):\n",
    "        # Create circuits, execute circuits on simulator, reconstruct density matrix of control qubits\n",
    "        random_idx = np.random.randint(12)\n",
    "        job_dict = execute_initialization(fidelity, random_idx, backend_QASM, shots=1024,\n",
    "                                          noise_model=noise_model, coupling_map=coupling_map, basis_gates=basis_gates)\n",
    "        rho_fit = rho_from_dict(job_dict, cond_tomo=False)\n",
    "        fitted_states.append(rho_fit)\n",
    "      \n",
    "    # Calculate average density matrix and add state to list\n",
    "    rho_average = sum(fitted_states)/N\n",
    "    simulated_initial_states.append(rho_average)\n",
    "    \n",
    "    # Calculate fidelity towards Psi^- and add fidelity to list\n",
    "    input_fidelities_noise.append(state_fidelity(rho_average, psi_minus))\n",
    "    \n",
    "    print(str(jj+1) + \"/\" + str(len(source_fidelities)) + \" simulation runs done.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If one compares this code with the corresponding code from **Sec. 3.2.2**, one notices on the one hand that we perform regular quantum state tomography here, i.e. we set `cond_tomo=False` in the function `rho_from_dict`. Conditional quantum state tomography is not required here because we want to find out what the state of the control qubits after the initialization process is without restriction. (We might as well reconstruct the state of the target qubits here or take an average of both.) On the other hand, we calculate the fidelity towards the Bell state $|\\Psi^-\\rangle$ because we have not yet exchanged the $|\\Psi^-\\rangle$ and $|\\Phi^+\\rangle$ components (see step 2 of **Sec. 1.** and **Sec. 2.3.4**).\n",
    "\n",
    "We can now plot the resulting input fidelities from the list `input_fidelities_noise` against the source fidelties `source_fidelities` to see the degree of noise degradation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "F_S = np.arange(0.5, 1.0, 0.01)\n",
    "plt.plot(F_S, F_S, 'k', label = \"$F = F_\\mathrm{S}$\")\n",
    "plt.plot(source_fidelities, input_fidelities_noise, 'ro', label='QASM simulator w/ noise')\n",
    "plt.legend()\n",
    "plt.xlabel('Source Fidelity $F_\\mathrm{S}$')\n",
    "plt.ylabel(\"Input Fidelity $F$\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For all source fidelities $F_\\mathrm{S}$ we observe in the simulation results that the corresponding input fidelity $F$ is less than $F_\\mathrm{S}$.\n",
    "\n",
    "Besides the magnitude of the initial fidelity, another crucial requirement for the functionality of the BBPSSW protocol is that the input states have Werner form $W_F$. \n",
    "\n",
    "> **Exercise 3**  \n",
    "Write a function `state_fidelity_Werner` that returns the fidelity of a two-qubit state towards the Werner state $W_F$. The function should use the following input arguments:\n",
    "> - `F_Werner`: Fidelity $F$ of a Werner state $W_F$.\n",
    "> - `state`: A two-qubit quantum state."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def state_fidelity_Werner(F_Werner, state):\n",
    "    \"\"\"Returns the fidelity of a two-qubit quantum state to a Werner state.\n",
    "    \n",
    "    Keyword arguments:\n",
    "    F_Werner -- fidelity of Werner state towards Psi^-, range [0, 1]\n",
    "    state -- two-qubit quantum state\n",
    "    \"\"\"\n",
    "    if F_Werner < 0 or F_Werner > 1:\n",
    "        raise Exception('Fidelity must be between 0 and 1.')\n",
    "    \n",
    "    ###\n",
    "    ### Enter your code below\n",
    "    ###\n",
    "    \n",
    "    fidelity = state_fidelity(Werner_state(F_Werner), state)\n",
    "    \n",
    "    return fidelity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> **Exercise 4**  \n",
    "We want to find out how well a Werner state can be initialized in the noisy simulation. For all states in `simulated_initial_states`, calculate the fidelity towards the corresponding Werner state using the function `state_fidelity_Werner`, which you defined above in exercise 5. Store the calculated fidelities in a list called `initial_fidelities_Werner`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "initial_fidelities_Werner = []\n",
    "\n",
    "###\n",
    "### Enter your code below\n",
    "###\n",
    "for fidelity, state in zip(input_fidelities_noise, simulated_initial_states):\n",
    "    fidelity = state_fidelity_Werner(fidelity, state)\n",
    "    initial_fidelities_Werner.append(fidelity)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can plot the fidelities `initial_fidelities_Werner` against the source fidelities $F_\\mathrm{S}$ and see that Werner states can be initialized with fidelities of about 99 % if the calculation of `initial_fidelities_Werner` is correct:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "F_S = np.arange(0.5, 1.0, 0.01)\n",
    "plt.plot(F_S, np.ones(len(F_S)), 'k', label = \"$Fidelity = 1$\")\n",
    "plt.plot(source_fidelities, initial_fidelities_Werner, 'ro', label='QASM simulator w/ noise')\n",
    "plt.legend()\n",
    "plt.xlabel('Source Fidelity $F_\\mathrm{S}$')\n",
    "plt.ylabel('Fidelity towards Werner state')\n",
    "plt.ylim(0.95, 1.001);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In summary, we have shown in this noisy simulation that after initialization by the quantum circuits `source_circ` and `twirl_circ` a Werner state $W_F$ can be generated, but with a fidelity $F$ that is less than the fidelity $F_\\mathrm{S}$ of the original source state.\n",
    "\n",
    "### 4.3.2 Entanglement Purification in the Presence of Noise\n",
    "\n",
    "In **Sec. 3.2.1** we have already defined the function `execute_purification`. To perform a noisy simulation of the BBPSSW protocol on the QASM simulator, we can simply pass the additional parameters `noise_model` and `coupling_map` to this function, which are then finally passed to the Qiskit function `execute`.\n",
    "\n",
    "Again we calculate the average of $N=100$ runs,\n",
    "\n",
    "$$\n",
    "\\rho'_{F'} \\approx \\frac{1}{N} \\sum_{n=1}^{N} \\rho_n',\n",
    "$$\n",
    "\n",
    "where $\\rho_n'$ now denotes the state of the control qubits after the application of the BBPSSW protocol in the execution run $n$, obtained by executing the function `execute_purification`. We perform this calculation for all source fidelities defined in `source_fidelities` and save the output fidelities in the list `output_fidelities_noise`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "N = 100\n",
    "\n",
    "# Empty list for storing all output fidelities\n",
    "output_fidelities_noise = []\n",
    "\n",
    "for jj, fidelity in enumerate(source_fidelities):\n",
    "    fitted_states = []\n",
    "    \n",
    "    for ii in range(N):\n",
    "        # Create circuits, execute circuits on simulator, reconstruct density matrix of control qubits\n",
    "        random_idx = np.random.randint(12)\n",
    "        job_dict = execute_purification(fidelity, random_idx, backend_QASM, shots=1024,\n",
    "                                        noise_model=noise_model, coupling_map=coupling_map, basis_gates=basis_gates)\n",
    "        rho_fit = rho_from_dict(job_dict, cond_tomo=True)\n",
    "        fitted_states.append(rho_fit)\n",
    "        \n",
    "    # Calculate average density matrix\n",
    "    rho_average = sum(fitted_states)/N\n",
    "    \n",
    "    # Calculate fidelity towards Phi^+ and add fidelity to list\n",
    "    output_fidelities_noise.append(state_fidelity(rho_average, phi_plus))\n",
    "    \n",
    "    print(str(jj+1) + \"/\" + str(len(source_fidelities)) + \" simulation runs done.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Comparing this code with the corresponding code from **Sec. 3.2.2**, the only difference is that the arguments `noise_model` and `coupling_map` are passed to the function `execute_purification`.\n",
    "\n",
    "Now let's see how well the purification worked if we include realistic noise parameters in the simulation. To do this, the output fidelities `output_fidelities_noise` are plotted as a function of the input fidelities `input_fidelities_noise`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "F_in = np.arange(0.5, 1.0, 0.01)\n",
    "plt.plot(F_in, F_in, 'k', label = \"$F' = F$\")\n",
    "plt.plot(F_in, output_fidelity(F_in), label = \"$F'(F)$\")\n",
    "plt.plot(input_fidelities, output_fidelities, 'ro', label='QASM simulator')\n",
    "plt.plot(input_fidelities_noise, output_fidelities_noise, color='orange', marker='s', linewidth=0,\\\n",
    "         label='QASM simulator w/ noise')\n",
    "plt.legend()\n",
    "plt.xlabel('Input Fidelity $F$')\n",
    "plt.ylabel(\"Output Fidelity $F'$\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that due to the presence of noise, the output fidelity (orange) is smaller than in the case of a unitary time evolution we have examined in **Part III** (red). Nevertheless, for fidelities around 0.7 to 0.8, where the increase in fidelity is largest, it can be observed that the output fidelity is slightly increased. Therefore, it may even be possible to successfully run the entanglement purification protocol on a real IBM Q device.\n",
    "\n",
    "And this is exactly what we will do in the following, last part of the tutorial."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Part V: Entanglement Purification on a Real IBM Q Device\n",
    "***\n",
    "\n",
    "In this last part of the tutorial, we will run the BBPSSW protocol for entanglement purification on a real IBM Q device using the IBM Q Experience.\n",
    "\n",
    "The first section explains the necessary details when using the IBM Q Experience. Then we proceed as in the noisy simulation, i.e. we first determine the true input fidelity into the purification protocol and then carry out the complete BBPSSW protocol.\n",
    "\n",
    "## 5.1 IBM Q Experience\n",
    "\n",
    "We have already defined the functions `execute_initialization` and `execute_purification` in such a way that we can simply pass the real IBM Q device as the target backend to execute the tomography circuits on it. The reason for this is that the Qiskit function `execute` itself sends a job to a real IBM Q device when this is selected as the backend.\n",
    "\n",
    "Furthermore, we have already defined the real device `ibmqx2` as the backend `Backend_IBMQ`, so we are basically ready to execute our quantum circuits on the `ibmqx2` machine.\n",
    "\n",
    "When working with real IBM Q devices via the IBM Q Experience, there are two points to consider when performing the entanglement purification experiment as in the previous sections. On the one hand we have to take into account a certain waiting time until the execution of a job and on the other hand we can only send a limited number of jobs to a real device. These two points will be explained in the next two subsections.\n",
    "\n",
    "### 5.1.1 Device Queue, Job Status, and Job Monitor\n",
    "\n",
    "We have already mentioned that a real device usually has a queue of pending jobs from other users. This means that, in addition to the actual execution time of a quantum circuit, the user must also wait until all pending jobs in the queue have been executed. A job has run successfully when it could be transmitted to the real device, passed the queue, and the quantum circuit could finally be executed. Only then can we obtain the results using the method `results` and perform all subsequent steps to reconstruct a density matrix.\n",
    "\n",
    "This is exactly why we store the job object itself in a dictionary in the functions `execute_initialization` and `execute_purification`, so that we can retrieve the results at a later time when the job has been executed.\n",
    "\n",
    "The current status of a job can be determined using the `status` method of the job class. Possible statuses are e.g. that the job is initialized, queued, actively running or done.\n",
    "\n",
    "Qiskit also offers the possibility to continuously monitor a job with the function `job_monitor`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiskit.tools.monitor import job_monitor"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, note that when using the `job_monitor` function, the notebook kernel remains busy until a job has been executed (or aborted for some reason). Since a job sent to a real device may be queued for a long period of time (up to hours or even days, depending on the workload of the device), the job monitor should be used carefully.\n",
    "\n",
    "The number of pending jobs on an IBM Q device can be checked either on your IBM Q Experience - Dashboard (https://quantum-computing.ibm.com/) or with a Qiskit Jupyter tool `backend_overview`. For this we also import the Jupyter Tools:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import qiskit.tools.jupyter"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This gives us a good overview of all IBM Q devices available to us:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "%qiskit_backend_overview"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5.1.2 Personal Credits\n",
    "\n",
    "The IBM Q Experience regulates the use of real devices with the help of a credit system. This means that an execution on a real device costs a certain number of credits and after a job has been executed, the used credits are replenished. The number of credits that an execution actually costs depends, for example, on the number of experiment shots and the size of the quantum circuits to be executed.\n",
    "\n",
    "As a private user you have a limited number of credits and therefore you can only send a maximum number of jobs to a real device. \n",
    "\n",
    "Remember that we need to execute a series of $N$ runs to implement the twirl operation, which on average creates a Werner state of $W_F$,\n",
    "\n",
    "$$\n",
    "W_F \\approx \\frac{1}{N} \\sum_{n=1}^N \\rho_n,\n",
    "$$\n",
    "\n",
    "where $\\rho_n$ denotes the state after the random bilateral rotation in run $n$ which can be obtained by executing the function `execute_initialization`. And we must also execute this series of $N$ runs to determine the average output state $\\rho'_{F'}$ of the BBPSSW protocol,\n",
    "\n",
    "$$\n",
    "\\rho'_{F'} \\approx \\frac{1}{N} \\sum_{n=1}^{N} \\rho_n',\n",
    "$$\n",
    "\n",
    "where $\\rho_n'$ now denotes the state of the control qubits after the application of the BBPSSW protocol in the execution run $n$, obtained by executing the function `execute_purification`. \n",
    "\n",
    "To make this averaging process work, we have chosen $N=100$ runs. If we select 1024 shots per experiment, the execution of all tomography circuits through `execute_initialization` or `execute_purification` in one job costs 3 credits. This means that we would need 300 credits to send all $N$ runs to the real device.\n",
    "\n",
    "As a private user you usually don't have that many credits. (In our case, for example, it is 15 credits.) The entire experiments must therefore be divided into smaller batches, which are then sent separately to the device and executed there.\n",
    "\n",
    "Furthermore, we want to reduce the number of runs $N$ as much as possible in order to keep the execution time as short as possible. To reduce the number of necessary runs, we will adjust the averaging procedure to obtain a Werner state $W_F$ or the final state $\\rho'_{F'}$. Instead of chosing randomly from the set of bilateral rotations, which would be the only way to go in a quantum communication network, we will apply each element of the set once to theoretically obtain an exact Werner state,\n",
    "\n",
    "Instead of randomly choosing from the set of bilateral rotations $\\{U_i\\}$, which would be the only practical possibility in a quantum communication network, we will use each element of the set $\\{U_i\\}$ once to theoretically obtain an exact Werner state at the beginning,\n",
    "\n",
    "$$\n",
    "W_F = \\frac{1}{12} \\sum_{i=1}^{12} U_i \\rho_\\textrm{S} U_i^\\dagger,\n",
    "$$\n",
    "\n",
    "where $\\rho_\\textrm{S}$ is the source state before the application of the twirl. So now we only have to execute 12 jobs to generate a Werner state on average as input for the purification protocol.\n",
    "\n",
    "In this tutorial we divide these 12 jobs into batches of 4. We then send only one batch at a time to the device and wait until the experiments have been completed, so that there are always enough credits available to execute all jobs of a batch:\n",
    "\n",
    "| Batch #  | Bilateral Rotations $i$  | Num. of Credits  |\n",
    "|:-:|:-:|:-:|\n",
    "| 1  |  1 - 4 | 12  | \n",
    "| 2  |  5 - 8 | 12  |\n",
    "| 3  |  9 - 12 | 12  |\n",
    "\n",
    "\n",
    "When all batches have been successfully executed, we can finally follow the procedures in **Part III** and **Part IV** to reconstruct the density matrices $W_F$ and $\\rho'_{F'}$ from the execution results with our predefined function `rho_from_dict`.\n",
    "\n",
    "In case you have a different number of credits available, you will see that you can easily adapt the following implementation to your own case.\n",
    "\n",
    "## 5.2 Initial State and Initial Fidelity\n",
    "\n",
    "Due to imperfections in the real device, the actual input fidelity $F$ to the BBPSSW protocol is less than the source fidelity $F_\\textrm{S}$ (cf. the noisy simulation and the plot in **Sec. 4.3.1**). As in the case of the noisy simulation, we must therefore first determine the actual input fidelity to the BBPSSW protocol.\n",
    "\n",
    "In the previous section, we explained why we cannot run the BBPSSW protocol on a real device exactly as we did in the simulations. As a private user you have a limited number of credits and we therefore only run the BBPSSW protocol on the real device for a single source fidelity of $F_\\textrm{S}=1$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "F_source = 0.8"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To execute the quantum circuits for the initialization of a Werner state $W_F$, we can use the function `execute_initialization` again. Each execution of this function returns a job dictionary and we store the dictionary of each run $i$ in a list `initialization_jobs`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "initialization_jobs = []"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use the list elements later to retrieve the job results when the execution of the corresponding quantum circuits is complete.\n",
    "\n",
    "So let's start with the execution of the first batch on the IBM Q device `ibmqx2`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for ii in range(0, 4):\n",
    "    job_dict = execute_initialization(F_source, ii, backend_IBMQ, shots=1024)\n",
    "    initialization_jobs.append(job_dict)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remember: Each element of this batch now consists of the tomography circuits belonging to the quantum circuit which first generates the source state $\\rho_\\textrm{S}$ and then implements the bilateral rotation $U_i$. And the whole thing is done for the indices $i \\in \\{1, 2, 3, 4\\}$.\n",
    "\n",
    "After running the above cell, we can check the status of each of the four jobs using the method `status`. E.g. for the first entry in the list `initialization_jobs`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "initialization_jobs[0][\"job\"].status()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have also mentioned another method to continuously monitor the status of a job by using the function `job_monitor`. We will do this for the latest job, i.e. the last entry in `initialization_jobs`, but remember that the kernel will be busy until the job is finished (or an error has occurred):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "job_monitor(initialization_jobs[-1][\"job\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the job was successfully executed, i.e. no error occurred and the job was not canceled, the output of the above cell should print `Job Status: job has successfully run`.\n",
    "\n",
    "If this is the case, we can the start the execution of the second batch by simpy changing the indices of the bilateral rotations to $i \\in \\{5,6,7,8\\}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for ii in range(4, 8):\n",
    "    job_dict = execute_initialization(F_source, ii, backend_IBMQ, shots=1024)\n",
    "    initialization_jobs.append(job_dict)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, too, we use the job monitor to watch the status of the latest job:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "job_monitor(initialization_jobs[-1][\"job\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If this job has also run successfully, we can execute the third and last batch for the indices $i \\in \\{9,10,11,12\\}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for ii in range(8, 12):\n",
    "    job_dict = execute_initialization(F_source, ii, backend_IBMQ, shots=1024)\n",
    "    initialization_jobs.append(job_dict)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We monitor the status of the job that was last sent again:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "job_monitor(initialization_jobs[-1][\"job\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To make sure that all jobs have actually run successfully, you can, for example, check the status of all jobs that we have now executed on the real IBM Q device:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "for i, entry in enumerate(initialization_jobs):\n",
    "    print('i = ' + str(i+1) + ': ' + str(entry[\"job\"].status()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the status of all jobs is `DONE`, we can now reconstruct the quantum state of the control qubits.\n",
    "\n",
    "In the noisy simulation (**Sec. 4.3.1**) we did this for each execution run with the previously defined function `rho_from_dict`, to which we passed the job dictionary of this run and set the value of the input parameter `cond_tomo=False`. (We perform regular quantum state tomography after initialization.)\n",
    "\n",
    "Here, in the case of a real device, we can use exactly the same procedure to reconstruct the density matrices of the control qubits using the function `rho_from_dict`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "experimental_initial_states = []\n",
    "\n",
    "for entry in initialization_jobs:\n",
    "    rho_fit = rho_from_dict(entry, cond_tomo=False)\n",
    "    experimental_initial_states.append(rho_fit)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We store the reconstructed initial state of each run $i$ in the list `experimental_initial_states` and then calculate the average:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rho_average = sum(experimental_initial_states)/12"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The state `rho_average` should now be close to a Werner state $W_F$ to be useful as an input state for the BBPSSW protocol. To check this, we must first calculate the fidelity $F$ towards the Bell state $|\\Psi^-\\rangle$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "experimental_input_fidelity = state_fidelity(rho_average, psi_minus)\n",
    "\n",
    "print('Source fidelity =', F_source)\n",
    "print('Experimental input fidelity = {:.4f}'.format(experimental_input_fidelity))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This fidelity now corresponds to the actual initial fidelity, if a real IBM Q device is used as backend and is therefore called `experimental_input_fidelity`.\n",
    "\n",
    "(For comparison: In experiments to test the tutorial we got fidelities `experimental_input_fidelity` of about 75 % using the IBM Q device `ibmqx2`.)\n",
    "\n",
    "As we could expect from the noisy simulation, the fidelity `experimental_input_fidelity` is smaller than the source fidelity `F_source`, i.e.:\n",
    "\n",
    "$$\n",
    "F < F_\\textrm{S}.\n",
    "$$\n",
    "\n",
    "We can now calculate the fidelity of the state `rho_average` towards the Werner state $W_F$ (where $F$ is now the `experimental_input_fidelity`):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fidelity_Werner = state_fidelity_Werner(experimental_input_fidelity, rho_average)\n",
    "\n",
    "print('Fidelity towards Werner state = {:.4f}'.format(fidelity_Werner))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(For comparison: In experiments to test the tutorial we got fidelities `fidelity_Werner` of about 95 % using the IBM Q device `ibmqx2`.)\n",
    "\n",
    "The results we obtained from the real IBM Q device `ibmqx2` can now be compared with the results of the noisy simulation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "F_S = np.arange(0.5, 1.0, 0.01)\n",
    "plt.plot(F_S, F_S, 'k', label = \"$F = F_\\mathrm{S}$\")\n",
    "plt.plot(source_fidelities, input_fidelities_noise, 'ro', label='QASM simulator w/ noise')\n",
    "plt.plot(F_source, experimental_input_fidelity, 'gx', label=backend_IBMQ.name())\n",
    "plt.legend()\n",
    "plt.xlabel('Source Fidelity $F_\\mathrm{S}$')\n",
    "plt.ylabel(\"Input Fidelity $F$\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Probably the experimentally observed input fidelity on the device `ibmqx2` is smaller than that of the noisy simulation. There may be various reasons why the results of the simulation and the use of the real device differ. For example, we used a very simple noise model for the simulation, which most likely does not take into account all the effects that lead to errors in the computation on the real device.\n",
    "\n",
    "An analysis of how good the simulation is should not be part of this tutorial. The focus of this section is basically to execute the quantum circuits correctly on a real IBM Q device using the IBM Q Experience.\n",
    "\n",
    "Now that we have obtained an experimental value for the input fidelity $F$, in the last step of this tutorial we will run the complete BBPSSW protocol on the `ibmqx2` device.\n",
    "\n",
    "## 5.3 Execution of the BBPSSW Protocol on a Real Device\n",
    "\n",
    "The execution of the complete purification experiment is similar to the steps in the previous section regarding the splitting of all necessary jobs into smaller batches. This is again due to the limited number of credits.\n",
    "\n",
    "To run the BBPSSW protocol on a real IBM Q device, we will use our function `execute_purification` and simply pass the IBM Q device `ibmqx2` as the target backend.\n",
    "\n",
    "We use each element of the set $\\{U_i\\}$ once again to implement an exact twirl. The output state of the control qubits $\\rho'_{F'}$ after to the BBPSSW protocol will then be the following average\n",
    "\n",
    "$$\n",
    "\\rho'_{F'} = \\frac{1}{12} \\sum_{i=1}^{N} \\rho_i',\n",
    "$$\n",
    "\n",
    "where $\\rho_i'$ is the control qubit state if the bilateral rotation $U_i$ was applied during the initialization process.\n",
    "\n",
    "This time we store the jobs of each run $i$ in the list `purification_jobs`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "purification_jobs = []"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We start with the execution of the first batch on the IBM Q device `ibmqx2`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for ii in range(0, 4):\n",
    "    job_dict = execute_purification(F_source, ii, backend_IBMQ, shots=1024)\n",
    "    purification_jobs.append(job_dict)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remember once again: Each element of this batch now consists of the tomography circuits to the quantum circuit, which uses the bilateral rotation $U_i$ in the twirl and then implements the entire BBPSSW protocol. And the whole thing is done for the indices $i \\in \\{1, 2, 3, 4\\}$.\n",
    "\n",
    "Here, too, we use the job monitor to watch the status of the latest job:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "job_monitor(purification_jobs[-1][\"job\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We unfortunately obtain an error here. These circuits cannot be executed because the barriers are still in our way (remember the explanation at the end of **Part II**). But we already took care of the by passing the argument `remove_barriers=True` to the function `execute_purification`. So let's start over again:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "purification_jobs = []"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for ii in range(0, 4):\n",
    "    job_dict = execute_purification(F_source, ii, backend_IBMQ, shots=1024, remove_barriers=True)\n",
    "    purification_jobs.append(job_dict)\n",
    "    \n",
    "job_monitor(purification_jobs[-1][\"job\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the job has successfully run, we can the start the execution of the second batch by simpy changing the indices of the bilateral rotations to $i \\in \\{5,6,7,8\\}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for ii in range(4, 8):\n",
    "    job_dict = execute_purification(F_source, ii, backend_IBMQ, shots=1024, remove_barriers=True)\n",
    "    purification_jobs.append(job_dict)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "job_monitor(purification_jobs[-1][\"job\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And if the latest job of the second batch has successfully run, the third batch can be sent to the device using the indices $i \\in \\{9,10,11,12\\}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for ii in range(8, 12):\n",
    "    job_dict = execute_purification(F_source, ii, backend_IBMQ, shots=1024, remove_barriers=True)\n",
    "    purification_jobs.append(job_dict)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "job_monitor(purification_jobs[-1][\"job\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To ensure that all jobs have been successfully executed, we will again check the status of all jobs that we have just executed on the real IBM Q device:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "for i, entry in enumerate(purification_jobs):\n",
    "    print('i = ' + str(i+1) + ': ' + str(entry[\"job\"].status()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When the status of all jobs is `DONE`, we can continue to reconstruct the state of the control qubits using conditional quantum state tomography. Again we follow the steps of the noisy simulation (**Sec. 4.3.2**) and do this  for each execution run with the previously defined function `rho_from_dict`, to which we pass the job dictionary of this run and set the value of the input parameter `cond_tomo=True`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "experimental_output_states = []\n",
    "\n",
    "for entry in purification_jobs:\n",
    "    rho_fit = rho_from_dict(entry, cond_tomo=True)\n",
    "    experimental_output_states.append(rho_fit)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we store the reconstructed output state of each run $i$ in the list `experimental_output_states`. We then calculate the average of all reconstructed states to obtain $\\rho'_{F'}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rho_average = sum(experimental_output_states)/12"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have the state $\\rho'_{F'}$, we are finally able to calculate the output fidelity $F'$ after executing the BBPSSW protocol on the real device,\n",
    "\n",
    "$$\n",
    "F' = \\langle \\Phi^+ | \\rho'_{F'}| \\Phi^+ \\rangle,\n",
    "$$\n",
    "\n",
    "and compare it to the experimental input fidelity $F$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "experimental_output_fidelity = state_fidelity(rho_average, phi_plus)\n",
    "\n",
    "print('Source fidelity =', F_source)\n",
    "print('Experimental input fidelity = {:.4f}'.format(experimental_input_fidelity))\n",
    "print('Theoretical output fidelity = {:.4f}'.format(output_fidelity(experimental_input_fidelity)))\n",
    "print('Experimental output fidelity = {:.4f}'.format(experimental_output_fidelity))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(For comparison: In experiments to test the tutorial we got output fidelities `experimental_output_fidelity` of about 70 - 75 % using the IBM Q device `ibmqx2`.)\n",
    "\n",
    "If the current errors in the system are not too large, the experimental output fidelity may be greater than the input fidelity. However, we will not be able to make any concrete statements about the result. This would require a much more thorough analysis, which cannot be part of this tutorial. For example, to compare input and output fidelity, two quantities observed in two different experiments, it would be necessary to prove that the experimental conditions were the same in both cases. Also a fundamental, detailed error analysis cannot be omitted.\n",
    "\n",
    "Nevertheless, at the end we summarize all the results we have received in this tutorial:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "F_in = np.arange(0.5, 1.0, 0.01)\n",
    "plt.plot(F_in, F_in, 'k', label = \"$F' = F$\")\n",
    "plt.plot(F_in, output_fidelity(F_in), label = \"$F'(F)$\")\n",
    "plt.plot(input_fidelities, output_fidelities, 'ro', label='QASM simulator')\n",
    "plt.plot(input_fidelities_noise, output_fidelities_noise, color='orange', marker='s', linewidth=0,\\\n",
    "         label='QASM simulator w/ noise')\n",
    "plt.plot(experimental_input_fidelity, experimental_output_fidelity, 'gx', label=backend_IBMQ.name())\n",
    "plt.legend()\n",
    "plt.xlabel('Input Fidelity $F$')\n",
    "plt.ylabel(\"Output Fidelity $F'$\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When we compare the results of the `ibmqx2` machine and the noisy simulation, we see a good agreement for the value of the output fidelity $F'$. However, we have already seen that the simulator significantly overestimates the value for the actual initial fidelity $F$. The errors that come from the part of the circuit that belongs exclusively to the BBPSSW protocol, i.e. gate errors in the circuit `purification_circ` and subsequent measurement errors, are therefore probably underestimated by the noise model used here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Part VI: Conclusion and Further Ideas\n",
    "\n",
    "In this tutorial we showed how to use the IBM framework Qiskit and the online platform IBM Q Experience for quantum information processing, using the BBPSSW protocol for entanglement purification as an example.\n",
    "\n",
    "The quantum circuit for the BBPSSW protocol consists of common single-qubit and two-qubit gates and we have explained in detail how this circuit is constructed in Qiskit. Furthermore, we have demonstrated how to use a library integrated in Qiskit to perform conditional quantum state tomography, which is also required in the BBPSSW protocol.\n",
    "\n",
    "Quantum circuits in Qiskit can be executed either on a simulator or on real IBM Q devices and we have shown how to use these two backends respectively. When using the simulator, we also explained how simple noise models can be created and integrated into the simulation.\n",
    "\n",
    "There are (at least) two possible topics that could usefully extend this tutorial. On the one hand, there are possibilities to improve the execution on real devices, and on the other hand, there are other protocols for entanglement purification in addition to the BBPSSW protocol (see e.g. [9] and the references therein).\n",
    "\n",
    "Some possibilities that could improve the quality of execution on a real device are:\n",
    "\n",
    "* **Error mitigation through calibration circuits**  \n",
    "The Qiskit element Ignis provides noise mitigation functions through specific calibration routines. The use of these functions could therefore lead to a better performance of the execution on the real device and thus to a more significant increase in the output fidelity $F'$.\n",
    "\n",
    "\n",
    "* **Circuit optimization with internal compiler**  \n",
    "A quantum circuit must be mapped to the physical layout of a real device, i.e. the coupling of the individual qubits with each other and the available set quantum gates is taken into account. This is done via a compiler and the result is a compiled quantum circuit which is actually executed on the real device. The Qiskit compiler knows four optimization levels and you can use this possibility to reduce the length of the compiled circuit, i.e. the number of operations required, and thus suppress the effects of errors.\n",
    "\n",
    "\n",
    "* **Circuit optimization with external compiler**  \n",
    "Another possibility is the use of external compilers, such as t|ket⟩™, to optimize quantum circuits. In addition, one could compare the results of different compilers with each other.\n",
    "\n",
    "As already mentioned, the BBPSSW protocol is one of many for entanglement purification. Some of these protocols work very similarly and differ only in the bilateral two-qubit interaction that Alice and Bob apply locally to their qubits. It is straightforward to implement such a protocol by simply changing the gates that make up the circuit `purification_circ`.\n",
    "\n",
    "But there are also protocols that work with three input pairs $\\rho_\\textrm{S}$, for example, and require a specific measurement outcome of two of these pairs for successful purification. In order to implement such a protocol in Qiskit according to the procedure presented in this tutorial, one would also have to adapt the creation of the ancilla tomography circuits and change the function `purification_shots` so that the conditions for a successful purification are fulfilled.\n",
    "\n",
    "Another provider has just been released that enables access to the trapped ion device of Alpine Quantum Technologies (AQT) [10,11]. If you have AQT credentials, it would be great to see how well the purification protocol works on the trapped ion device. The results could then be compared with earlier experimental results demonstrating entanglement purification for ion qubits [12].\n",
    "\n",
    "Finally, we would like to emphasize that by providing these ideas, we hope to encourage other users to continue learning with this tutorial notebook by trying out the implementation of the ideas."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# References\n",
    "***\n",
    "\n",
    "[1] C. H. Bennett, H. J. Bernstein, S. Popescu, and B. Schumacher, *Concentrating partial entanglement by local operations*, Phys. Rev. A **53**, 2046 (1996).\n",
    "\n",
    "[2] C. H. Bennett, G. Brassard, C. Crépeau, R. Jozsa, A. Peres, and W. K. Wootters, *Teleporting an unknown quantum state via dual classical and Einstein-Podolsky-Rosen channels*, Phys. Rev. Lett. **70**, 1895 (1993).\n",
    "\n",
    "[3] A. K. Ekert, *Quantum cryptography based on Bell’s theorem*, Phys. Rev. Lett. **67**, 661 (1991).\n",
    "\n",
    "[4] H.-J. Briegel, W. Dür, J. I. Cirac, and P. Zoller, *Quantum Repeaters: The Role of Imperfect Local Operations in Quantum Communication*, Phys. Rev. Lett. **81**, 5932 (1998).\n",
    "\n",
    "[5] C. H. Bennett, D. P. DiVincenzo, J. A. Smolin, and W. K. Wootters, *Mixed-state entanglement\n",
    "and quantum error correction*, Phys. Rev. A **54**, 3824 (1996).\n",
    "\n",
    "[6] *Quantum State Tomography*, Retrieved from https://github.com/Qiskit/qiskit-iqx-tutorials/blob/master/qiskit/advanced/ignis/6a_state_tomography.ipynb (2019, December 1).\n",
    "\n",
    "[7] R. LaRose, A. Mari, P. J. Karalekas, N. Shammah, and W. J. Zeng, *Mitiq: A software package for error mitigation on noisy quantum computers*, arXiv:2009.04417 (2020)\n",
    "\n",
    "[8] *The IBM Q Account*, Retrieved from https://github.com/Qiskit/qiskit-iqx-tutorials/blob/5c18c22465d19b22ca19a915518f051de78ba4bc/qiskit/fundamentals/3_the_ibmq_account.ipynb (2019, December 1).\n",
    "\n",
    "[9] W. Dür and H. J. Briegel, *Entanglement purification and quantum error correction*, Rep. Prog. Phys. **70**, 1381 (2007).\n",
    "\n",
    "[10] A. Javadi-Abhari, P. Nation, and J. Gambetta, *Qiskit – Write once, target multiple architectures*, Retrieved from https://www.ibm.com/blogs/research/2019/11/qiskit-for-multiple-architectures/ (2019, December 1)\n",
    "\n",
    "[11] Retrived from https://www.aqt.eu/ (2019, December 1).\n",
    "\n",
    "[12] R. Reichle, D. Leibfried, E. Knill, J. Britton, R. B. Blakestad, J. D. Jost, C. Langer, R. Ozeri, S. Seidelin, and D. J. Wineland, *Experimental purification of two-atom entanglement*, Nature **443**, 838 (2006)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%qiskit_version_table"
   ]
  },
  {
   "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.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}