samples/algorithms/noisy-amp-est/NoisyAmpEst.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "da44ff43-687f-4d4c-bb4f-c44ef21cd188",
   "metadata": {
    "nteract": {
     "transient": {
      "deleting": false
     }
    }
   },
   "source": [
    "# Noisy Amplitude Estimation with low-depth algorithms\n",
    "\n",
    "## Introduction \n",
    "This notebook deals with a low-depth algorithm to estimate the amplitude of a unitary operation $U$ on a noisy system following these articles:\n",
    "* [[Tiron et al.]](https://arxiv.org/abs/2012.03348v1): Tudor Giurgica-Tiron, Iordanis Kerenidis, Farrokh Labib, Anupam Prakash, and William Zeng (2022), \"Low depth algorithms for quantum amplitude estimation\", _Quantum,_ Volume 6, page 745;  arXiv:2012.03348\n",
    "\n",
    "* [[Tanaka et al.]](https://link.springer.com/article/10.1007/s11128-021-03215-9) Tomoki Tanaka, Yohichi Suzuki, Shumpei Uno, Rudy Raymond, Tamiya Onodera & Naoki Yamamoto (2021), \"Amplitude estimation via maximum likelihood on noisy quantum computer\", _Quantum Information Processing_, Volume 20, Issue 9, Article 293; arXiv:2006.16223\n",
    "\n",
    "* [[Suzuki et al.]](http://arxiv.org/abs/1904.10246v2): \n",
    "Yohichi Suzuki, Shumpei Uno, Rudy Raymond, Tomoki Tanaka, Tamiya Onodera, and Naoki Yamamoto (2020): \"Amplitude estimation without phase estimation\", _Quantum Information Processing,_ Volume 19, Issue 2, Article 75; arXiv:1904.10246\n",
    "\n",
    "* [[Brassard et al.]](https://arxiv.org/abs/quant-ph/0005055v1): Gilles Brassard, Peter Hoyer, Michele Mosca, and Alain Tapp (2002), \"Quantum amplitude amplification and estimation\", _Contemporary Mathematics_, Volume 305, pp. 53-74; arXiv:0005055v1"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "0c762ab1",
   "metadata": {},
   "source": [
    "## Background\n",
    "\n",
    "### 1. Amplitude Estimation\n",
    "\n",
    "Given a unitary $U$ that acts on $(n+1)$ qubits such that $U:\\ket{0}^n \\ket{0} \\mapsto \\cos \\theta\\ \\ket{\\Psi_0}\\ket{0} + \\sin \\theta \\ket{\\Psi_1} \\ket{1}$ as an oracle with an unknown angle $\\theta\\in[0,\\pi/2]$, the amplitude estimation problem aims to find an angle $\\hat{\\theta}$ such that $|\\hat{\\theta} - \\theta| \\leq \\varepsilon$ for some precision parameter $\\varepsilon$. Doing this classically, i.e. based on the probability $\\sin^2(\\theta)$,  would require calling the oracle $O(\\frac{1}{\\varepsilon^2})$ times. The quantum amplitude estimation algorithm uses the Quantum Fourier Transform as a sub-routine and estimates the amplitude using only $O(\\frac{1}{\\varepsilon})$ queries. However, obtaining this quadratic speedup requires a high circuit depth of $O(\\frac{1}{\\varepsilon})$ and in the presence of noise, the accuracy this circuit cannot be guaranteed.\n",
    "\n",
    "On the other hand, by settling for a sub-quadratic speedup, it is possible to design quantum-classical hybrid algorithms that require only a constant depth on noisy systems. This class of algorithms use amplitude amplification followed by maximum likelihood estimation (MLE) to compute $\\hat{\\theta}$. This Q\\# + Python notebook demonstrates one such low-depth algorithm following the work of [Tiron et al.](https://arxiv.org/abs/2012.03348v1), [Tanaka et al.](https://link.springer.com/article/10.1007/s11128-021-03215-9) and [Suzuki et al.](http://arxiv.org/abs/1904.10246v2). By adjusting the various parameters in this notebook it is possible to reproduce the Power Law based algorithm in the last reference as well as the Linear and Exponential schedule algorithms in the first one."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "858bf92d-697c-4aa2-9fd2-2dc3f992a352",
   "metadata": {
    "nteract": {
     "transient": {
      "deleting": false
     }
    }
   },
   "source": [
    "### 2. Maximum Likelihood Estimation\n",
    "\n",
    "As the name suggests, Maximum Likelihood Estimation (MLE) aims to answer the question: *\"Given a sequence of i.i.d samples from an unknown distribution parameterized by a variable, say $\\alpha$, what value of $\\alpha$ is the most likely to produce the observed samples?\"*\n",
    "\n",
    "Formally, for a random variable $X$, let $\\text{Pr}(X, \\alpha)$ be the density function parametrized by an unknown parameter $\\alpha$. Now, given $n$ independent samples from this distribution such that $X_i = x_i$ for $i \\in \\{1,\\dots,n\\}$, the likelihood of $\\alpha$ taking a value $\\theta$ is given by the *likelihood function*:\n",
    "$$ L(\\theta|x) = \\Pi_{i=1}^n \\text{Pr} (X_i=x_i | \\alpha = \\theta). $$\n",
    "For ease of notation, we will denote the likelihood function as $L(\\theta)$ and the dependency on the samples $x_i$ is clear from context.\n",
    "MLE suggests that the $\\theta$ that maximizes the above equation is a good estimate for the ideal value $\\theta^*$ actually taken by $\\alpha$. To avoid working with small probabilities, it is common to instead maximize the *log-likelihood* function $\\ell(\\theta) = \\log L(\\theta)$. Then, a good estimate for alpha will be:\n",
    "$$ \\hat{\\theta} = \\arg \\max_{\\theta} L(\\theta) = \\arg \\max_{\\theta} \\ell(\\theta).$$\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "bbd8abdd",
   "metadata": {},
   "source": [
    "## Setup\n",
    "\n",
    "First, we import the requisite packages, load the Q# environment, and the necessary Q# libraries needed to run the code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5506ce60",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-08-02T19:22:30.807085Z",
     "start_time": "2022-08-02T19:22:30.801082Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from numpy import *\n",
    "import math\n",
    "import traceback\n",
    "import sys"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eccfd963",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-08-02T19:22:52.847157Z",
     "start_time": "2022-08-02T19:22:31.680257Z"
    }
   },
   "outputs": [],
   "source": [
    "import qsharp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a3c22d05",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-08-02T19:22:55.217194Z",
     "start_time": "2022-08-02T19:22:52.849155Z"
    },
    "microsoft": {
     "language": "qsharp"
    }
   },
   "outputs": [],
   "source": [
    "%%qsharp\n",
    "open Microsoft.Quantum.Random;\n",
    "open Microsoft.Quantum.Arrays;\n",
    "open Microsoft.Quantum.Measurement;\n",
    "open Microsoft.Quantum.Intrinsic;\n",
    "open Microsoft.Quantum.Canon;\n",
    "open Microsoft.Quantum.Convert;\n",
    "open Microsoft.Quantum.Math;\n",
    "open Microsoft.Quantum.Diagnostics;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0cf3ce3a-4c20-4bca-8099-a711a0584f81",
   "metadata": {
    "nteract": {
     "transient": {
      "deleting": false
     }
    }
   },
   "source": [
    "## Amplitude Amplification: a refresher\n",
    "\n",
    "For an unknown angle $\\theta^* \\in [0,\\pi/2]$, define the unitary $U$ acting on $(n+1)$ qubits with the action\n",
    "$$ U: \\ket{0}^n \\ket{0} \\mapsto \\ket{\\Psi} := \\sin \\theta^* \\ket{\\Psi_1}\\ket{1} + \\cos \\theta^* \\ket{\\Psi_0}\\ket{0}$$\n",
    "where $\\ket{\\Psi_1}$ labeled as a *good* state. \n",
    "In this sample, we set the default angle as $\\theta^* = \\pi / (2 \\cdot \\Phi) \\approx 0.9708$ where $\\Phi = (1 + \\sqrt{5})/2$ is the Golden Ratio."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a9f8b863",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%qsharp\n",
    "/// # Summary\n",
    "/// Returns the value of the unknown angle.\n",
    "/// Angle value should be in [0, PI()/2] \n",
    "function UnknownTHETA(): Double {\n",
    "    return PI()/(1.0 + Sqrt(5.0));\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6ceb46a1",
   "metadata": {},
   "source": [
    "📝 **Remark:** Modify the value of the angle in the previous cell to try the sample for different values of $\\theta^*$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "330f5ff7",
   "metadata": {},
   "source": [
    "\n",
    "The following operation describes one such unitary acting on $2$ qubits applying an arbitrary rotation for a fixed angle $\\theta^*$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "298805b0",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-08-02T19:22:56.432212Z",
     "start_time": "2022-08-02T19:22:55.219157Z"
    },
    "microsoft": {
     "language": "qsharp"
    }
   },
   "outputs": [],
   "source": [
    "%%qsharp\n",
    "/// # Summary \n",
    "/// The U operation. \n",
    "/// Unitary takes |0>|0> -> cos theta* |0>|0> + sin theta* |1>|1>\n",
    "/// This operation is accomplished by applying Ry(-2*theta)\n",
    "///\n",
    "/// # Inputs\n",
    "/// ## qs\n",
    "/// The 2 qubits register on which the oracle is applied.\n",
    "operation ArbitRotation (qs : Qubit[]) : Unit is Adj + Ctl {\n",
    "    \n",
    "    Ry(2.0*UnknownTHETA(), qs[0]); // Takes |0>|0> -> cos theta* |0>|0> + sin theta* |1>|0>\n",
    "    CNOT(qs[0], qs[1]); // Takes cos theta* |0>|0> + sin theta* |1>|0> -> cos theta* |0>|0> + sin theta* |1>|1>\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1919279a",
   "metadata": {},
   "source": [
    "📝 **Remark:** To run the sample for a different unitary operation, modify the above cell accordingly."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc494d1c-25ca-4a03-99c2-7bf61b4da5cd",
   "metadata": {
    "nteract": {
     "transient": {
      "deleting": false
     }
    }
   },
   "source": [
    "By repeatedly measuring the state $\\ket{\\Psi}$ and counting the fraction of $1$ outcomes would give an estimate for $\\sin^2 \\theta^*$. While this is a low-depth circuit, it does not provide any speedup over classical methods. \n",
    "\n",
    "However, using quantum amplitude amplification, it is possible to *amplify* the probability of seeing a $1$ outcome. Then, fewer queries can be used to obtain a good estimate for $\\sin^2 \\theta^*$ and, by extension, $\\theta^*$. This is done by applying the operator \n",
    "$$ \\mathbf{Q} = - U S_0 U^{\\dagger} S_{\\chi} $$\n",
    "where $S_{\\chi}$ marks the *good* states and $S_0$ marks the all-zeroes state. In other words, \n",
    "$ S_{\\chi} \\ket{\\Psi_1}\\ket{1} = -  \\ket{\\Psi_1}\\ket{1}$ and other states remain unchanged while $S_{0} \\ket{0}^{n+1} = - \\ket{0}^{n+1}$ and other states remain unchanged. The following operations describe how to construct the amplification operator $\\mathbf{Q}$ by first constructing the $S_0$ and $S_{\\chi}$ operations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dc160c26",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-08-02T19:22:57.767209Z",
     "start_time": "2022-08-02T19:22:56.434213Z"
    },
    "microsoft": {
     "language": "qsharp"
    }
   },
   "outputs": [],
   "source": [
    "%%qsharp\n",
    "/// #Summary \n",
    "/// The S_0 operation. Applies -1 to |0>|0> state and nothing to others. \n",
    "/// Can be easily extended to act on qubit registers of any length.\n",
    "///\n",
    "/// #Inputs\n",
    "/// ## register\n",
    "/// The 2-qubit register on which the operation acts.  \n",
    "operation PhaseFlip_AllZeros(register : Qubit[]) : Unit is Adj + Ctl {\n",
    "    use aux = Qubit();\n",
    "    (ControlledOnInt(0, X))(register, aux);\n",
    "    Z(aux);\n",
    "    (ControlledOnInt(0, X))(register, aux);\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0ae2ada7",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-08-02T19:22:59.657322Z",
     "start_time": "2022-08-02T19:22:57.770215Z"
    },
    "microsoft": {
     "language": "qsharp"
    }
   },
   "outputs": [],
   "source": [
    "%%qsharp\n",
    "/// #Summary\n",
    "/// The S_chi operation. Applies -1 to the good state i.e., |1>|1> and nothing to others.\n",
    "/// Can be extended to act on registers of arbitrary length and a specific good state |\\Psi>\n",
    "///\n",
    "/// #Inputs\n",
    "/// ## register\n",
    "/// The 2-qubits register on which the marking operation acts\n",
    "operation PhaseFlip_Marking(register : Qubit[]) : Unit is Adj + Ctl {\n",
    "    \n",
    "    use aux = Qubit();\n",
    "    within {\n",
    "        //convert target to |0> - |1>\n",
    "        X(aux);\n",
    "        H(aux);\n",
    "    }\n",
    "    apply {\n",
    "        //Get phase kickback for |1>|1> state\n",
    "        Controlled X (register, aux);\n",
    "        \n",
    "    }\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9c11b0c2",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-08-02T19:23:01.352210Z",
     "start_time": "2022-08-02T19:22:59.665224Z"
    },
    "microsoft": {
     "language": "qsharp"
    }
   },
   "outputs": [],
   "source": [
    "%%qsharp\n",
    "/// #Summary\n",
    "/// The amplification iterate Q = - U S_0 U^dag S_chi\n",
    "///\n",
    "/// #Inputs\n",
    "/// ## register\n",
    "/// The 2-qubit register on which the amplification iterate is applied.\n",
    "operation AmplificationIterate(register : Qubit[]) : Unit is Adj + Ctl {\n",
    "    PhaseFlip_Marking(register); //S_chi\n",
    "    Adjoint ArbitRotation(register); // U^dag    \n",
    "    PhaseFlip_AllZeros(register); //S_0    \n",
    "    ArbitRotation(register); // U    \n",
    "    R(PauliI, 2.0*PI(), register[0]);    // Applies -1 global phase\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ed1acfac-712a-4be7-a42f-e9f4265229c0",
   "metadata": {
    "nteract": {
     "transient": {
      "deleting": false
     }
    }
   },
   "source": [
    "For $\\theta^* \\in [0, \\frac{\\pi}{2}]$, and $\\ket{\\Psi}$ obtained by applying $U$ once, [Brassard et al.](http://arxiv.org/abs/quant-ph/0005055v1) showed that $m$ applications of $\\mathbf{Q}$ gives: \n",
    "$$ \\mathbf{Q}^m :\\ket{\\Psi} \\mapsto \\sin ((2m+1) \\theta^*) \\ket{\\Psi_1} \\ket{1} + \\cos ((2m+1) \\theta^*) \\ket{\\Psi_0} \\ket{0} \\tag{1}$$\n",
    "Measuring the above state gives $1$ with a probability that is about $4m^2$ larger than measuring $\\ket{\\Psi}$ with only $2m$ more queries to $U$. This intuitively leads to a possible quadratic quantum speedup. The following operation describes how to raise a given unitary to the $m^{\\text{th}}$ power."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7b5f4439",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-08-02T19:23:02.732215Z",
     "start_time": "2022-08-02T19:23:01.354211Z"
    },
    "microsoft": {
     "language": "qsharp"
    }
   },
   "outputs": [],
   "source": [
    "%%qsharp\n",
    "/// #Summary\n",
    "/// Given a multi-qubit Unitary V, and a power m apply V^m\n",
    "///\n",
    "/// #Inputs\n",
    "/// ## power\n",
    "/// The integer power to raise the unitary to.\n",
    "/// ## target\n",
    "/// The qubit register on which the unitary is to be applied.\n",
    "/// ## V\n",
    "/// The unitary operation to be applied.\n",
    "operation PowUnitary (power : Int, V : (Qubit[] => Unit is Adj + Ctl), target : Qubit[]) : Unit is Adj + Ctl {\n",
    "    Fact(power >= 0, \"power should be a positive integer.\");\n",
    "    for i in 1..power {\n",
    "        V(target);\n",
    "    }\n",
    "}"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "7befb68c",
   "metadata": {},
   "source": [
    "## Simulating noise"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "0a0f4ccd",
   "metadata": {},
   "source": [
    "We consider a commonly used noise model: *the depolarizing noise* as considered in [Tanaka et al.](https://link.springer.com/article/10.1007/s11128-021-03215-9) and [Tiron et al.](https://arxiv.org/abs/2012.03348v1). Its action is described as a channel $\\mathcal{D}$ on an $n$-qubit density matrix $\\rho$ as: \n",
    "$$ \\mathcal{D} (\\rho) = p \\rho + (1 - p) \\frac{\\mathbb{I}}{2^n}$$\n",
    "Intuitively this means that with probability $(1-p)$, the qubits are depolarized and remained unchanged otherwise. Similarly, the `AmplificationIterate` operation $\\mathbf{Q}$ can be described as a channel $\\mathcal{Q}$ with:\n",
    "$$ \\mathcal{Q} (\\rho) = \\mathbf{Q} \\rho \\mathbf{Q}^\\dag$$\n",
    "\n",
    "As the most commonly repeated part of the amplitude amplification algorithm is the `AmplificationIterate` operation, we consider the interaction of the noise channel $\\mathcal{D}$ with the $\\mathcal{Q}$ channel to create the noisy amplification iterate $\\mathcal{DQ}$ where\n",
    "$$ \\mathcal{DQ}(\\rho) = \\mathcal{QD}(\\rho) = p \\mathbf{Q} \\rho \\mathbf{Q}^\\dag + (1 - p) \\frac{\\mathbb{I}}{2^n}.$$\n",
    "Extending this to $m$ repetitions of the noisy amplification process gives\n",
    "$$ (\\mathcal{QD})^m(\\rho) = \\mathcal{D}^m\\mathcal{Q}^m(\\rho) = p^m \\mathbf{Q}^m \\rho \\mathbf{Q}^{\\dag m} + (1 - p^m) \\frac{\\mathbb{I}}{2^n}. \\tag{2}$$\n",
    "In other words, with probability $p^m$, apply $m$ repetitions of $\\mathbf{Q}$, otherwise, the qubits depolarize. This action is obtained by defining the `AmplitudeAmplification` procedure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e396a4bd",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%qsharp\n",
    "/// #Summary\n",
    "/// The quantum operations to apply the noisy amplification channel mvar times and measure the output\n",
    "///\n",
    "/// #Inputs\n",
    "/// ## mvar\n",
    "/// The power to which the amplification channel DQ is raised to.\n",
    "///\n",
    "/// ## depol_prob\n",
    "/// The probability with which the qubits are depolarized\n",
    "///\n",
    "/// #Outputs\n",
    "/// ## result\n",
    "operation AmplitudeAmplification(mvar : Int, depol_prob : Double) : Bool { \n",
    "    \n",
    "    let p = 1.0 - depol_prob;    // qubits depolarize with probability depol_prob\n",
    "    let prob = DrawRandomDouble(0.0, 1.0);\n",
    "\n",
    "    use qs = Qubit[2]; // initialize qubits\n",
    "    \n",
    "    if prob <= PowD(p, IntAsDouble(mvar)) {\n",
    "        ArbitRotation(qs); // Output: cos theta/2 |00> + sin theta/2 |11>\n",
    "        PowUnitary(mvar, AmplificationIterate, qs); // Output: cos (2m+1)theta/2 |00> + sin (2m+1)theta/2 |11>\n",
    "    }\n",
    "    else {    // if the qubits are depolarized, reset all the qubits. \n",
    "        ResetAll(qs);\n",
    "        H(qs[1]);    // Set the that will be measured to output 1 with probability 0.5 mimicking depolarized behavior\n",
    "    }\n",
    "    \n",
    "    let res = M(qs[1]); //Measure qubit\n",
    "    \n",
    "    ResetAll(qs); //release qubits cleanly\n",
    "    \n",
    "    let result = ResultAsBool(res);\n",
    "    return result;\n",
    "}"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "2e949146-932a-439b-b0ff-9d797cbb6fb9",
   "metadata": {
    "nteract": {
     "transient": {
      "deleting": false
     }
    }
   },
   "source": [
    "## Quantum-classical Amplitude Estimation in noisy systems\n",
    "\n",
    "The quantum-classical approach to Amplitude Estimation, applies the noisy amplification channel $\\mathcal{DQ}$ $m$ times and then measures the output qubit to record a $0$ and $1$ measurement. For angle $\\theta^*$, setting $\\gamma = - \\ln p$, the probability of obtaining a $1$ outcome after $m$ repetitions of the noisy amplitude amplification procedure is $ \\frac{1}{2} - \\frac{1}{2} e^{-\\gamma m} \\cos(2 (2m+1)\\theta^*)$. Repeating this procedure for $N_{\\text{shots}}$ shots where $1$ is observed $N_1$ times and $0$ is observed $N_0$ times (such that $N_0 + N_1 = N_{\\text{shots}}$), the likelihood of the unknown angle being $\\theta$ is\n",
    "$$ \\Pr(X | \\alpha = \\theta) := \\left[\\frac{1}{2} - \\frac{1}{2} e^{-\\gamma m} \\cos(2 (2m+1)\\theta)\\right]^{N_1} \\left[\\frac{1}{2} + \\frac{1}{2} e^{-\\gamma m} \\cos(2 (2m+1)\\theta)\\right]^{N_0}.$$\n",
    "\n",
    "Taking $K$ independent samples for different values of $m_k$, for $k = 1, \\ldots, K$ the likelihood function becomes\n",
    "$$ L(\\theta) = \\Pi_{k=1}^K   \\left[\\frac{1}{2} - \\frac{1}{2} e^{-\\gamma m_k} \\cos(2 (2m_k+1)\\theta)\\right]^{N_{k,1}} \\left[\\frac{1}{2} + \\frac{1}{2} e^{-\\gamma m_k} \\cos(2 (2m_k+1)\\theta)\\right]^{N_{k,0}}.$$\n",
    "Equivalently, we can use the log-likelihood function and compute \n",
    "$$ \\ell(\\theta) = \\sum_{k=1}^K \\left[N_{k,1} \\log \\left(\\frac{1}{2} - \\frac{1}{2} e^{-\\gamma m_k} \\cos(2 (2m_k+1)\\theta)\\right) + N_{k,0} \\log \\left(\\frac{1}{2} + \\frac{1}{2} e^{-\\gamma m_k} \\cos(2 (2m_k+1)\\theta)\\right)\\right].$$\n",
    "\n",
    "Computing this log-likelihood function required \n",
    "$$\n",
    "N_{\\text{shots}}\\sum_{k=1}^K (2m_k + 1)\n",
    "$$\n",
    "usages of the operator $U$. \n",
    "\n",
    "As there is no prior information, initially we assume that any choice for $\\theta \\in [0, \\pi/2]$ could be equally likely. \n",
    "So, the log-likelihood function is initialized to be the the all-zero function. \n",
    "\n",
    "Additionally, while $\\theta$ is a continuous parameter, as we only need to provide an estimate up to a precision of $\\varepsilon$, the interval $[0, \\pi/2]$ is divided into $\\lceil \\frac{1}{\\varepsilon} \\rceil$ discrete parts that then form the possible choices for $\\theta$.\n",
    "This log-likelihood function for these choices is initialized below. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c056408f",
   "metadata": {},
   "source": [
    "📝 **Note:** As subsequent code in this sample now focuses on the classical part of the quantum-classical hybrid, it is expressed in Python."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b5cc2658",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-08-02T19:23:04.307215Z",
     "start_time": "2022-08-02T19:23:04.294213Z"
    }
   },
   "outputs": [],
   "source": [
    "# #Summary\n",
    "# Initialize the log likelihood function to log of uniform probability over choices for theta.\n",
    "#\n",
    "# #Inputs\n",
    "# ## eps\n",
    "# The precision parameter used to discretize the choice of angles in [0, pi/2] \n",
    "#\n",
    "# #Outputs\n",
    "# ## init_llhood\n",
    "# The initial log-likelihood function set to a uniform value over choices for theta\n",
    "def init_log_likelihood(eps):\n",
    "    if(eps <= 0 or eps > 1):\n",
    "        raise ValueError(\"Precision eps should be a number in (0, 1]\")\n",
    "    \n",
    "    arr_len = math.ceil(1/eps) \n",
    "    norm = arr_len*eps\n",
    "    init_llhood = np.full(arr_len, 0.0)\n",
    "    \n",
    "    return init_llhood"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "40366e90-7809-407b-a7d6-66b8e87bd263",
   "metadata": {
    "nteract": {
     "transient": {
      "deleting": false
     }
    }
   },
   "source": [
    "\n",
    "The Bayesian update for the log-likelihood function at each iteration is computed as\n",
    "$$ \\ell_{\\text{new}}(\\theta) = \\ell_{\\text{old}}(\\theta) + N_{k,1} \\log \\left(\\frac{1}{2} - \\frac{1}{2} e^{-\\gamma m_k} \\cos(2 (2m_k+1)\\theta)\\right) + N_{k,0} \\log \\left(\\frac{1}{2} + \\frac{1}{2} e^{-\\gamma m_k} \\cos(2 (2m_k+1)\\theta)\\right) \\tag{3}$$\n",
    "which is described in the method below.\n",
    "\n",
    "**Note:** We have to be careful to deal with the exceptions when we deal with the logarithm of $0$ here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0f438dcf",
   "metadata": {},
   "outputs": [],
   "source": [
    "# #Summary\n",
    "# Compute the Bayesian update of the likelihood function for the current iteration as per Equation (2)\n",
    "#\n",
    "# #Inputs\n",
    "# ## eps \n",
    "# The precision parameter\n",
    "# \n",
    "# ## m \n",
    "# The power to which the amplitude amplification iterate was raised to\n",
    "# \n",
    "# ## N0 \n",
    "# The number of 0 outcomes observed for this iteration\n",
    "# \n",
    "# ## N1 \n",
    "# The number of 1 outcomes observed for this iteration\n",
    "# \n",
    "# ## llhood \n",
    "# An array holding the previous values of the log likelihood function for each choice of theta where for index t we have theta_t = pi*t*eps/2\n",
    "#\n",
    "# ## depol\n",
    "# Probability with which depolarizing noise is applied to the qubits\n",
    "#\n",
    "# #Output\n",
    "# ## llhood \n",
    "# Return the updated value of the log-likelihood function\n",
    "\n",
    "def update_log_likelihood(eps, m, N0, N1, llhood, depol):\n",
    "    length = len(llhood)\n",
    "    one_val=0\n",
    "    zero_val=0\n",
    "\n",
    "    p = 1 - depol \n",
    "    \n",
    "    if(length != math.ceil(1/eps)):\n",
    "        raise IndexError(\"Array length mismatch with expected value.\")\n",
    "\n",
    "    for t in range(length):\n",
    "        val = llhood[t]\n",
    "        angle = 2*(2*m+1)*t*math.pi*eps/2\n",
    "        cosval = math.cos(angle)    # cos (2 (2m+1) theta)\n",
    "        if m == 0:    # to deal with the classical schedule when m is 0\n",
    "            expval = 1\n",
    "        else:\n",
    "            expval = math.pow(p, m)\n",
    "        if cosval == expval:    # when the second term in eq. 3 could be -infinity\n",
    "            if N1 == 0:    # to deal with edge case of 0 * -math.inf\n",
    "                one_val = 0\n",
    "            else:\n",
    "                one_val = -math.inf\n",
    "        elif cosval == -expval:    # when the third term in eq. 3 could be -infinity\n",
    "            if N0 == 0:    # to deal with edge case of 0 * -math.inf\n",
    "                zero_val = 0\n",
    "            else:\n",
    "                zero_val = -math.inf\n",
    "        else:\n",
    "            one_val = N1*math.log(0.5 - 0.5*cosval*expval)\n",
    "            zero_val = N0*math.log(0.5 + 0.5*cosval*expval)\n",
    "        \n",
    "        llhood[t] = val+one_val+zero_val\n",
    "    \n",
    "    return llhood"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f0d80251-ade8-430c-a4e0-bfc97b7d63ea",
   "metadata": {
    "nteract": {
     "transient": {
      "deleting": false
     }
    }
   },
   "source": [
    "Finally, the estimate for the angle is determined as \n",
    "$$ \\hat{\\theta} = \\arg\\max_{\\theta} L(\\theta) = \\arg\\max_{\\theta} \\ell(\\theta).$$\n",
    "\n",
    "This naturally leads to the question of choosing an appropriate $K$ and values for $m_k$ to minimize the depth of the circuits while generating a good estimate for the angle to within a given precision of $\\varepsilon$. This is because the circuit depth scales as as the number of oracle calls for one run of the algorithm which is given by $(2m_k + 1)$ for iteration $k$. Since $m_k$ monotonically increases with $k$, the maximum number of oracle calls in the circuits is $(2 m_K + 1)$ [Tiron et al.](https://arxiv.org/abs/2012.03348v1) use a power-law schedule that for parameter $\\beta \\in (0, 1]$:\n",
    "1. Makes $N = O(1/\\varepsilon^{1+\\beta})$ total oracle calls\n",
    "2. Uses circuits of depth at most $D = O(1/\\varepsilon^{1 - \\beta})$\n",
    "3. Sets number of amplifications $m_k = \\left\\lfloor k^{\\tfrac{1 - \\beta}{2\\beta}} \\right\\rfloor$ for each iteration $k$\n",
    "\n",
    "Observe that the possible values $\\beta \\in (0,1]$ covers the gamut of all possible quantum-classical strategies for amplitude estimation from \n",
    "- fully quantum (without Quantum Fourier Transform) at $\\beta \\leftarrow 0$; to\n",
    "- purely classical sampling at $\\beta = 1$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "99dbfd9b",
   "metadata": {},
   "source": [
    "The Python method below interacts with the Q# code from the previous sections of the notebook to generate one sample for the the power law amplitude estimation algorithm and updates the likelihood function accordingly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e18a09b1-8d15-4798-bd64-485d77908512",
   "metadata": {
    "jupyter": {
     "outputs_hidden": false,
     "source_hidden": false
    },
    "nteract": {
     "transient": {
      "deleting": false
     }
    }
   },
   "outputs": [],
   "source": [
    "# #Summary\n",
    "# Execute one run of the Hybrid Amplitude Estimation algorithm.\n",
    "#\n",
    "# #Inputs\n",
    "# ## angle\n",
    "# The ideal angle used for the oracle and used here for verification purposes\n",
    "# \n",
    "# ## K\n",
    "# The number of independent samples used to generate the log-likelihood function\n",
    "# \n",
    "# ## power\n",
    "# The power to which iteration index k is raised to obtain circuit depth m_k.\n",
    "# Related to power law parameter beta and defined as (1 - beta)/2*beta\n",
    "# \n",
    "# ## log_likelihood\n",
    "# An array holding the log-likelihood function for each possible choice for theta\n",
    "# \n",
    "# ## eps\n",
    "# Accuracy parameter\n",
    "# \n",
    "# ## Nshots\n",
    "# Number of times Amplitude Amplification is repeated to estimate the likelihood function for each iteration\n",
    "#\n",
    "# ## depol\n",
    "# Probability with which depolarizing noise is applied to the qubits\n",
    "def amp_est_iteration(K, power, llhood, eps, Nshots, depol):\n",
    "\n",
    "    queries = 0    # Computed total number of queries made\n",
    "\n",
    "    for k in range(1, K+1):\n",
    "    \n",
    "        N0 = 0\n",
    "        N1 = 0\n",
    "\n",
    "        if power == math.inf:    # To execute an exponential schedule i.e., m_k = 2^k\n",
    "            m = math.floor(pow(2, k))\n",
    "        elif power == 0:    # To execute a classical schedule i.e., m_k = 0\n",
    "            m = 0\n",
    "        else:\n",
    "            m = math.floor(pow(k, power))\n",
    "        \n",
    "        for i in range(Nshots):\n",
    "            res = AmplitudeAmplification.simulate(mvar=m, depol_prob=depol)    # Apply quantum operations and collect measurement result in res\n",
    "            if res == False:\n",
    "                N0 = N0 + 1\n",
    "            else:\n",
    "                N1 = N1 + 1\n",
    "        \n",
    "        llhood = update_log_likelihood(eps, m, N0, N1, llhood, depol)    # Update log-likelihood function\n",
    "        \n",
    "        queries = queries + (2*m+1)*Nshots;    # Calculates total number of queries made in this iteration\n",
    "\n",
    "        if K > 10:\n",
    "            interval = math.floor(K/10)\n",
    "        else:    # Account for K < 10\n",
    "            interval = 1\n",
    "\n",
    "        if(k % interval == 0):    # Displaying intermediate estimates at every {K/10}^th iteration\n",
    "            print(f\"Iteration k = {k}:\")\n",
    "            print(f\"m_k = {m}, N0_k = {N0}, N1_k = {N1}\")\n",
    "            temp_t = np.argmax(llhood)\n",
    "            temp_angle = math.pi*temp_t*eps/2\n",
    "            print(f\"Current estimate for \\u03B8*: {round(temp_angle, 4)}\")    # Printing theta to 4 bits of precision\n",
    "            \n",
    "            # Computing additive error of current estimate\n",
    "            err = (UnknownTHETA() - temp_angle)\n",
    "            print(f\"Additive error of current estimate: {round(err, 5)}\\n\")    # Printing error to 5 bits of precision\n",
    "\n",
    "    return llhood, queries"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d7cf5f94-33d5-424c-bba7-35e63293080d",
   "metadata": {
    "nteract": {
     "transient": {
      "deleting": false
     }
    }
   },
   "source": [
    "The main method that executes the power law amplitude estimation algorithm is described below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0baf5506",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-08-02T19:23:28.774051Z",
     "start_time": "2022-08-02T19:23:28.748232Z"
    }
   },
   "outputs": [],
   "source": [
    "# #Summary\n",
    "# Execute the power law amplitude estimation algorithm for a specific choice of parameters\n",
    "#\n",
    "# #Inputs\n",
    "# ## beta\n",
    "# Power law parameter\n",
    "#\n",
    "# ## Nshots\n",
    "# Number of amplitude amplification repetitions performed to accurately estimate the likelihood in each iteration\n",
    "#\n",
    "# ## eps\n",
    "# Accuracy parameter\n",
    "#\n",
    "# ## depol\n",
    "# Probability with which depolarizing noise is applied to the qubits\n",
    "#\n",
    "# #Outputs\n",
    "# ## theta_hat\n",
    "# The final estimate for the unknown angle output by the algorithm\n",
    "\n",
    "def power_law_amp_est(beta, Nshots, eps, depol):\n",
    "\n",
    "    if (Nshots <= 0 or not type(Nshots) is int):\n",
    "        raise ValueError(\"Nshots should be an integer > 0\")\n",
    "    if(beta < 0 or beta > 1):\n",
    "        raise ValueError(\"Power law parameter beta should be a number in [0, 1]\")\n",
    "    if(eps <= 0 or eps > 1):\n",
    "        raise ValueError(\"Precision eps should be a number in (0, 1]\")\n",
    "    if(depol < 0 or depol > 1):\n",
    "        raise ValueError(\"Probability of depolarizing qubits should be in [0, 1]\")\n",
    "    if(UnknownTHETA() < 0 or UnknownTHETA() > math.pi/2):\n",
    "        raise ValueError(\"Unknown angle in UnknownTHETA() should be in [0, pi/2]\")\n",
    "    \n",
    "    log_likelihood = init_log_likelihood(eps)    # Initialize distribution\n",
    "    N0=0\n",
    "    N1=0\n",
    "    queries=0\n",
    "    queries_expected=0\n",
    "    \n",
    "    K = math.ceil(max(math.log(1/eps), 1/pow(eps, 2*beta)))    # Compute K \n",
    "\n",
    "    if beta == 0:    # Accounting for the edge case of executing an exponential schedule\n",
    "        power = math.inf\n",
    "        max_depth = 2*math.floor(pow(2, K))+1\n",
    "        for k in range(1, K+1):\n",
    "            queries_expected= queries_expected + Nshots*(2*math.floor(pow(2, k))+1)\n",
    "    elif beta == 1:    # Accounting for the edge case of executing a classical schedule\n",
    "        max_depth = 1\n",
    "        power = 0\n",
    "        for k in range(1, K+1):\n",
    "            queries_expected= queries_expected + Nshots*(max_depth)\n",
    "    else:\n",
    "        power = (1 - beta)/(2*beta)\n",
    "        max_depth = math.floor(pow(K, power)) \n",
    "        for k in range(1, K+1):\n",
    "            queries_expected= queries_expected + Nshots*(2*math.floor(pow(k, power))+1)\n",
    "    \n",
    "    length = len(log_likelihood)\n",
    "    print(\"\\033[1mAlgorithm Metrics \\033[0m\")\n",
    "    print(f\"Unknown angle \\u03B8*: {round(UnknownTHETA(), 4)}\")\n",
    "    print(f\"Number of samples K: {K}\")\n",
    "    \n",
    "    if(beta == 0):\n",
    "        print(\"Schedule: Exponential, i.e. m_k = 2^k, for k = 1,...,K\")\n",
    "    elif (beta == 1):\n",
    "        print(\"Schedule: Classical i.e. m_k = 0 for k = 1,...,K\")\n",
    "    elif round(beta-1/3, 4) == 0:\n",
    "        print(\"Schedule: Linear, i.e. m_k = k, for k = 1,...,K\")\n",
    "    else:\n",
    "        print(f\"Schedule: Power law with \\u03B2 = {round(beta, 4)}, i.e. m_k = \\u230Ak^((1 - \\u03B2)/2\\u03B2)\\u230B = \\u230Ak^{round((1 - beta)/(2*beta),4)}\\u230B, for k = 1,...,{K}\")\n",
    "\n",
    "    print(f\"Maximum depth: O(2 m_K + 1); and (2 m_K + 1) = {max_depth}\")\n",
    "    print(f\"Resolution of choices for \\u03B8*: {length}\")\n",
    "    print(f\"Number of oracle calls: {queries_expected}\")\n",
    "    print(f\"Probability of depolarizing qubits: {depol}\\n\")\n",
    "    \n",
    "    print(\"\\033[1mIntermediate Estimates \\033[0m\")\n",
    "    log_likelihood, queries = amp_est_iteration(K, power, log_likelihood, eps, Nshots, depol)\n",
    "    \n",
    "    # Compute final estimate    \n",
    "    print(\"\\n\\033[1mFinal estimates \\033[0m\")\n",
    "    t = np.argmax(log_likelihood)\n",
    "    theta_hat = math.pi*t*eps/2\n",
    "    print(f\"Unknown angle \\u03B8*: {round(UnknownTHETA(), 4)}\") \n",
    "    print(f\"Estimate for \\u03B8*: {round(theta_hat, 4)}\")    # Printing to 4 bits of precision.\n",
    "    print(f\"Total number of oracle calls: {queries}\")\n",
    "    \n",
    "    # Compute additive error \n",
    "    # The estimate is within precision limits if additive error is in [-eps, +eps].\n",
    "    err = (UnknownTHETA() - theta_hat)\n",
    "    \n",
    "    print(f\"Additive error of final estimate of \\u03B8*: {round(err, 5)}\")    # Printing error to 5 bits of precision\n",
    "\n",
    "    return theta_hat"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6dee5a1a",
   "metadata": {},
   "source": [
    "## Examples"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "55c29739",
   "metadata": {},
   "source": [
    "First, we execute the algorithm for an exponential schedule with $\\beta = 0$, $N_{\\text{shots}} = 50$, $\\varepsilon = 0.005$ and $20\\%$ depolarizing noise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f60bdf91",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Executing the exponential schedule with 20% depolarizing noise\n",
    "theta_hat = power_law_amp_est(beta=0, Nshots=50, eps=0.005, depol = 0.2)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "d4ff8edc",
   "metadata": {},
   "source": [
    "The key observation here is that the depolarizing noise is such that the exponential schedule does not provide an accurate estimate for $\\theta^*$. \n",
    "\n",
    "📝 **Remark:** It is possible that one run of the cell above does give an accurate estimate of $\\theta^*$. However, running the cell multiple times, one can observe that this is a rare occurrence. This can be attributed to the probabilistic nature of the noise model considered. In particular, in the *rare case* when there is almost no noise present in a run, the exponential schedule gives a good estimate."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "f12bfa0f",
   "metadata": {},
   "source": [
    "To overcome this uncertainty in getting a good estimate, we instead try the power law schedule with $\\beta = 0.5$."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "0a9377bb",
   "metadata": {},
   "source": [
    "🛑 The computation below can take up to a minute to complete."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "174aae9e",
   "metadata": {},
   "outputs": [],
   "source": [
    "theta_hat = power_law_amp_est(beta=0.5, Nshots=50, eps=0.005, depol=0.2)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "8014e674",
   "metadata": {},
   "source": [
    "📝 **Remark:** By running the above cell multiple observe that, with very high probability, each run gives an accurate estimate of $\\theta^*$."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "56bbddae",
   "metadata": {},
   "source": [
    "From these two runs, we can draw the following conclusion. Although the exponential schedule has fewer samples, the higher depth circuit in the presence of noise in unable to provide an accurate estimate for $\\theta^*$. On the other hand, although the power law schedule require many more samples and oracle calls, the relatively low depth of the circuits helps one obtain an accurate estimate for $\\theta^*$."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "6b16425e",
   "metadata": {},
   "source": [
    "Further, we can back up the need for the hybrid approach of the power law schedule by contrasting its requirements (in terms of samples and oracle calls) with that of a fully classical sampling approach where $\\beta = 1$. As the full execution of this classical schedule (with all other parameters unchanged) will take more than 100 times longer than the previous power law schedule, we only give the algorithm metrics for this case below."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "2a2d6b68",
   "metadata": {},
   "source": [
    "**Algorithm Metrics**\n",
    "```\n",
    "Unknown angle θ*:0.9708\n",
    "Number of samples K: 40000\n",
    "Schedule: Classical i.e. m_k = 0 for k = 1,...,K\n",
    "Maximum depth: O(2 m_K + 1); and (2 m_K + 1) = 1\n",
    "Resolution of choices for θ*: 200\n",
    "Number of oracle calls: 2000000\n",
    "Probability of depolarizing qubits: 0.2\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4f738fdc",
   "metadata": {},
   "source": [
    "## Next steps\n",
    "\n",
    "This notebook uses a simple oracle to demonstrate the effectiveness of hybrid classical and quantum computing to overcome the challenges faced in noisy intermediate scale systems. Feel free to use this as a starting point for your experiments and further exploration. For example, you can \n",
    "\n",
    "- explore how the results change for different oracles and different input parameters;\n",
    "- deploy this code on a QPU to explore if the hardware results match the above simulations;\n",
    "- perform resource estimation to get a break down of resources needed under different assumptions; "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4231c5f5",
   "metadata": {},
   "source": []
  }
 ],
 "metadata": {
  "kernel_info": {
   "name": "python3"
  },
  "kernelspec": {
   "display_name": "qsharp-env",
   "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.8"
  },
  "nteract": {
   "version": "nteract-front-end@1.0.0"
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  },
  "vscode": {
   "interpreter": {
    "hash": "270cbec9a4b9716ac32a0616dd36b8ed00b906828156309cf54cec975b25649d"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}