examples/advanced_circuits_algorithms/Randomness/Randomness_Generation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Robust randomness generation on quantum processing units"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-29T21:32:38.617278Z",
     "start_time": "2023-08-29T21:32:38.389075Z"
    }
   },
   "outputs": [],
   "source": [
    "# Use Braket SDK Cost Tracking to estimate the cost to run this example\n",
    "from braket.tracking import Tracker\n",
    "\n",
    "t = Tracker().start()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Random numbers are a ubiquitous resource in computation and cryptography. For example, in security, random numbers are crucial to creating keys for encryption. Quantum random number generators (QRNGs), that make use of the inherent unpredictability in quantum physics, promise enhanced security compared to standard cryptographic pseudo-random number generators (CPRNGs) based on classical technologies.\n",
    "\n",
    "In this notebook, we implement our own QRNG. Namely, we program two separate quantum processor units (QPUs) from different suppliers in Amazon Braket to supply two streams of weakly random bits. We then show how to generate physically secure randomness from these two weak sources by means of classical post-processing based on randomness extractors. The prerequisites for the tutorial are a basic understanding of quantum states, quantum measurements, and quantum channels. For a detailed explanation of these concepts we refer to the Amazon Braket notebook [Simulating noise on Amazon Braket](https://github.com/amazon-braket/amazon-braket-examples/blob/main/examples/braket_features/Simulating_Noise_On_Amazon_Braket.ipynb).\n",
    "\n",
    "We believe that randomness generation is a practical application of nowadays available noisy and intermediate scale quantum (NISQ) technologies.\n",
    "\n",
    "## Table of contents\n",
    "\n",
    "* [Amazon Braket](#amazon_braket)\n",
    "* [Quantum circuit for randomness generation](#quantum_circuit)\n",
    "* [Quick implementation](#quick_implementation)\n",
    "* [Critical assessment](#critical_assessment)\n",
    "* [Interlude randomness extractors](#interlude)\n",
    "    * [Quantum information](#quantum_information)\n",
    "* [Extractor construction](#extractor_construction)\n",
    "    * [Example](#example)\n",
    "    * [Implementation](#implementation)\n",
    "* [Unpredictability of physical sources](#physical_sources)\n",
    "    * [Noise as leakage](#noise_leakage)\n",
    "    * [Numerical evaluation](#numerical_evaluation)\n",
    "* [Putting everything together](#putting-together)\n",
    "* [Beyond current implementation](#beyond_current)\n",
    "* [Literature](#literature)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Amazon Braket <a name=\"amazon_braket\"></a>\n",
    "\n",
    "We start out with some general Amazon Braket imports, as well as some mathemtical tools needed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-29T21:34:00.351045Z",
     "start_time": "2023-08-29T21:34:00.015307Z"
    }
   },
   "outputs": [],
   "source": [
    "# AWS imports: Import Braket SDK modules\n",
    "import math\n",
    "\n",
    "# import convex solver\n",
    "import cvxpy as cp  # general math imports\n",
    "import numpy as np\n",
    "from scipy.fft import fft, ifft\n",
    "\n",
    "from braket.aws import AwsDevice, AwsQuantumTask\n",
    "from braket.circuits import Circuit\n",
    "from braket.devices import Devices, LocalSimulator\n",
    "\n",
    "# set up local simulator device\n",
    "device = LocalSimulator()\n",
    "\n",
    "# magic word for producing visualizations in notebook\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-29T21:34:06.520777Z",
     "start_time": "2023-08-29T21:34:03.869787Z"
    }
   },
   "outputs": [],
   "source": [
    "# set up Rigetti quantum device\n",
    "rigetti = AwsDevice(Devices.Rigetti.Ankaa3)\n",
    "\n",
    "# set up IonQ quantum device\n",
    "ionq = AwsDevice(Devices.IonQ.Aria1)\n",
    "\n",
    "# set up IQM quantum device\n",
    "iqm = AwsDevice(Devices.IQM.Garnet)\n",
    "\n",
    "# simulator alternative: set up the on-demand simulator SV1\n",
    "simulator = AwsDevice(Devices.Amazon.SV1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quantum circuit for randomness generation <a name=\"quantum_circuit\"></a>\n",
    "\n",
    "Arguably the simplest way of generating a random bit on a quantum computer is as follows:\n",
    "* Prepare the basis state vector $|0\\rangle$\n",
    "* Apply the Hadamard gate $H=\\frac{1}{\\sqrt{2}}\\begin{pmatrix} 1 & 1\\\\ 1 & 1 \\end{pmatrix}$ leading to the state vector $|H\\rangle=\\frac{1}{\\sqrt{2}}\\big(|0\\rangle+|1\\rangle\\big)$\n",
    "* Measure in the computational basis $\\big\\{|0\\rangle\\langle0|,|1\\rangle\\langle1|\\big\\}$.\n",
    "\n",
    "By the laws of quantum physics, the post-measurement probability distribution is then the uniformly distributed $(1/2,1/2)$ and leads to one random bit $0/1$.\n",
    "\n",
    "In the following, we discuss how above protocol is conceptually different from randomness obtained from classical sources and show in detail how it is implemented reliable even when the underlying quantum processing units employed are noisy. By the end of this tutorial, you will be able to create your own random bits from the quantum processing units available on Amazon Braket.\n",
    "\n",
    "## Quick implementation <a name=\"quick_implementation\"></a>\n",
    "\n",
    "The Hadmard gate based quantum circuit for generating one random bit can be repeated or run in parallel $n$ times, leading to a random bit string of length $n$. The corresponding circuit is easily implemented in Amazon Braket:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-29T21:34:08.544420Z",
     "start_time": "2023-08-29T21:34:08.523880Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "T  : │  0  │\n",
      "      ┌───┐ \n",
      "q0 : ─┤ H ├─\n",
      "      └───┘ \n",
      "      ┌───┐ \n",
      "q1 : ─┤ H ├─\n",
      "      └───┘ \n",
      "      ┌───┐ \n",
      "q2 : ─┤ H ├─\n",
      "      └───┘ \n",
      "      ┌───┐ \n",
      "q3 : ─┤ H ├─\n",
      "      └───┘ \n",
      "      ┌───┐ \n",
      "q4 : ─┤ H ├─\n",
      "      └───┘ \n",
      "T  : │  0  │\n"
     ]
    }
   ],
   "source": [
    "# function for Hadamard cirquit\n",
    "def hadamard_circuit(n_qubits):\n",
    "    \"\"\"\n",
    "    function to apply Hadamard gate on each qubit\n",
    "    input: number of qubits\n",
    "    \"\"\"\n",
    "\n",
    "    # instantiate circuit object\n",
    "    circuit = Circuit()\n",
    "\n",
    "    # apply series of Hadamard gates\n",
    "    for i in range(n_qubits):\n",
    "        circuit.h(i)\n",
    "\n",
    "    return circuit\n",
    "\n",
    "\n",
    "# define circuit\n",
    "n_qubits = 5\n",
    "state = hadamard_circuit(n_qubits)\n",
    "\n",
    "# print circuit\n",
    "print(state)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us run this Hadamard circuit with $n=5$ qubits in the local quantum simulator for $m=1$ shots:\n",
    "\n",
    "(Note: We will work on actual QPUs towards the end of this tutorial.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-29T21:34:10.235692Z",
     "start_time": "2023-08-29T21:34:10.052502Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The output bit string is:  ['10010']\n"
     ]
    }
   ],
   "source": [
    "# run circuit\n",
    "m_shots = 1\n",
    "result = device.run(state, shots=m_shots).result()\n",
    "\n",
    "# get measurement shots\n",
    "counts = result.measurement_counts.keys()\n",
    "\n",
    "# print counts\n",
    "list_one = list(counts)[0]\n",
    "array_one = np.array([list_one])\n",
    "print(\"The output bit string is: \", array_one)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Critical assessment <a name=\"critical_assessment\"></a>\n",
    "\n",
    "The advantage of such quantum random number generators over implementations based on classical technologies is that the outcomes are intrinsically random. That is, according to the laws of quantum physics, the outcome of the measurement is not only hard to predict, but rather impossible to know before the measurement has taken place.\n",
    "\n",
    "However, since current quantum processing units are noisy to a certain degree, there are at least three potential problems that need to be addressed:\n",
    "* First, the noise acting on all states and operations performed might lead to a systematic bias towards the probability of getting the measurement outcomes $0$ or $1$, respectively.\n",
    "* Second, even if aforementioned noise is not biased towards certain measurement outcomes, the generated randomness is no longer solely based on intrinsically random quantum effects, but rather partly on the noise present.\n",
    "* Third, whereas by the laws of quantum physics a pure quantum state cannot be correlated to the outside world, any noisy acting on the system corresponds to information leaking to the environment. This is because no information is destroyed in quantum physics and hence, a malicious third party knowing about the noise occurring will be able to guess the generated bits (at least up to a certain degree).\n",
    "\n",
    "When the noise model acting on the quantum processor units is characterized to some degree (e.g., by means of previous benchmarking), these shortcomings can be overcome by employing two independent quantum processor units, together with an appropriate classical post-processing. The latter is based on classical algorithms from the theory of pseudo-randomness, so-called two-source extractors. This is what we discuss next.\n",
    "\n",
    "In the following, we refer a few times to the theory paper [1] that features formal cryptographic security definitions, together with mathematical proofs, as well as some statistical methods tailored to intermediate scale quantum devices. These pointers can be safely ignored when only interested in the implementation of our QRNG.\n",
    "\n",
    "## Interlude randomness extractors <a name=\"interlude\"></a>\n",
    "\n",
    "Two-source extractors allow distillation of physically secure random bits from two independent weak sources of randomness whenever they are sufficiently unpredictable to start with. The relevant measure of unpredictability is thereby given by the min-entropy of the respective sources, defined for the probability distribution $\\{p_x\\}_{x\\in X}$ as\n",
    "\n",
    "$$ \\text{$H_{\\min}(X)=-\\log_2p_{\\text{guess}(X)}$ with $p_{\\text{guess}(X)}=\\max_{x\\in X}p_x$.} $$\n",
    "\n",
    "That is, the min-entropy exactly quantifies how well we can guess the value of the source, or in other words, how unpredictable the source is. For example, for $n$-bit distributions $X$ we have $H_{\\min}(X)\\in[0,n]$, where $0$ corresponds to a deterministic distribution containing no randomness and $n$ to the perfectly random, uniform distribution.\n",
    "\n",
    "A two-source extractor is a function $\\text{Ext}:\\{0,1\\}^{n_1}\\times\\{0,1\\}^{n_2}\\to\\{0,1\\}^m$ such that for any two independent sources with min-entropy at least $H_{\\min}(X_1)\\geq k_1$ and $H_{\\min}(X_2)\\geq k_2$, respectively, the output of length $m$ is $\\epsilon\\in[0,1]$ close in variational distance to the perfectly random, uniform distribution $U_M$ of size $m$:\n",
    "\n",
    "$$ \\frac{1}{2}\\left\\|\\text{Ext}(X_1,X_2)-U_M\\right\\|_1\\leq\\epsilon. $$\n",
    "\n",
    "So, two inpedendent sources that are only weakly random get condensed by these algorithms to one output that is (nearly) perfectly random! Importantly, the output becomes truly physically random with no computational assumptions introduced.\n",
    "\n",
    "### Quantum information <a name=\"quantum_information\"></a>\n",
    "\n",
    "For our setting we need an extension of this concept, as a potentially malicious third party can collect quantum information $Q$ about the weak source of randomness $X$. The corresponding conditional min-entropy is defined as\n",
    "\n",
    "$$ H_{\\min}(X|Q)=-\\log_2p_{\\text{guess}(X|Q)}, $$\n",
    "\n",
    "where $p_{\\text{guess}(X|Q)}$ denotes the maximal probability allowed by quantum physics to guess the classical value $X$ by applying any measurements on the quantum information $Q$. We notice that even though $p_{\\text{guess}(X|Q)}$ does not have a closed form expression, it is efficiently computed by means of a semidefinite program (see the theory notes [1] for details). This is also how we will evaluate the conditional min-entropy quantity later on.\n",
    "\n",
    "Accordingly, a quantum-proof two-source extractor is then function $\\text{Ext}:\\{0,1\\}^{n_1}\\times\\{0,1\\}^{n_2}\\to\\{0,1\\}^m$ such that for any two independent sources with quantum conditional min-entropy at least $H_{\\min}(X_1|Q_1)\\geq k_1$ and $H_{\\min}(X_2|Q_2)\\geq k_2$, respectively, we have for $\\epsilon\\in[0,1]$ in quantum variational distance\n",
    "\n",
    "$$ \\frac{1}{2}\\left\\|\\rho_{\\text{Ext}(X_1,X_2)Q_1Q_2}-\\tau_{M}\\otimes\\rho_{Q_1}\\otimes\\rho_{Q_2}\\right\\|_1\\leq\\epsilon,$$\n",
    "\n",
    "where $\\tau_M$ denotes the fully mixed $m$ qubit state. That is, the extractor should not only make the $m$ output bits perfectly random, but also decouple them from any outside correlations initially present - up to the security parameter $\\epsilon\\in[0,1]$.\n",
    "\n",
    "For more details about these concepts, we refer to the theory notes [1] and references therein. All that is important to us here, is that there exist quantum-proof two-source extractors with good parameters. Next, we discuss one particular such construction that we subsequently implement in an efficient manner.\n",
    "\n",
    "## Extractor construction <a name=\"extractor_construction\"></a>\n",
    "\n",
    "In this paragraph, we provide an explicit construction of a quantum-proof two-source extractor that efficiently provides non-zero output $M$ for a wide range of sizes of the inputs $X_1$ and $X_2$. Namely, we employ a Toeplitz matrices based construction originally discussed in [2]:\n",
    "\n",
    "For the security parameter $\\epsilon\\in(0,1]$ and inputs $X_1,X_2$ of size $n$ and $n-1$, respectively, the function $\\text{Ext}:\\{0,1\\}^n\\times\\{0,1\\}^{n-1}\\to\\{0,1\\}^m$ defined below is a quantum-proof two-source randomness extractor with output size\n",
    "\n",
    "$$ m=\\left\\lfloor(k_1+k_2-n)+1-2\\log\\left(1/\\epsilon\\right)\\right\\rfloor. $$\n",
    "\n",
    "The function is explicitly given via the vector-matrix multiplication\n",
    "\n",
    "$$ (x,y)\\mapsto \\text{Ext}(x,y)=x\\cdot(T(y)|1_m)^T \\mod{2}$$\n",
    "\n",
    "featuring the Toeplitz matrix\n",
    "\n",
    "$$ T(y)=\\begin{pmatrix}\n",
    "y_0 & y_1 & \\ldots & y_{n-m-1}\\\\\n",
    "y_{-1} & y_0 & \\ldots & y_{n-m-2}\\\\\n",
    "\\vdots & \\vdots & \\vdots & \\vdots\\\\\n",
    "y_{1-m} & y_{2-m} & \\ldots & y_{n-2m}\n",
    "\\end{pmatrix}\\;\\text{from}\\;y=(y_{1-m},\\ldots,y_0,\\ldots,y_{n-m-1})\\in\\{0,1\\}^{n-1}, $$\n",
    "\n",
    "The quantum-proof property of this construction, as well as its complexity, is explicitly disccused in the theory notes [1].\n",
    "\n",
    "For our setting, we will have sources with linear min-entropy rates $k_i=\\alpha_i\\cdot n$ for $i=1,2$. The Toeplitz construction then works whenever $\\alpha_1+\\alpha_2-1>0$ and we can compute the required input size for fixed output size as\n",
    "\n",
    "$$ n=\\left\\lfloor\\frac{m-1+2\\log(1/\\epsilon)}{\\alpha_1+\\alpha_2-1}\\right\\rfloor.$$\n",
    "\n",
    "A simple example shows that these numbers work well in practice, even for very small input sizes around $n\\approx100$.\n",
    "\n",
    "### Example <a name=\"example\"></a>\n",
    "\n",
    "Let us set the security parameter to $\\epsilon=10^{-8}$, say that we have min-entropy sources with linear rates $k_1=n\\cdot0.8$ and $k_2=(n-1)\\cdot0.8$, and ask for $m=100$ fully random bits. According to the formulae above, $n=253$ together with $n-1=252$ weakly random bits can be condensed into $100$ fully random bits (up to the security parameter $\\epsilon=10^{-8}$).\n",
    "\n",
    "Next, we give an efficient implementation of this Toeplitz based construction.\n",
    "\n",
    "### Implementation <a name=\"implementation\"></a>\n",
    "\n",
    "The vector-matrix multiplication $x\\cdot(T(y)|1_m)^T$ a priori has asymptotic complexity $O(n^2)$ in big O-notation, which is prohibitive for larger input sizes $n\\geq10^4$. However, we discuss in the theory notes [1] that the operation is actually implemented in asymptotic complexity $O(n\\log n)$ by first embedding the problem into circulant matrices and then making use of the Fast Fourier Transform (FFT). The corresponding code then performs well up to input sizes $n\\geq10^7$. The following example demonstrates this implementation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-29T21:34:12.002044Z",
     "start_time": "2023-08-29T21:34:11.948208Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Security parameter: 1e-08.\n",
      "Desired output length: 10 bits.\n",
      "Min-entropy rate of first source: 0.8.\n",
      "Min-entropy rate of second source: 0.8\n",
      "Required length of each input source: 103 bits.\n",
      "The 10 random output bits are:\n",
      "[[1 0 1 0 0 1 1 1 0 1]].\n"
     ]
    }
   ],
   "source": [
    "# work with local simulator for testing Toeplitz construction\n",
    "device = LocalSimulator()\n",
    "\n",
    "# set security parameter\n",
    "power = 8\n",
    "eps = 10 ** (-power)\n",
    "print(f\"Security parameter: {eps}.\")\n",
    "\n",
    "# set number of output bits\n",
    "m = 10\n",
    "print(f\"Desired output length: {m} bits.\")\n",
    "\n",
    "# set min-entropy rates for sources\n",
    "k_1 = 0.8\n",
    "print(f\"Min-entropy rate of first source: {k_1}.\")\n",
    "k_2 = 0.8\n",
    "print(f\"Min-entropy rate of second source: {k_2}\")\n",
    "\n",
    "# required number of input bits (for each source)\n",
    "n = math.floor((m - 1 - 2 * math.log2(eps)) / (k_1 + k_2 - 1))\n",
    "print(f\"Required length of each input source: {n} bits.\")\n",
    "\n",
    "# quantum circuit for generating weakly random bit string one\n",
    "n1_qubits = 1\n",
    "m1_shots = n\n",
    "state1 = hadamard_circuit(n1_qubits)\n",
    "result1 = device.run(state1, shots=m1_shots).result()\n",
    "array_one = result1.measurements.reshape(1, m1_shots * n1_qubits)\n",
    "# print(array_one)\n",
    "\n",
    "# quantum circuit for generating weakly random bit string two\n",
    "n2_qubits = 1\n",
    "m2_shots = n\n",
    "state2 = hadamard_circuit(n2_qubits)\n",
    "result2 = device.run(state2, shots=m2_shots).result()\n",
    "array_two = result2.measurements.reshape(1, m2_shots * n2_qubits)\n",
    "# print(array_two)\n",
    "\n",
    "###\n",
    "# alternative for generating two bit strings when no quantum source is available:\n",
    "\n",
    "# create first list of pseudo-random bits\n",
    "# alternative when no quantum source is available\n",
    "# list_one = []\n",
    "# for number in range(n):\n",
    "#    b = int(random.randint(0, 1))\n",
    "#    list_one.append(b)\n",
    "# array_one = np.array([list_one])\n",
    "\n",
    "# create second list of pseudo-random bits\n",
    "# list_two = []\n",
    "# for number in range(n):\n",
    "#    b = int(random.randint(0, 1))\n",
    "#    list_two.append(b)\n",
    "# array_two = np.array([list_two])\n",
    "###\n",
    "\n",
    "# computing output of Toeplitz extractor by vector-matrix multiplication\n",
    "# via efficient Fast Fourier Transform (FFT) as discussed in [1]\n",
    "\n",
    "# setting up arrays for FFT implementation of Toeplitz\n",
    "array_two_under = np.array(array_two[0, 0 : n - m])[np.newaxis]\n",
    "zero_vector = np.zeros((1, n + m - 3), dtype=int)\n",
    "array_two_zeros = np.hstack((array_two_under, zero_vector))\n",
    "array_two_over = array_two[0, n - m : n][np.newaxis]\n",
    "array_one_merged = np.zeros((1, 2 * n - 3), dtype=int)\n",
    "for i in range(m):\n",
    "    array_one_merged[0, i] = array_one[0, m - 1 - i]\n",
    "for j in range(n - m - 1):\n",
    "    array_one_merged[0, n + m - 2 + j] = array_one[0, n - 2 - j]\n",
    "\n",
    "# FFT multplication output of Toeplitz\n",
    "output_fft = np.around(ifft(fft(array_one_merged) * fft(array_two_zeros)).real)\n",
    "output_addition = output_fft[0, 0:m] + array_two_over\n",
    "output_final = output_addition.astype(int) % 2\n",
    "print(f\"The {m} random output bits are:\\n{output_final}.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As an alternative, we note that efficient implementations of other quantum-proof two-source extractors are discussed in [3].\n",
    "\n",
    "## Unpredictability of physical sources <a name=\"physical_sources\"></a>\n",
    "\n",
    "Given above methods on randomness extraction, the next step is to give lower bounds on the min-entropy present in the output distributions generated from our $n$-fold Hadamard circuit. For that, we need to model the noise present in the quantum processing units.\n",
    "\n",
    "Generally, for any given quantum processing unit, the supplier typically publishes some type of noise specification with it. This includes both, the noise characterization of state preparation, as well as the read-out measurements. In case such specifications are not available, or if one wants to double check them, it is in principle possible to benchmark the device. We refer to the theory notes [1] for more details on this and just mention here that we do not need a full characterization of the device, but rather only conservative upper bounds on the noise strength. This then translates into lower bounds on the min-entropy present in the system.\n",
    "\n",
    "For our case, since we are only applying single qubit gates for our Hadamard circuit, the noise is captured well by single qubit noise models. Moreover, for the state preparation step via Hadamard gates, the typical noise in quantum architectures is uniform, depolarizing noise of some strength $\\lambda\\in[0,1]$. That is, all possible single qubit errors such as bit flip or phase flip errors are equally likely, leading to the effective evolution\n",
    "\n",
    "$$ \\psi=|\\psi\\rangle\\langle\\psi|\\mapsto\\text{Dep}^\\lambda(\\psi)=(1-\\lambda)\\cdot\\psi+\\lambda\\cdot\\frac{|0\\rangle\\langle0|+|1\\rangle\\langle1|}{2}, $$\n",
    "\n",
    "mapping any input qubit state $\\psi$ onto a linear combination of itself and the maximally mixed qubit state $\\frac{|0\\rangle\\langle0|+|1\\rangle\\langle1|}{2}$. Note that we now work with general mixed states in order to model classical statistical uncertainty coming from the noise model.\n",
    "\n",
    "So effectively, before the measurement step, instead of the perfect pure state $|H\\rangle\\langle H|_A$ as defined by the vector $|H\\rangle_A=\\frac{1}{\\sqrt{2}}\\big(|0\\rangle_A+|1\\rangle_A\\big)$, we have the mixed state\n",
    "\n",
    "$$ \\rho_A^\\lambda=\\frac{1}{2}\\big(|0\\rangle\\langle0|_A+(1-\\lambda)|0\\rangle\\langle1|_A+(1-\\lambda)|1\\rangle\\langle0|_A+|1\\rangle\\langle1|_A\\big) $$\n",
    "\n",
    "at hand.\n",
    "\n",
    "Next, we note that instead of the ideal measurement device given by $\\mathcal{M}=\\big\\{|0\\rangle\\langle0|_A,|1\\rangle\\langle1|_A\\big\\}$, the typical noisy measurement device is described by \n",
    "\n",
    "$$ \\mathcal{N}^\\mu=\\big\\{1_A-\\mu|1\\rangle\\langle1|_A,\\mu|1\\rangle\\langle1|_A\\big\\} $$\n",
    "\n",
    "with some bias $\\mu\\in(0,1)$ towards reading-out the ground state $|0\\rangle\\langle0|_A$ over $|1\\rangle\\langle1|_A$. The post-measurement probability distribution is then given as\n",
    "\n",
    "$$ Q^{\\lambda,\\mu}=(q_0=1-\\mu/2,q_1=\\mu/2) $$\n",
    "\n",
    "instead of the perfectly uniform distribution $P=(p_0=1/2,p_1=1/2)$. Note that the former distribution has non-maximal min-entropy\n",
    "\n",
    "$$ H_{\\min}(X)_{Q^{\\lambda,\\mu}}=1-\\log\\mu\\leq1=H_{\\min}(X)_P. $$\n",
    "\n",
    "More generally, as measurement devices are the most sensitive element of quantum randomness generation, we discuss in the theory notes [1] methods for benchmarking them (even when only noisy state preparation is available).\n",
    "\n",
    "### Noise as leakage <a name=\"noise_leakage\"></a>\n",
    "\n",
    "It is, however, crucial to realize that $H_{\\min}(X)_{Q^{\\lambda,\\mu}}$ is not yet the quantity relevant for secure quantum randomness generation. As information is never lost in quantum mechanics, all the noise carries information to the environment, where it can in principle be picked up by an attacker. That is, we need to estimate the conditional min-entropy of the post-measurement probability distribution $X$ given any complementary information that leaked into the environment [5].\n",
    "\n",
    "This is worked out in detail in the theory notes [1] by means of so-called purifications of both, the noisy qubit state $\\rho_A^\\lambda$, as well as the noisy measurement device $\\mathcal{N}^\\mu$, leading to the additional qubit registers $A'$ and $E_2$, respectively. The relevant so-called classical-quantum state to consider then takes the form\n",
    "\n",
    "$$ \\omega_{XA'E_2}^{\\lambda,\\mu}= q_0\\cdot|0\\rangle\\langle0|_X\\otimes\\omega_{A'E_2}^{\\lambda,\\mu}(0) + q_1\\cdot|1\\rangle\\langle1|_X\\otimes\\omega_{A'E_2}^{\\lambda,\\mu}(1) $$\n",
    "\n",
    "with a classing part $X$ and the quantum parts $A'E_2$ depending on both noise parameters $\\lambda,\\mu$. The corresponding conditional min-entropy is\n",
    "\n",
    "$$ H_{\\min}(X|A'E_2)_{\\omega^{\\lambda,\\mu}}\\leq H_{\\min}(X)_{Q^{\\lambda,\\mu}} $$\n",
    "\n",
    "with the gap between the conditional and the unconditional term typically being strict.\n",
    "\n",
    "We reiterate that the reason for mathematically introducing the purification registers $AE_2$ and working with the conditional min-entropy, is to make sure that the output bits generated are random even conditioned on the knowledge of the noise. This ensures that the randomness created is from a purely quantum origin and is secure against any eavesdropping from a malicious third party.\n",
    "\n",
    "We mention that a more detailed discussion of this point is given in the work [5].\n",
    "\n",
    "### Numerical evaluation <a name=\"numerical_evaluation\"></a>\n",
    "\n",
    "In the following, we use the noise parameters $\\lambda=0.02$ and $\\mu=0.98$ for the quantum processing units and this immediately gives\n",
    "\n",
    "$$H_{\\min}(X)\\approx0.944. $$\n",
    "\n",
    "The next step is to numerically evaluate the conditional min-entropy\n",
    "\n",
    "$$ H_{\\min}(X|A'E_2)=-\\log p_{\\text{guess}}(X|A'E_2). $$\n",
    "\n",
    "This is done by means of a standard semidefinite program (sdp) solver as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-29T21:34:14.333898Z",
     "start_time": "2023-08-29T21:34:13.802934Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "===============================================================================\n",
      "                                     CVXPY                                     \n",
      "                                     v1.6.0                                    \n",
      "===============================================================================\n",
      "(CVXPY) Feb 07 04:31:36 PM: Your problem has 16 variables, 48 constraints, and 0 parameters.\n",
      "(CVXPY) Feb 07 04:31:36 PM: It is compliant with the following grammars: DCP, DQCP\n",
      "(CVXPY) Feb 07 04:31:36 PM: (If you need to solve this problem multiple times, but with different data, consider using parameters.)\n",
      "(CVXPY) Feb 07 04:31:36 PM: CVXPY will first compile your problem; then, it will invoke a numerical solver to obtain a solution.\n",
      "(CVXPY) Feb 07 04:31:36 PM: Your problem is compiled with the CPP canonicalization backend.\n",
      "-------------------------------------------------------------------------------\n",
      "                                  Compilation                                  \n",
      "-------------------------------------------------------------------------------\n",
      "(CVXPY) Feb 07 04:31:36 PM: Compiling problem (target solver=SCS).\n",
      "(CVXPY) Feb 07 04:31:36 PM: Reduction chain: Complex2Real -> Dcp2Cone -> CvxAttr2Constr -> ConeMatrixStuffing -> SCS\n",
      "(CVXPY) Feb 07 04:31:36 PM: Applying reduction Complex2Real\n",
      "(CVXPY) Feb 07 04:31:36 PM: Applying reduction Dcp2Cone\n",
      "(CVXPY) Feb 07 04:31:36 PM: Applying reduction CvxAttr2Constr\n",
      "(CVXPY) Feb 07 04:31:36 PM: Applying reduction ConeMatrixStuffing\n",
      "(CVXPY) Feb 07 04:31:36 PM: Applying reduction SCS\n",
      "(CVXPY) Feb 07 04:31:36 PM: Finished problem compilation (took 1.280e-02 seconds).\n",
      "-------------------------------------------------------------------------------\n",
      "                                Numerical solver                               \n",
      "-------------------------------------------------------------------------------\n",
      "(CVXPY) Feb 07 04:31:36 PM: Invoking solver SCS  to obtain a solution.\n",
      "------------------------------------------------------------------\n",
      "\t       SCS v3.2.5 - Splitting Conic Solver\n",
      "\t(c) Brendan O'Donoghue, Stanford University, 2012\n",
      "------------------------------------------------------------------\n",
      "problem:  variables n: 32, constraints m: 108\n",
      "cones: \t  s: psd vars: 108, ssize: 3\n",
      "settings: eps_abs: 1.0e-05, eps_rel: 1.0e-05, eps_infeas: 1.0e-07\n",
      "\t  alpha: 1.50, scale: 1.00e-01, adaptive_scale: 1\n",
      "\t  max_iters: 100000, normalize: 1, rho_x: 1.00e-06\n",
      "\t  acceleration_lookback: 10, acceleration_interval: 10\n",
      "lin-sys:  sparse-direct-amd-qdldl\n",
      "\t  nnz(A): 168, nnz(P): 0\n",
      "------------------------------------------------------------------\n",
      " iter | pri res | dua res |   gap   |   obj   |  scale  | time (s)\n",
      "------------------------------------------------------------------\n",
      "     0| 3.89e+00  1.00e+00  1.56e+01 -7.07e+00  1.00e-01  7.52e-04 \n",
      "   250| 9.44e-05  2.51e-05  1.63e-05  6.07e-01  8.20e+00  2.74e-03 \n",
      "   375| 2.45e-07  1.99e-05  1.38e-07  6.08e-01  2.59e+01  3.86e-03 \n",
      "------------------------------------------------------------------\n",
      "status:  solved\n",
      "timings: total: 3.86e-03s = setup: 4.98e-04s + solve: 3.36e-03s\n",
      "\t lin-sys: 4.77e-04s, cones: 2.51e-03s, accel: 2.19e-05s\n",
      "------------------------------------------------------------------\n",
      "objective = 0.607508\n",
      "------------------------------------------------------------------\n",
      "-------------------------------------------------------------------------------\n",
      "                                    Summary                                    \n",
      "-------------------------------------------------------------------------------\n",
      "(CVXPY) Feb 07 04:31:36 PM: Problem status: optimal\n",
      "(CVXPY) Feb 07 04:31:36 PM: Optimal value: 6.075e-01\n",
      "(CVXPY) Feb 07 04:31:36 PM: Compilation took 1.280e-02 seconds\n",
      "(CVXPY) Feb 07 04:31:36 PM: Solver (including time spent in interface) took 4.227e-03 seconds\n",
      "\u001b[1mThe coditional min-entropy is:  0.7190235372097399\n",
      "As a comparison, the unconditional min-entropy is:  0.9439714610772487\n"
     ]
    }
   ],
   "source": [
    "# fix noise parameters\n",
    "lamb = 0.02\n",
    "mu = 0.98\n",
    "\n",
    "# purification of rho input state\n",
    "rho = 0.5 * np.array([[1, 1 - lamb], [1 - lamb, 1]])\n",
    "eigvals, eigvecs = np.linalg.eig(rho)\n",
    "\n",
    "rho_vector = (\n",
    "    math.sqrt(eigvals[0]) * np.kron(eigvecs[:, 0], eigvecs[:, 0])[np.newaxis]\n",
    "    + math.sqrt(eigvals[1]) * np.kron(eigvecs[:, 1], eigvecs[:, 1])[np.newaxis]\n",
    ")\n",
    "rho_pure = np.kron(rho_vector, rho_vector.T)\n",
    "\n",
    "# sigma state of noisy measurement device\n",
    "sigma_vector = np.array([[math.sqrt(1 - mu), 0, 0, math.sqrt(mu)]])\n",
    "sigma_pure = np.kron(sigma_vector, sigma_vector.T)\n",
    "\n",
    "# omega state relevant for conditional min-entropy\n",
    "rho_sigma = np.kron(rho_pure, sigma_pure)\n",
    "\n",
    "id_2 = np.identity(2)\n",
    "zero = np.array([[1, 0]])\n",
    "one = np.array([[0, 1]])\n",
    "\n",
    "zerozero = np.kron(np.kron(zero, id_2), np.kron(zero, id_2))\n",
    "zeroone = np.kron(np.kron(zero, id_2), np.kron(one, id_2))\n",
    "onezero = np.kron(np.kron(one, id_2), np.kron(zero, id_2))\n",
    "oneone = np.kron(np.kron(one, id_2), np.kron(one, id_2))\n",
    "\n",
    "omega_0 = (\n",
    "    zerozero @ rho_sigma @ zerozero.T\n",
    "    + zeroone @ rho_sigma @ zeroone.T\n",
    "    + onezero @ rho_sigma @ onezero.T\n",
    ")\n",
    "omega_1 = oneone @ rho_sigma @ oneone.T\n",
    "\n",
    "omega = []\n",
    "omega.append(omega_0)\n",
    "omega.append(omega_1)\n",
    "\n",
    "# sdp solver\n",
    "m = 4  # dimension of quantum side information states\n",
    "c = 2  # number of classical measurement outcomes\n",
    "\n",
    "sigma = cp.Variable((m, m), complex=True)  # complex variable\n",
    "constraints = [sigma >> 0]  # positive semi-definite\n",
    "constraints += [sigma >> omega[i] for i in range(c)]  # min-entropy constraints\n",
    "obj = cp.Minimize(cp.real(cp.trace(sigma)))  # objective function\n",
    "\n",
    "prob = cp.Problem(obj, constraints)  # set up sdp problem\n",
    "prob.solve(solver=cp.SCS, verbose=True)  # solve sdp problem using splitting conic solver (SCS)\n",
    "guess = prob.value\n",
    "qmin_entropy = (-1) * math.log2(guess)\n",
    "min_entropy = 1 - math.log2(2 - mu * (1 - lamb))\n",
    "\n",
    "print(\"\\033[1m\" + \"The coditional min-entropy is: \", qmin_entropy)\n",
    "print(\"As a comparison, the unconditional min-entropy is: \", min_entropy)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That is, for the chosen noise parameters $\\lambda=0.02$ and $\\mu=0.98$ we find\n",
    "\n",
    "$$ H_{\\min}(X|A'E_2)\\approx0.719<0.944\\approx H_{\\min}(X).$$\n",
    "\n",
    "By varying the noise parameter to other values, one also sees by inspection that the conditional min-entropy is monotone in both $\\lambda$ and $\\mu$. Importantly, this ensures that the outputted randomness will be safe to use whenever we put a conservative estimate on the noisy strength, even in the absence of an exact characterization of the underlying quantum processing units.\n",
    "\n",
    "Finally, whenever we run our Hadamard circuit $n$ times or on $n$ qubits in parallel, the overall conditional min-entropy is just given by\n",
    "\n",
    "$$ H_{\\min}(X|RE_2)_{\\omega^{\\otimes n}} = n\\cdot H_{\\min}(X|RE_2)_\\omega. $$\n",
    "\n",
    "This reason is that we only consider single qubit product noise together with the min-entropy being additive on product states. This then leads to the promised linear min-entropy rates $k_i=\\alpha_i\\cdot n$ for $i=1,2$ going into the Toeplitz two-source extractor.\n",
    "\n",
    "As an added bonus, the Toeplitz two-source extractor used has the so-called strong property. That is, even conditioned on the knowledge of one input string of weakly random bits, the output bits are still fully random. We refer to the technical notes [1] for discussion, a consequence being that if one provider builds in some backdoors in the provided unit, the random generation scheme is still not broken unless the two providers come together to cooperate.\n",
    "\n",
    "\n",
    "## Putting everything together <a name=\"putting-together\"></a>\n",
    "\n",
    "Now that we have determined the conditional min-entropy of our physical sources and have an efficient quantum-proof two-source extractor in place, all that remains is to put the two pieces together.\n",
    "\n",
    "1. First, we specify how many random bits we want to generate, the desired security parameter, and the conditional min-entropy of our weak sources of randomness from the quantum processing units:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-29T21:34:16.379045Z",
     "start_time": "2023-08-29T21:34:16.360499Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Security parameter: 1e-08.\n",
      "Desired output length: 10 bits.\n",
      "Min-entropy rate of first source: 0.8.\n",
      "Min-entropy rate of second source: 0.8.\n",
      "Required length of each input source: 141 bits.\n"
     ]
    }
   ],
   "source": [
    "# set security parameter\n",
    "power = 8\n",
    "eps = 10 ** (-power)\n",
    "print(f\"Security parameter: {eps}.\")\n",
    "\n",
    "# set number of output bits\n",
    "m = 10\n",
    "print(f\"Desired output length: {m} bits.\")\n",
    "\n",
    "# set min-entropy rates for sources - qmin_entropy from above\n",
    "k_one = qmin_entropy\n",
    "print(f\"Min-entropy rate of first source: {k_1}.\")\n",
    "k_two = qmin_entropy\n",
    "print(f\"Min-entropy rate of second source: {k_2}.\")\n",
    "\n",
    "# required number of input bits (for each source)\n",
    "n = math.floor((m - 1 - 2 * math.log2(eps)) / (k_one + k_two - 1))\n",
    "print(f\"Required length of each input source: {n} bits.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. At the beginning of the notebook, we loaded two separate QPUs as available in Amazon Braket. We now run on the respective QPUs the Hadamard circuit followed by measurements in the computational basis:\n",
    "\n",
    "    (Note: If preferred, one can alternatively use the pre-loaded on-demand simulator SV1. In this case, please comment out the QPU code in the cell below and instead uncomment the provided SV1 code.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-29T21:34:22.584748Z",
     "start_time": "2023-08-29T21:34:18.293952Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The first raw bit string is:  [[1 1 0 1 1 0 1 1 1 0 0 0 1 0 0 1 1 1 0 1 0 1 0 1 1 1 1 1 0 0 1 0 1 0 1 1\n",
      "  1 0 0 0 1 1 0 1 1 1 1 0 1 1 1 0 0 1 1 1 1 1 1 1 1 1 1 0 0 1 0 1 0 0 0 1\n",
      "  1 1 0 1 1 1 0 0 1 1 0 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 1 1 1 0]]\n",
      "The second raw bit string is:  [[0 0 0 1 0 1 1 1 0 0 0 1 1 1 0 1 0 1 0 0 0 0 1 1 1 1 1 1 1 0 0 1 1 1 0 1\n",
      "  1 0 0 0 1 1 0 1 1 1 0 0 1 1 1 1 1 0 1 0 1 1 0 0 1 0 1 1 0 1 1 0 0 1 1 1\n",
      "  1 0 0 1 0 0 1 1 0 0 0 0 1 1 1 1 0 0 1 1 0 1 0 1 0 0 1 1 1 1 1]]\n"
     ]
    }
   ],
   "source": [
    "# Use on-demand simulator SV1\n",
    "\n",
    "# quantum circuit for generating weakly random bit string one (simulate Rigetti source)\n",
    "n1_q = 1 # alternatively run multiple qubits in parallel\n",
    "m1_s = int(math.ceil(n/n1_q))\n",
    "state1 = hadamard_circuit(n1_q)\n",
    "task1 = simulator.run(state1, shots=m1_s)\n",
    "result1 = task1.result()\n",
    "rigetti_task_id = task1.id\n",
    "array_rigetti = result1.measurements.reshape(1,m1_s*n1_q)\n",
    "print(\"The first raw bit string is: \",array_one)\n",
    "\n",
    "# quantum circuit for generating weakly random bit string two (simulate IonQ source)\n",
    "n2_q = 1 # alternatively run multiple qubits in parallel\n",
    "m2_s = int(math.ceil(n/n2_q))\n",
    "state2 = hadamard_circuit(n2_q)\n",
    "task2 = simulator.run(state2, shots=m2_s)\n",
    "result2 = task2.result()\n",
    "ionq_task_id = task2.id\n",
    "array_ionq = result2.measurements.reshape(1,m2_s*n2_q)\n",
    "print(\"The second raw bit string is: \",array_two)\n",
    "\n",
    "###\n",
    "\n",
    "# [Uncomment the following lines to use the Rigetti and IonQ devices. Be aware that a usage cost will incur when you submit the circuit.]\n",
    "\n",
    "# Rigetti: quantum circuit for generating weakly random bit string one\n",
    "# n1_q = (\n",
    "#     1  # alternatively use more than one qubit (attention: 32 max + lower bounds on number of shots)\n",
    "# )\n",
    "# m1_s = int(math.ceil(n / n1_q))\n",
    "# state_rigetti = hadamard_circuit(n1_q)\n",
    "# rigetti_task = rigetti.run(state_rigetti, shots=m1_s, poll_timeout_seconds=5 * 24 * 60 * 60)\n",
    "\n",
    "# rigetti_task_id = rigetti_task.id\n",
    "# rigetti_status = rigetti_task.state()\n",
    "# print(\"Status of Rigetti task:\", rigetti_status)\n",
    "\n",
    "# IonQ: quantum circuit for generating weakly random bit string two\n",
    "# n2_q = (\n",
    "#     1  # alternatively use more than one qubit (attention: 11 max + lower bounds on number of shots)\n",
    "# )\n",
    "# m2_s = int(math.ceil(n / n2_q))\n",
    "# state_ionq = hadamard_circuit(n2_q)\n",
    "# ionq_task = ionq.run(state_ionq, shots=m2_s, poll_timeout_seconds=5 * 24 * 60 * 60)\n",
    "\n",
    "# ionq_task_id = ionq_task.id\n",
    "# ionq_status = ionq_task.state()\n",
    "# print(\"Status of IonQ task:\", ionq_status)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3. The quantum tasks have now been sent to the respective QPUs and we can recover the results at any point in time once the status is completed:\n",
    "\n",
    "    (Note: In case you opted to use the on-demand simulator SV1 instead of the QPUs, please do not run the following cell.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-29T21:34:23.471216Z",
     "start_time": "2023-08-29T21:34:22.586872Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Status of Rigetti task: COMPLETED\n",
      "The first raw bit string is:  [[0 0 0 0 1 1 1 0 0 0 0 0 0 1 0 1 1 1 1 0 1 0 1 0 0 1 0 1 1 1 0 0 0 1 1 1\n",
      "  1 0 1 0 1 1 0 0 0 1 1 0 1 0 1 1 1 1 1 0 1 0 1 1 0 1 1 0 0 1 0 0 0 0 1 0\n",
      "  0 1 0 1 1 0 0 0 1 1 1 1 1 1 0 0 0 1 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0\n",
      "  0 0 0 0 1 1 1 0 1 1 1 0 1 0 1 1 1 1 1 0 1 0 0 0 1 1 1 0 1 1 1 1 0]]\n",
      "Status of IonQ quantum task: COMPLETED\n",
      "The second raw bit string is:  [[0 1 1 1 0 0 0 0 0 0 0 1 1 0 1 1 0 0 0 0 0 0 0 1 0 1 1 1 1 1 0 1 0 1 1 0\n",
      "  1 0 0 0 0 1 0 1 1 0 0 0 1 1 1 1 1 0 1 1 0 0 0 1 1 1 1 0 0 0 1 1 0 0 1 0\n",
      "  0 0 0 1 1 1 0 1 0 1 0 1 1 0 1 0 1 0 0 0 1 1 1 1 1 1 0 0 1 0 1 1 0 1 0 1\n",
      "  1 1 1 1 1 1 0 0 1 0 1 0 1 0 1 0 0 1 1 1 0 1 0 0 0 1 0 1 1 0 0 1 1]]\n"
     ]
    }
   ],
   "source": [
    "# recover Rigetti quantum task\n",
    "task_load_rigetti = AwsQuantumTask(arn=rigetti_task_id)\n",
    "\n",
    "# print status\n",
    "status_rigetti = task_load_rigetti.state()\n",
    "print(\"Status of Rigetti task:\", status_rigetti)\n",
    "# wait for hybrid job to complete\n",
    "# terminal_states = ['COMPLETED', 'FAILED', 'CANCELLED']\n",
    "if status_rigetti == \"COMPLETED\":\n",
    "    # get results\n",
    "    rigetti_results = task_load_rigetti.result()\n",
    "\n",
    "    # array\n",
    "    array_rigetti = rigetti_results.measurements.reshape(1, m1_s * n1_q)\n",
    "    print(\"The first raw bit string is: \", array_rigetti)\n",
    "\n",
    "elif status_rigetti in [\"FAILED\", \"CANCELLED\"]:\n",
    "    # print terminal message\n",
    "    print(\"Your Rigetti quantum task is in terminal status, but has not completed.\")\n",
    "\n",
    "else:\n",
    "    # print current status\n",
    "    print(\n",
    "        \"Sorry, your Rigetti quantum task is still being processed and has not been finalized yet.\"\n",
    "    )\n",
    "\n",
    "# recover IonQ task\n",
    "task_load_ionq = AwsQuantumTask(arn=ionq_task_id)\n",
    "\n",
    "# print status\n",
    "status_ionq = task_load_ionq.state()\n",
    "print(\"Status of IonQ quantum task:\", status_ionq)\n",
    "# wait for hybrid job to complete\n",
    "# terminal_states = ['COMPLETED', 'FAILED', 'CANCELLED']\n",
    "if status_ionq == \"COMPLETED\":\n",
    "    # get results\n",
    "    ionq_results = task_load_ionq.result()\n",
    "\n",
    "    # array\n",
    "    # print(m1_shots,n1_qubits)\n",
    "    # print(m2_shots,n2_qubits)\n",
    "    # print(ionq_results.measurements)\n",
    "    array_ionq = ionq_results.measurements.reshape(1, m2_s * n2_q)\n",
    "    print(\"The second raw bit string is: \", array_ionq)\n",
    "\n",
    "elif status_ionq in [\"FAILED\", \"CANCELLED\"]:\n",
    "    # print terminal message\n",
    "    print(\"Your IonQ quantum task is in terminal status, but has not completed.\")\n",
    "\n",
    "else:\n",
    "    # print current status\n",
    "    print(\"Sorry, your IonQ quantum task is still being processed and has not been finalized yet.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "4. We run the Toeplitz two-source extractor on the two sequences of raw random input bits:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-29T21:34:23.483829Z",
     "start_time": "2023-08-29T21:34:23.479594Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The 10 random output bits are:\n",
      "[[1 0 0 1 0 1 1 1 1 0]].\n"
     ]
    }
   ],
   "source": [
    "# setting up arrays for fft implementation of Toeplitz\n",
    "if status_ionq == \"COMPLETED\":\n",
    "    array_two_under = np.array(array_ionq[0, 0 : n - m])[np.newaxis]\n",
    "    zero_vector = np.zeros((1, n + m - 3), dtype=int)\n",
    "    array_two_zeros = np.hstack((array_two_under, zero_vector))\n",
    "    array_two_over = array_ionq[0, n - m : n][np.newaxis]\n",
    "    array_one_merged = np.zeros((1, 2 * n - 3), dtype=int)\n",
    "    if status_rigetti == \"COMPLETED\":\n",
    "        for i in range(m):\n",
    "            array_one_merged[0, i] = array_rigetti[0, m - 1 - i]\n",
    "        for j in range(n - m - 1):\n",
    "            array_one_merged[0, n + m - 2 + j] = array_rigetti[0, n - 2 - j]\n",
    "\n",
    "        # fft multplication output of Toeplitz\n",
    "        output_fft = np.around(ifft(fft(array_one_merged) * fft(array_two_zeros)).real)\n",
    "        output_addition = output_fft[0, 0:m] + array_two_over\n",
    "        output_final = output_addition.astype(int) % 2\n",
    "        print(f\"The {m} random output bits are:\\n{output_final}.\")\n",
    "    else:\n",
    "        print(f\"Your Rigetti quantum task is in {status_rigetti} state.\")\n",
    "else:\n",
    "    print(f\"Your IonQ quantum task is in {status_ionq} state.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "5. That is it, above bits are random up to the chosen security parameter!\n",
    "\n",
    "## Beyond current implementation <a name=\"beyond_current\"></a>\n",
    "\n",
    "Different quantum processing units are potentially governed by different noise models. Correspondingly, this will lead to different conditional min-entropy rates for the raw sources of randomness. In this notebook and following the theory notes [1], we can change the noise model for the numerical evaluation of the min-entropy (with all other parts remaining the same).\n",
    "\n",
    "From a cryptographic viewpoint, the quantum hardware providers as well as Amazon Braket have to be trusted in order to end up with secure random numbers. More broadly, we note that there are more intricate ways of generating randomness from quantum sources than presented here, in particular for the use case when one is not ready to trust the underlying quantum hardware because of potential backdoors inbuilt. Such a so-called (semi) device-independent scheme is for example given in [3], with a corresponding implementation in Amazon Braket [7]. For a general in-depth discussion of QRNGs, we refer to the review article [4], as well as the extensive security report [6].\n",
    "\n",
    "\n",
    "## Literature <a name=\"literature\"></a>\n",
    "\n",
    "[1] M. Berta and F. Brandão, Robust randomness generation on quantum computers, [available online](http://marioberta.info/wp-content/uploads/2021/07/randomness-theory.pdf).\n",
    "\n",
    "[2] M. Hayashi and T. Tsurumaru, More Efficient Privacy Amplification With Less Random Seeds via Dual Universal Hash Function, IEEE Transactions on Information Theory, [10.1109/TIT.2016.2526018](https://ieeexplore.ieee.org/document/7399404).\n",
    "\n",
    "[3] C. Foreman, S. Wright, A. Edgington, M. Berta, F. Curchod, Practical randomness and privacy amplification, [arXiv:2009.06551](https://arxiv.org/abs/2009.06551).\n",
    "\n",
    "[4] M. Herrero-Collantes and J.-C. Garcia-Escartin, Quantum random number generators, Review of Modern Physics, [10.1103/RevModPhys.89.015004](https://doi.org/10.1103/RevModPhys.89.015004).\n",
    "\n",
    "[5] D. Frauchiger, R. Renner, M. Troyer, True randomness from realistic quantum devices, [arXiv:1311.4547](https://arxiv.org/abs/1311.4547).\n",
    "\n",
    "[6] M. Piani, M. Mosca, B. Neill, Quantum random-number generators: practical considerations and use cases, [evolutionQ](https://evolutionq.com/quantum-safe-publications/qrng-report-2021-evolutionQ.pdf).\n",
    "\n",
    "[7] Quantum-Proof Cryptography with IronBridge, TKET and Amazon Braket, A. Edgington, C. Foreman, D. Jones, [Cambridge Quantum Computing](https://medium.com/cambridge-quantum-computing/quantum-proof-cryptography-with-ironbridge-tket-and-amazon-braket-e8e96777cacc)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-08-29T21:34:23.950740Z",
     "start_time": "2023-08-29T21:34:23.487815Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Quantum Task Summary\n",
      "{<_Amazon.SV1: 'arn:aws:braket:::device/quantum-simulator/amazon/sv1'>: {'shots': 282, 'tasks': {'COMPLETED': 2}, 'execution_duration': datetime.timedelta(microseconds=4000), 'billed_execution_duration': datetime.timedelta(seconds=6)}}\n",
      "Note: Charges shown are estimates based on your Amazon Braket simulator and quantum processing unit (QPU) task usage. Estimated charges shown may differ from your actual charges. Estimated charges do not factor in any discounts or credits, and you may experience additional charges based on your use of other services such as Amazon Elastic Compute Cloud (Amazon EC2).\n",
      "Estimated cost to run this example: 0.01 USD\n"
     ]
    }
   ],
   "source": [
    "print(\"Quantum Task Summary\")\n",
    "print(t.quantum_tasks_statistics())\n",
    "print(\n",
    "    \"Note: Charges shown are estimates based on your Amazon Braket simulator and quantum processing unit (QPU) task usage. Estimated charges shown may differ from your actual charges. Estimated charges do not factor in any discounts or credits, and you may experience additional charges based on your use of other services such as Amazon Elastic Compute Cloud (Amazon EC2).\"\n",
    ")\n",
    "print(\n",
    "    f\"Estimated cost to run this example: {t.qpu_tasks_cost() + t.simulator_tasks_cost():.2f} USD\"\n",
    ")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.15"
  },
  "vscode": {
   "interpreter": {
    "hash": "590fab68195cf107911461461f81d5c472d3d6127f579badfcfad30f03e5cab2"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}