{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a56dec3d",
   "metadata": {},
   "source": [
    "# Variational Quantum Algorithms in Q\\#"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "96818c17",
   "metadata": {},
   "source": [
    "## Abstract"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c55ddd55",
   "metadata": {},
   "source": [
    "In this sample, we will explore how the rich classical control provided by Q# can be used to efficiently write out variational quantum algorithms. In particular, we'll focus on the example of the _variational quantum eigensolver_, also known as VQE."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4fc6c47b",
   "metadata": {},
   "source": [
    "## Preamble\n",
    "\n",
    "$\n",
    "    \\renewcommand{\\ket}[1]{\\left|#1\\right\\rangle}\n",
    "$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "11fa63fc",
   "metadata": {},
   "outputs": [],
   "source": [
    "from itertools import product\n",
    "import qutip as qt\n",
    "from qdk import qsharp\n",
    "\n",
    "qsharp.init(project_root=\"SPSA\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9a80f41f",
   "metadata": {},
   "source": [
    "For Q# code embedded in this notebook, it will be helpful to open a few namespaces before we proceed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "373cb753",
   "metadata": {
    "vscode": {
     "languageId": "qsharp"
    }
   },
   "outputs": [],
   "source": [
    "%%qsharp\n",
    "import Std.Arrays.*;\n",
    "import Std.Convert.*;\n",
    "import Std.Diagnostics.*;\n",
    "import Std.Random.*;\n",
    "import Std.Math.*;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dce084a9",
   "metadata": {},
   "source": [
    "## Introducing the Variational Quantum Eigensolver"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "983fbce7",
   "metadata": {},
   "source": [
    "In variational quantum algorithms, rather than performing all of our computation on the quantum device itself, we use a quantum program to estimate some quantity, then use a classical optimizer to find inputs to our quantum program that minimize or maximize that quantity. That is, rather than thinking of classical computation purely as a pre-processing or post-processing step, our computation uses both classical and quantum computation together."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "150b8540",
   "metadata": {},
   "source": [
    "> **💡 TIP:** We could also consider using classical computation while qubits are still alive, rather than returning an estimate to a classical optimizer. This approach is indeed very useful, for instance in iterative phase estimation. For this sample, however, we'll focus on the variational case."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8f8c35bf",
   "metadata": {},
   "source": [
    "\n",
    "For example, consider the problem of finding the _ground state energy_ of a given operator $H$, often called a _Hamiltonian_. The ground state energy of $H$ is defined as its smallest eigenvalue $E_0$, as the Hamiltonian represents the possible energy configurations of a given physical system.\n",
    "\n",
    "As stated, even though this is a minimization, we need to know the eigenvalues of $H$ to make any progress. It's pretty straightforward to rephrase the problem in terms of arbitrary states, however. In particular, the expectation value $\\left\\langle H \\right\\rangle H = \\left\\langle \\psi | H | \\psi \\right\\rangle$ must be at least as large as $E_0$. To see this, we can expand the expectation value in terms of the eigenvectors $\\{\\ket{\\phi_i}\\}$ of $H$, such that $H\\ket{\\phi_i} = E_i$.\n",
    "\n",
    "Since the decomposition of $H$ into eigenvectors gives us a basis, we know that $\\left\\langle \\phi_i | \\phi_j \\right\\rangle$ is zero whenever $i \\ne j$, allowing us to expand $\\ket{\\psi}$ into the eigenbasis of $H$ as $\\ket{\\psi} = \\sum_i \\alpha_i \\ket{\\phi_i}$ for some complex coefficients $\\{\\alpha_i\\}$. Using this decomposition, we can expand the expectation $\\left\\langle \\psi | H | \\psi \\right\\rangle$ in terms of the eigenvalues of $H$:\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "    \\left\\langle \\psi | H | \\psi \\right\\rangle\n",
    "        & = \\sum_i |\\alpha_i|^2 \\left\\langle \\phi_i | H | \\phi_i \\right\\rangle \\\\\n",
    "        & = \\sum_i |\\alpha_i|^2 \\left\\langle \\phi_i | E_i \\phi_i \\right\\rangle \\\\\n",
    "        & = \\sum_i |\\alpha_i|^2 E_i.\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "Since each $|\\alpha_i|^2$ is nonnegative, and since $\\sum_i |\\alpha_i|^2 = 1$, the above gives us that\n",
    "$$\n",
    "    E_0 = \\min_i E_i \\le \\left\\langle \\psi | H | \\psi \\right\\rangle \\le \\max_i E_i.\n",
    "$$\n",
    "\n",
    "Thus, we can rephrase the original problem as a minimization not just over eigenstates, but of all arbitrary states,\n",
    "$$\n",
    "    E_0 = \\min_{\\ket{\\psi}} \\left\\langle \\psi | H | \\psi \\right\\rangle,\n",
    "$$\n",
    "where the minimum is achieved when $\\ket{\\psi} = \\ket{\\phi_0}$.\n",
    "\n",
    "Using this rephrasing, the _variational quantum eigensolver_ algorithm is just the variational algorithm that we get by using a classical optimizer to find $\\min_{\\ket{\\psi}} \\left\\langle \\psi | H | \\psi \\right\\rangle$. In pseudocode:\n",
    "\n",
    "```\n",
    "operation FindMinimumEigenvalue {\n",
    "    Pick an initial guess |ψ⟩.\n",
    "    until target accuracy reached {\n",
    "        Estimate E = ⟨ψ | H | ψ⟩ using a quantum operation.\n",
    "        Use a classical optimizer to pick the next |ψ⟩.\n",
    "    }\n",
    "    Return the minimum E that we found and the state |ψ⟩ that achieved that minimum.\n",
    "}\n",
    "```"
   ]
  },
  {
   "attachments": {
    "image.png": {
     "image/png": "[encoded data, 19700 characters]"
    }
   },
   "cell_type": "markdown",
   "id": "41037e71",
   "metadata": {},
   "source": [
    "![image.png](attachment:image.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "68e94097",
   "metadata": {},
   "source": [
    "To turn the above pseudocode into a real quantum program, we still need to figure out two things:\n",
    "\n",
    "- How to estimate $\\left\\langle \\psi | H | \\psi \\right\\rangle$ for a given state $\\ket{\\psi}$.\n",
    "- How to optimize over all quantum states $\\ket{\\psi}$.\n",
    "\n",
    "In the rest of this sample, we'll see how to do each of these, and how you can use Q# to write a VQE implementation."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2f7b662b",
   "metadata": {},
   "source": [
    "## Estimating expectation values of $H$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c3bb1ba6",
   "metadata": {},
   "source": [
    "Before proceeding to see how to estimate $\\left\\langle \\psi | H | \\psi \\right\\rangle$, however, it helps to consider what that expectation value _is_. To do so, let's consider a concrete example of a two-qubit Hamiltonian, using QuTiP to construct a Python object representing $H$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a992fe41",
   "metadata": {},
   "outputs": [],
   "source": [
    "H = qt.Qobj([\n",
    "    [-2, 0, 0, 3],\n",
    "    [0, 7, 0, 1],\n",
    "    [0, 0, -4, 0],\n",
    "    [3, 1, 0, -1]\n",
    "], dims=[[2, 2]] * 2)\n",
    "H"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4b09b667",
   "metadata": {},
   "source": [
    "Here, `H` is a Python object representing the 4 × 4 matrix $H$. We can use the `eigenstates` method provided by QuTiP to quickly find the eigenvalues and eigenvectors of $H$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "33551a8c",
   "metadata": {},
   "outputs": [],
   "source": [
    "H.eigenstates()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba4c191e",
   "metadata": {},
   "source": [
    "In particular, we can use the `min` function provided by Python to minimize over the eigenvalues of $H$ and find its ground state $\\ket{\\phi_0}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3199b955",
   "metadata": {},
   "outputs": [],
   "source": [
    "min_energy, ground_state = min(zip(*H.eigenstates()), key=lambda eig: eig[0])\n",
    "min_energy, ground_state"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2ce8638c",
   "metadata": {},
   "source": [
    "Here, `ground_state` represents the eigenvector $\\ket{\\phi_0}$ corresponding to the smallest eigenvalue $E_0$ of $H$. By the above argument, we would expect that $\\left\\langle \\phi_0 | H | \\phi_0 \\right\\rangle = E_0$. We can check that using QuTiP as well:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a2140642",
   "metadata": {},
   "outputs": [],
   "source": [
    "(ground_state.dag() * H * ground_state)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "853df02f",
   "metadata": {},
   "source": [
    "What's going on here? Effectively, for any operator $O$ such that $O^{\\dagger} = O$, expressions like $\\left\\langle \\psi | O | \\psi \\right\\rangle$ represent another way of thinking about quantum measurement. In particular, if we think of the eigenvalues of $O$ as labels for its corresponding various eigenvectors, then $\\left\\langle \\psi | O | \\psi \\right\\rangle$ is the average label that we get when we measure $O$ in the basis of its eigenvectors."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8537e34c",
   "metadata": {},
   "source": [
    "> **💡 TIP:** This way of thinking about quantum measurement is sometimes called the _observable framework_, with $O$ being called an _observable_. We avoid using this terminology here to avoid confusion, however, as the expectation value of $O$ cannot be observed directly, but only inferred from repeated measurements."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "817f94ec",
   "metadata": {},
   "source": [
    "For example, if $O = Z$, then the eigenstate $\\ket{0}$ is labeled by its eigenvalue $+1$, while $\\ket{1}$ is labeled by $-1$. The expectation value $\\left\\langle \\psi | Z | \\psi \\right\\rangle$ is then the probability of getting a `Zero` minus the probability of getting a `One`.\n",
    "\n",
    "More generally, finding the expectation value of an arbitrary Pauli operator for an arbitrary input state is straightforward using the `Measure` and `EstimateFrequencyA` operations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d6558958",
   "metadata": {
    "vscode": {
     "languageId": "qsharp"
    }
   },
   "outputs": [],
   "source": [
    "%%qsharp\n",
    "\n",
    "import Chemistry.JordanWigner.VQE.EstimateFrequency;\n",
    "\n",
    "operation EstimatePauliExpectation(pauli : Pauli[], preparation : (Qubit[] => Unit is Adj), nShots : Int) : Double {\n",
    "    return 2.0 * EstimateFrequency(\n",
    "        preparation,\n",
    "        Measure(pauli, _),\n",
    "        Length(pauli),\n",
    "        nShots\n",
    "    ) - 1.0;\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2d07f00c",
   "metadata": {},
   "source": [
    "Here, we've represented the state $\\ket{\\psi}$ by an operation `preparation` that prepares that state. Since each Pauli operator other than the identity operator has exactly two eigenvalues, $+1$ and $-1$, corresponding to `Zero` and `One` respectively, we can turn the estimate of the probability $p_0$ with which `Measure(pauli, _)` returns `Zero` into an expectation value $p_0 - p_1 = p_0 - (1 - p_0) = 2 p_0 - 1$.\n",
    "\n",
    "For example, doing nothing prepares the state $\\ket{0}$, so we get an expectation value of $1$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "22a23913",
   "metadata": {
    "vscode": {
     "languageId": "qsharp"
    }
   },
   "outputs": [],
   "source": [
    "%%qsharp\n",
    "\n",
    "operation EstimateExpectationOfZero() : Double {\n",
    "    return EstimatePauliExpectation(\n",
    "        [PauliZ],\n",
    "        qs => I(qs[0]),\n",
    "        100\n",
    "    );\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1135713f",
   "metadata": {},
   "outputs": [],
   "source": [
    "from qdk import code\n",
    "code.EstimateExpectationOfZero()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aef4f547",
   "metadata": {},
   "source": [
    "Similarly, using `X` to prepare $\\ket{1}$ gives us an expectation of $-1$, while using `H` to prepare $\\ket{+} = \\frac{1}{\\sqrt{2}} (\\ket{0} + \\ket{1})$ gives an expectation value of $0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8dbf0842",
   "metadata": {
    "vscode": {
     "languageId": "qsharp"
    }
   },
   "outputs": [],
   "source": [
    "%%qsharp\n",
    "\n",
    "operation EstimateExpectationOfOne() : Double {\n",
    "    return EstimatePauliExpectation(\n",
    "        [PauliZ],\n",
    "        ApplyToEachCA(X, _),\n",
    "        100\n",
    "    );\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2cab1fca",
   "metadata": {},
   "outputs": [],
   "source": [
    "code.EstimateExpectationOfOne()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "408d1da3",
   "metadata": {
    "vscode": {
     "languageId": "qsharp"
    }
   },
   "outputs": [],
   "source": [
    "%%qsharp\n",
    "\n",
    "operation EstimateExpectationOfPlus() : Double {\n",
    "    return EstimatePauliExpectation(\n",
    "        [PauliZ],\n",
    "        ApplyToEachCA(H, _),\n",
    "        100\n",
    "    );\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fc6d8273",
   "metadata": {},
   "outputs": [],
   "source": [
    "code.EstimateExpectationOfPlus()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6d023ef6",
   "metadata": {},
   "source": [
    "Note that in practice, we won't always get exactly 0 due to using a finite number of measurements."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b2d9316d",
   "metadata": {},
   "source": [
    "In any case, to recap, it's easy to use a quantum program to find $\\left\\langle \\psi | H | \\psi \\right\\rangle$ in the special case that $H$ is a multi-qubit Pauli operator. What about the more general case? The linearity of quantum mechanics saves us again! It turns out that we can expand the expectation of any other operator in terms of expectations of Pauli operators. To see how that works, suppose that $H = 2 Z - X$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ae9236c7",
   "metadata": {},
   "outputs": [],
   "source": [
    "2 * qt.sigmaz() - qt.sigmax()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8d7db2d3",
   "metadata": {},
   "source": [
    "We can expand $\\left\\langle \\psi | (2Z - X) | \\psi \\right\\rangle$ by using linearity:\n",
    "\n",
    "$$\n",
    "    \\left\\langle \\psi | (2Z - X) | \\psi \\right\\rangle = 2 \\left\\langle \\psi | Z | \\psi \\right\\rangle - \\left\\langle \\psi | X | \\psi \\right\\rangle.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "797bd006",
   "metadata": {},
   "source": [
    "Each of the terms in this expansion is something that we can estimate easily using our `EstimatePauliExpectation` operation from above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d10b40ab",
   "metadata": {
    "vscode": {
     "languageId": "qsharp"
    }
   },
   "outputs": [],
   "source": [
    "%%qsharp\n",
    "\n",
    "operation EstimateExpectation(terms : (Double, Pauli[])[], preparation : (Qubit[] => Unit is Adj), nShots : Int) : Double {\n",
    "    mutable sum = 0.0;\n",
    "    for (coefficient, pauli) in terms {\n",
    "        sum += coefficient * EstimatePauliExpectation(\n",
    "            pauli, preparation, nShots\n",
    "        );\n",
    "    }\n",
    "    return sum;\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "49bf9fcc",
   "metadata": {},
   "source": [
    "With this, we almost have everything we need to estimate expectations $\\left\\langle \\psi | H | \\psi \\right\\rangle$. We just need a way of finding the decomposition of $H$ into Pauli operators. Thankfully, QuTiP can help here as well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5496baa1",
   "metadata": {},
   "outputs": [],
   "source": [
    "def pauli_to_qobj(pauli):\n",
    "    if pauli == qsharp.Pauli.I:\n",
    "        return qt.Qobj(qt.identity(2))\n",
    "    elif pauli == qsharp.Pauli.X:\n",
    "        return qt.Qobj(qt.sigmax())\n",
    "    elif pauli == qsharp.Pauli.Y:\n",
    "        return qt.Qobj(qt.sigmay())\n",
    "    elif pauli == qsharp.Pauli.Z:\n",
    "        return qt.Qobj(qt.sigmaz())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "731556db",
   "metadata": {},
   "outputs": [],
   "source": [
    "def pauli_basis(n_qubits):\n",
    "    scale = 2 ** n_qubits\n",
    "    return [\n",
    "        tuple([P, qt.tensor(*(pauli_to_qobj(p) for p in P)) / scale])\n",
    "        for P in product([qsharp.Pauli.I, qsharp.Pauli.X, qsharp.Pauli.Y, qsharp.Pauli.Z], repeat=n_qubits)\n",
    "        ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c8e0b922",
   "metadata": {},
   "outputs": [],
   "source": [
    "def expand_in_pauli_basis(op):\n",
    "    return [\n",
    "        (coeff.real, list(label))\n",
    "        for label, coeff in [\n",
    "            tuple([label, (P.dag() * op).tr()])\n",
    "            for label, P in pauli_basis(n_qubits=len(op.dims[0]))\n",
    "        ]\n",
    "        if abs(coeff) >= 1e-10\n",
    "    ]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ff5c27bc",
   "metadata": {},
   "source": [
    "For example, we can use the two functions above to find that $H = -3 𝟙 \\otimes Z + 0.5 X \\otimes 𝟙 + 1.5 X \\otimes X - 0.5 X \\otimes Z - 1.5 Y \\otimes Y + 2.5 Z \\otimes 𝟙 - 1.5 Z \\otimes Z$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cd8e5438",
   "metadata": {},
   "outputs": [],
   "source": [
    "H_decomposition = expand_in_pauli_basis(H)\n",
    "H_decomposition"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8879b674",
   "metadata": {},
   "source": [
    "This decomposition is exactly what we need to pass as `terms` above. For example, to estimate the expectation of the $\\ket{++}$ state, $\\left\\langle ++ | H | ++ \\right\\rangle$, we can pass `terms` to `EstimateExpectation`. In doing so, we'll use the name \"energy\" for our new operation, pointing to that expectation values of Hamiltonian operators $H$ represent the average energy of a system given the quantum state of that system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fad87b56",
   "metadata": {
    "vscode": {
     "languageId": "qsharp"
    }
   },
   "outputs": [],
   "source": [
    "%%qsharp\n",
    "\n",
    "operation EstimateEnergyOfPlus(terms : (Double, Pauli[])[]) : Double {\n",
    "    return EstimateExpectation(\n",
    "        terms,\n",
    "        ApplyToEachCA(H, _),\n",
    "        100\n",
    "    );\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1f0e22df",
   "metadata": {},
   "outputs": [],
   "source": [
    "code.EstimateEnergyOfPlus(H_decomposition)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "29a48a9f",
   "metadata": {},
   "source": [
    "This is pretty far from the minimum $E_0$ we found from the eigenvalue decomposition above, so in the next part we'll see one more trick we can use to write our VQE implementation."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9bfad16a",
   "metadata": {},
   "source": [
    "## Preparing ansatz states"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3ae1ddfe",
   "metadata": {},
   "source": [
    "Recall that our goal is to use a classical optimizer to find the minimum energy $E_0$ of some Hamiltonian operator $H$:\n",
    "$$\n",
    "    E_0 = \\min_{\\ket{\\psi}} \\left\\langle \\psi | H | \\psi \\right\\rangle.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fe83f982",
   "metadata": {},
   "source": [
    "In practice, we can't reasonably optimize over all possible quantum states $\\ket{\\psi}$, so that the above optimization problem is intractable for all but the smallest systems. Instead, we note that we can find an upper-bound for $E_0$ by solving a simpler problem. Suppose that there is a set $A$ of states that are easier to prepare. Then,\n",
    "$$\n",
    "    E_0 = \\min_{\\ket{\\psi}} \\left\\langle \\psi | H | \\psi \\right\\rangle \\le \\min_{\\ket{\\psi} \\in A} \\left\\langle \\psi | H | \\psi \\right\\rangle.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3926dbbf",
   "metadata": {},
   "source": [
    "If we choose $A$ carefully so that $\\ket{\\phi_0} \\in A$, this inequality will be tight. More often, however, $\\ket{\\phi_0}$ won't actually be in $A$, but will be close to some other state in $A$, giving us a reasonable approximation to $E_0$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c6bce67e",
   "metadata": {},
   "source": [
    "We say that $A$ is our _ansatz_ for running VQE; in this way, $A$ plays a similar role to a model in an ML training loop, representing what we believe to be a reasonable set of guesses for $\\ket{\\phi_0}$. One can get pretty involved with their choice of ansatz, but for this sample, we'll choose a really simple one, and set $A$ to be those states that can be prepared by a small number of Pauli rotations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "950a5184",
   "metadata": {
    "vscode": {
     "languageId": "qsharp"
    }
   },
   "outputs": [],
   "source": [
    "%%qsharp\n",
    "\n",
    "operation PrepareAnsatz(axes : Pauli[][], angles : Double[], register : Qubit[]) : Unit is Adj {\n",
    "    for i in IndexRange(angles) {\n",
    "        Exp(axes[i], angles[i], register);\n",
    "    }\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b377850f",
   "metadata": {},
   "source": [
    "Our minimization $\\min_{\\ket{\\psi} \\in A} \\left\\langle \\psi | H | \\psi \\right\\rangle$ is now a minimization over `angles` for some fixed list of Pauli rotation axes. Using `DumpMachine` and `DumpRegister`, we can visualize how this ansatz works for a few different choices of parameters `angles`, convincing ourselves that $A$ contains enough interesting states to find $E_0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ec65a772",
   "metadata": {
    "vscode": {
     "languageId": "qsharp"
    }
   },
   "outputs": [],
   "source": [
    "%%qsharp\n",
    "\n",
    "operation DumpAnsatz(axes : Pauli[][], angles : Double[]) : Unit {\n",
    "    use register = Qubit[Length(angles)];\n",
    "    within {\n",
    "        PrepareAnsatz(axes, angles, register);\n",
    "    } apply {\n",
    "        DumpRegister(register);\n",
    "    }\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b99c1d78",
   "metadata": {},
   "outputs": [],
   "source": [
    "ansatz_axes = [[qsharp.Pauli.Z, qsharp.Pauli.Z], [qsharp.Pauli.X, qsharp.Pauli.Y]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7df99a83",
   "metadata": {},
   "outputs": [],
   "source": [
    "code.DumpAnsatz(ansatz_axes, [1.2, 1.9])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fc79ef91",
   "metadata": {},
   "outputs": [],
   "source": [
    "code.DumpAnsatz(ansatz_axes, [1.2, -0.7])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6d2f1e81",
   "metadata": {},
   "source": [
    "At this point, it's helpful to write some new operations that directly estimate the energy of a state given a parameterization of our ansatz, rather than an opaque operation that prepares the state."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "11583999",
   "metadata": {
    "vscode": {
     "languageId": "qsharp"
    }
   },
   "outputs": [],
   "source": [
    "%%qsharp\n",
    "\n",
    "operation EstimateEnergyAtAnsatz(terms : (Double, Pauli[])[], axes : Pauli[][], angles : Double[], nShots : Int) : Double {\n",
    "    return EstimateExpectation(\n",
    "        terms,\n",
    "        PrepareAnsatz(axes, angles, _),\n",
    "        nShots\n",
    "    );\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "418ace1a",
   "metadata": {},
   "outputs": [],
   "source": [
    "code.EstimateEnergyAtAnsatz(\n",
    "    H_decomposition,\n",
    "    ansatz_axes,\n",
    "    [1.2, 1.9],\n",
    "    1000\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6f07737a",
   "metadata": {},
   "source": [
    "## Running VQE in Q\\#"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b0970ee4",
   "metadata": {},
   "source": [
    "Using this new operation, we have something that looks much more like classical optimization problems that we may be used to! Indeed, we can simply use our favorite optimization algorithm to pick different values for `angles` until we find a good approximation for $E_0$. In this sample, we'll use the _SPSA algorithm_ to optimize `angles`, as this algorithm works especially well for objective functions that have some amount of noise, such as that added by taking a finite number of shots above.\n",
    "\n",
    "We won't go through the details of SPSA here, but if you're interested to learn more, check out [`SPSA.qs`](./SPSA/src/SPSA.qs) in this folder to see how we implemented SPSA in Q#."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7f3049df",
   "metadata": {
    "vscode": {
     "languageId": "qsharp"
    }
   },
   "outputs": [],
   "source": [
    "%%qsharp\n",
    "\n",
    "import SPSA.*;\n",
    "\n",
    "operation FindMinimumEnergy(terms : (Double, Pauli[])[], axes : Pauli[][], initialGuess : Double[], nShots : Int) : (Double[], Double) {\n",
    "    let oracle = EstimateEnergyAtAnsatz(terms, axes, _, nShots);\n",
    "    let options = DefaultSpsaOptions();\n",
    "    return FindMinimumWithSpsa(\n",
    "        oracle,\n",
    "        initialGuess,\n",
    "        options\n",
    "    );\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "337c5318",
   "metadata": {},
   "outputs": [],
   "source": [
    "code.FindMinimumEnergy(H_decomposition, ansatz_axes, [1.2, 1.9], 1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d6276b06",
   "metadata": {},
   "outputs": [],
   "source": [
    "min_energy"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "464ea5c5cae2889b0c59ecc40909802f0909ef10120e6598fbb774b74247ba9f"
  },
  "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.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
