{
    "cells": [
     {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
       "# Iterative Phase Estimation via the Cloud\n",
       "\n",
       "This sample code and notebook was written by members of KPMG Quantum team in Australia. It aims to demonstrate expanded capabilities of targets that support integrated hybrid computing. It makes use of bounded loops, classical function calls, nested conditional if statements, mid-program measurements and qubit reuse.\n",
       "\n",
       "## Connect to Azure Quantum and set up target\n",
       "\n",
       "First, we must configure the qsharp module to connect the azure workspace and specify a target. The Quantinuum H1 emulator target `quantinuum.sim.h1-1e` is chosen by default, which performs hardware-modeled noisy simulation. Running the experiments in this notebook against that target will consume approximately 25 HQC.\n",
       "\n",
       "We also configure the target profile as `Adaptive_RI` to indicate we use the QIR profile with support for mid-circuit measurement, measurement-based control flow, and classical integer computation as part of compilation.\n",
       "\n",
       "Replace the `subscription_id`, `resource_group`, `name`, and `location` connection parameters with the values for your configured Azure Quantum Workspace, or try copying the Python code for connecting to your workspace from the VS Code Azure Quantum Extension."
      ]
     },
     {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
       "import qsharp\n",
       "\n",
       "qsharp.init(target_profile=qsharp.TargetProfile.Adaptive_RI)"
      ]
     },
     {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
       "import azure.quantum\n",
       "\n",
       "workspace = azure.quantum.Workspace(\n",
       "    subscription_id = \"SUBSCRIPTION ID\",\n",
       "    resource_group = \"RESOURCE GROUP\",\n",
       "    name = \"WORKSPACE NAME\",\n",
       "    location = \"LOCATION\",\n",
       ")"
      ]
     },
     {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
       "target = workspace.get_targets(\"quantinuum.sim.h1-1e\")"
      ]
     },
     {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
       "## Two Dimensional Inner Product Calculation Using Iterative Phase Estimation on Three Qubits\n",
       "\n",
       "This notebook demonstrates an iterative phase estimation within Q#. It will use iterative phase estimation to calculate an inner product between two 2-dimensional vectors encoded on a target qubit and an ancilla qubit. An additional control qubit is also initialised which will be the only qubit used for measurement.\n",
       "\n",
       "The circuit begins by encoding the pair of vectors on the target qubit and the ancilla qubit. It then applies an Oracle operator to the entire register, controlled off the control qubit (which is set up in the $\\ket +$ state). The controlled Oracle operator generates a phase on the $\\ket 1$ state of the control qubit. This can then be read by applying a H gate to the control qubit to make the phase observable when measuring."
      ]
     },
     {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
       "## Encoding vectors\n",
       "\n",
       "The vectors v and c are to be encoded onto the target qubit and the ancilla qubit. The vector $v = (cos(\\frac{\\theta_1}{2}),sin(\\frac{\\theta_1}{2}))$ can be represented by the quantum state $\\ket v = cos(\\frac{\\theta_1}{2})\\ket 0 + sin(\\frac{\\theta_1}{2})\\ket 1$, similarly $c$ can be constructed using $\\theta_2$. \n",
       "\n",
       "A Y rotation applied to a target qubit in the $\\ket 0$ state:\n",
       "\n",
       "$$RY(\\theta)\\ket 0 = e^{iY\\theta/2}\\ket 0 = cos(\\frac{\\theta}{2})\\ket 0 + sin(\\frac{\\theta}{2})\\ket 1$$\n",
       "\n",
       "**Note**: <u> A factor of 2 </u> is present here on theta. An application of a $RY(2\\pi)$ gate on $\\ket 0$ gives the state $-\\ket 0$ and would encode the vector $(-1,0)$. This phase cannot be considered a global phase and removed as the entire register will be entangled.\n",
       "\n",
       "The register of the target qubit and ancilla qubit is,\n",
       "\n",
       "$$\\ket  \\Psi = \\ket {\\Psi_\\text{Target qubit}}\\ket {\\Psi_\\text{Ancilla qubit}}$$\n",
       "The state to be created is on the target qubit and the ancilla qubit is,\n",
       "\n",
       "$$\\ket{\\Psi}=\\frac{1}{\\sqrt{2}}(\\ket{v}\\ket{+}+\\ket{c}\\ket{-}),$$\n",
       "\n",
       "\n",
       "which also takes the form,\n",
       "\n",
       "$$\\ket{\\Psi} = \\frac{1}{2}(\\ket{v}+\\ket{c})\\ket{0}+\\frac{1}{2}(\\ket{v}-\\ket{c})\\ket{1}.$$\n"
      ]
     },
     {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
       "vscode": {
        "languageId": "qsharp"
       }
      },
      "outputs": [],
      "source": [
       "%%qsharp\n",
       "\n",
       "operation StateInitialisation(\n",
       "    TargetReg : Qubit,\n",
       "    AncilReg : Qubit,\n",
       "    theta1 : Double,\n",
       "    theta2 : Double\n",
       ") : Unit is Adj + Ctl {\n",
       "    H(AncilReg);\n",
       "    // Arbitrary controlled rotation based on theta. This is vector v.\n",
       "    Controlled R([AncilReg], (PauliY, theta1, TargetReg));\n",
       "    // X gate on ancilla to change from |+〉 to |-〉.\n",
       "    X(AncilReg);\n",
       "    // Arbitrary controlled rotation based on theta. This is vector c.\n",
       "    Controlled R([AncilReg], (PauliY, theta2, TargetReg));\n",
       "    X(AncilReg);\n",
       "    H(AncilReg);\n",
       "}"
      ]
     },
     {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
       "## The Oracle\n",
       "\n",
       "An oracle G needs to be constructed such that it generates an eigenphase on the state encoded on the target qubit and the ancilla qubit. The construction of this oracle is unimportant to the demonstration within this notebook, but the operation it applies is,\n",
       "\n",
       "$$G\\ket \\Psi = e^{2\\pi i\\phi} \\ket \\Psi.$$\n",
       "\n",
       "where the inner product $\\braket {v|c}$ is contained within the phase $\\phi$, which is bound between [0,1]. When applied controlled on the control qubit which begins in that state $\\ket{\\Psi_\\text{Control Qubit}} = \\ket +$,\n",
       "\n",
       "$$\\begin{aligned}\n",
       "    \\text{Controlled }G \\ket{\\Psi_\\text{Control Qubit}} \\ket \\Psi  & = \\frac {1}{\\sqrt{2}} (\\ket 0 \\ket \\Psi + e^{2\\pi i\\phi}\\ket 1 \\ket \\Psi )\\\\\n",
       "    & =\\frac {1}{\\sqrt{2}} (\\ket 0 + e^{2\\pi i\\phi}\\ket 1) \\ket \\Psi\n",
       "\\end{aligned}$$\n",
       "\n",
       "Now the control qubit contains the phase $\\phi$ which relates to the inner product $\\braket {v|c}$\n",
       "\n",
       "$$\\ket{\\Psi_\\text{Control Qubit}} = \\frac {1}{\\sqrt{2}} (\\ket 0 + e^{2\\pi i\\phi}\\ket 1)$$"
      ]
     },
     {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
       "vscode": {
        "languageId": "qsharp"
       }
      },
      "outputs": [],
      "source": [
       "%%qsharp\n",
       "\n",
       "operation GOracle(\n",
       "    TargetReg : Qubit,\n",
       "    AncilReg : Qubit,\n",
       "    theta1 : Double,\n",
       "    theta2 : Double\n",
       ") : Unit is Adj + Ctl {\n",
       "    Z(AncilReg);\n",
       "    within {\n",
       "        Adjoint StateInitialisation(TargetReg, AncilReg, theta1, theta2);\n",
       "        X(AncilReg);\n",
       "        X(TargetReg);\n",
       "    } apply {\n",
       "        Controlled Z([AncilReg], TargetReg);\n",
       "    }\n",
       "}"
      ]
     },
     {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
       "## Iteration\n",
       "\n",
       "Now for the iterative part of the circuit. For n measurements, consider that the phase can be represented as a binary value $\\phi$, and that applying $2^n$ oracles makes the nth binary point of the phase observable (through simple binary multiplication, and modulus $2\\pi$). The value of the control qubit can be readout, placed in a classical register and the qubit reset for use in the next iteration. The next iteration applies $2^{n-1}$ oracles, correcting phase on the control qubit dependent on the nth measurement. The state on the control qubit can be represented as,\n",
       "\n",
       "$$ \\ket {\\Psi_{\\text{Control Qubit}}} = \\ket 0 + e^{2\\pi i\\phi}\\ket 1 $$\n",
       "\n",
       "where $\\phi = 0.\\phi_0\\phi_1\\phi_2\\phi_3$...\n",
       "\n",
       "Applying $2^n$ controlled oracles gives the state on the control qubit,\n",
       "\n",
       "$$ G^{2^n}\\ket {\\Psi_{\\text{Control Qubit}}} = \\ket 0 + e^{2\\pi i 0.\\phi_n\\phi_{n+1}\\phi_{n+2}\\phi_{n+3}...}\\ket 1 $$\n",
       "\n",
       "Consider that the phase has no terms deeper than $\\phi_n$ (ie, terms $\\phi_{n+1},\\phi_{n+2}, \\text{etc}$),\n",
       "\n",
       "$$ G^{2^n}\\ket {\\Psi_{\\text{Control Qubit}}} = \\ket 0 + e^{2\\pi i 0.\\phi_n}\\ket 1 $$\n",
       "\n",
       "Now the value $\\phi_n$ can be observed with a H gate and a measurement projecting along the Z axis. Resetting the control qubit and applying the oracle $2^{n-1}$ times,\n",
       "\n",
       "$$ G^{2^{n-1}}\\ket {\\Psi_{\\text{Control Qubit}}} = \\ket 0 + e^{2\\pi i 0.\\phi_{n-1}\\phi_n}\\ket 1 $$\n",
       "\n",
       "Using the previous measured value for $\\phi_n$, the additional binary point can be rotated out.\n",
       "\n",
       "$$ RZ(-2\\pi \\times 0.0\\phi_n)G^{n-1}\\ket {\\Psi_{\\text{Control Qubit}}} = \\ket 0 + e^{2\\pi i 0.\\phi_{n-1}}\\ket 1 $$\n",
       "\n",
       "This process is iteratively applied for some bit precision n to obtain the phase $0.\\phi_0\\phi_1\\phi_2...\\phi_{n}$. The value is stored as a binary value $x = \\phi_0\\phi_1\\phi_2...\\phi_{n}$ as only integers are manipulatable at runtime currently.\n",
       "\n",
       "As the readout tells nothing of either vector, only the inner product between them, the states on the target qubit and ancilla qubit <u>remain in the same state</u> throughout the process!"
      ]
     },
     {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
       "vscode": {
        "languageId": "qsharp"
       }
      },
      "outputs": [],
      "source": [
       "%%qsharp\n",
       "\n",
       "open Microsoft.Quantum.Math;\n",
       "open Microsoft.Quantum.Convert;\n",
       "\n",
       "operation IterativePhaseEstimation(\n",
       "    TargetReg : Qubit,\n",
       "    AncilReg : Qubit,\n",
       "    theta1 : Double,\n",
       "    theta2 : Double,\n",
       "    Measurements : Int\n",
       ") : Int {\n",
       "    use ControlReg = Qubit();\n",
       "    mutable MeasureControlReg = [Zero, size = Measurements];\n",
       "    mutable bitValue = 0;\n",
       "    //Apply to initialise state, this is defined by the angles theta1 and theta2\n",
       "    StateInitialisation(TargetReg, AncilReg, theta1, theta2);\n",
       "    for index in 0..Measurements - 1 {\n",
       "        H(ControlReg);\n",
       "        //Don't apply rotation on first set of oracles\n",
       "        if index > 0 {\n",
       "            //Loop through previous results\n",
       "            for index2 in 0..index - 1 {\n",
       "                if MeasureControlReg[Measurements - 1 - index2] == One {\n",
       "                    //Rotate control qubit dependent on previous measurements and number of measurements\n",
       "                    let angle = -IntAsDouble(2^(index2)) * PI() / (2.0^IntAsDouble(index));\n",
       "                    R(PauliZ, angle, ControlReg);\n",
       "                }\n",
       "            }\n",
       "        }\n",
       "\n",
       "        let powerIndex = (1 <<< (Measurements - 1 - index));\n",
       "        //Apply a number of oracles equal to 2^index, where index is the number or measurements left\n",
       "        for _ in 1..powerIndex {\n",
       "            Controlled GOracle([ControlReg], (TargetReg, AncilReg, theta1, theta2));\n",
       "        }\n",
       "        H(ControlReg);\n",
       "        //Make a measurement mid circuit\n",
       "        set MeasureControlReg w/= (Measurements - 1 - index) <- MResetZ(ControlReg);\n",
       "        if MeasureControlReg[Measurements - 1 - index] == One {\n",
       "            //Assign bitValue based on previous measurement\n",
       "            set bitValue += 2^(index);\n",
       "        }\n",
       "    }\n",
       "    return bitValue;\n",
       "}"
      ]
     },
     {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
       "Finally to calculate the inner product from the measured value,\n",
       "\n",
       "$$\\braket {v|c} = -cos(2\\pi x / 2^n)$$\n",
       "\n",
       "where $x = \\phi_0\\phi_1\\phi_2...\\phi_{n}$. The denominator within the cosine function is to shift the binary point to match the original value $\\phi$.\n",
       "\n",
       "**Note**: For inner product that are not -1 or 1, the solutions are paired with a value difference of $2^{n-1}$. For example for n=3 measurements, the measured bit value of 2 would also have a pair solution of 6. Either of these values produce the same value of the inner product when input as the variable to the even function cosine (resulting in an inner product of 0 in this example).\n",
       "\n",
       "**Note**: For inner product solutions between the discrete bit precision, a distribution of results will be produced based on where the inner product lies between the discrete bit value. "
      ]
     },
     {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
       "vscode": {
        "languageId": "qsharp"
       }
      },
      "outputs": [],
      "source": [
       "%%qsharp\n",
       "\n",
       "operation QuantumInnerProduct(theta1 : Double, theta2 : Double, iterationCount : Int) : Int {\n",
       "    //Create target register\n",
       "    use TargetReg = Qubit();\n",
       "    //Create ancilla register\n",
       "    use AncilReg = Qubit();\n",
       "    //Run iterative phase estimation\n",
       "    let Results = IterativePhaseEstimation(TargetReg, AncilReg, theta1, theta2, iterationCount);\n",
       "    Reset(TargetReg);\n",
       "    Reset(AncilReg);\n",
       "    return Results;\n",
       "}"
      ]
     },
     {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
       "vscode": {
        "languageId": "qsharp"
       }
      },
      "outputs": [],
      "source": [
       "%%qsharp\n",
       "\n",
       "operation PerformMeasurements(theta1 : Double, theta2 : Double, n : Int) : Int {\n",
       "    let measurementCount = n + 1;\n",
       "    return QuantumInnerProduct(theta1, theta2, measurementCount);\n",
       "}"
      ]
     },
     {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
       "Set the parameters and call the `PerformMeasurements` operation to simulate it."
      ]
     },
     {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
       "vscode": {
        "languageId": "qsharp"
       }
      },
      "outputs": [],
      "source": [
       "%%qsharp\n",
       "\n",
       "PerformMeasurements(PI() / 7.0, PI() / 5.0, 4)"
      ]
     },
     {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
       "Submit the job to the target. The circuit is approximately 0.4 EHQC each shot. The number of shots specified is 128, but this can be increased to reduce the variance of the result, up to some stable distribution. The total EHQCs for a job can be viewed within the Azure portal under \"Job management\". Selecting the desired job, the cost estimation can be viewed."
      ]
     },
     {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
       "vscode": {
        "languageId": "qsharp"
       }
      },
      "outputs": [],
      "source": [
       "job = target.submit(qsharp.compile(\"PerformMeasurements(PI() / 7.0, PI() / 5.0, 4)\"), shots=128)\n",
       "job.wait_until_completed()\n",
       "job.get_results()"
      ]
     }
    ],
    "metadata": {
     "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
     },
     "language_info": {
      "codemirror_mode": {
       "name": "ipython",
       "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.11.0"
     }
    },
    "nbformat": 4,
    "nbformat_minor": 2
   }
   