{
  "cells": [
    {
      "cell_type": "markdown",
      "id": "e0cf0747",
      "metadata": {},
      "source": [
        "---\n",
        "title: Variational quantum algorithms\n",
        "description: This tutorial provides an overview of a hybrid quantum-classical algorithm, VQE, and the QAOA\n",
        "---\n",
        "\n",
        "# Quantum algorithms: Variational quantum algorithms\n",
        "\n",
        "<Admonition type=\"note\">\n",
        "  Takashi Imamichi (24 May 2024)\n",
        "\n",
        "  [Download the pdf](https://ibm.ent.box.com/s/blnffu0pd7yzxarq3zc3w0jv90365ny2) of the original lecture. Note that some code snippets might become deprecated since these are static images.\n",
        "\n",
        "  *Approximate QPU time to run this experiment is 9 minutes (tested on an Eagle processor).*\n",
        "\n",
        "  (this notebook might not evaluate in the time allowed on the Open Plan. Please use quantum computing resources wisely.)\n",
        "</Admonition>\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "59fbaadc",
      "metadata": {},
      "source": [
        "## 1. Introduction\n",
        "\n",
        "This tutorial provides an overview of a hybrid quantum-classical algorithm, specifically focusing on the variational quantum eigensolver (VQE) and the quantum approximate optimization algorithm (QAOA). The primary objective of these algorithms is to tackle optimization problems by employing quantum circuits with parameterized quantum gates.\n",
        "\n",
        "Despite the advancements in quantum computing, the presence of noise in current quantum devices makes it challenging to extract meaningful results from deep quantum circuits. To overcome this challenge, VQE and QAOA adopt a hybrid quantum-classical approach, which involves iteratively executing relatively short quantum circuits using quantum computation and optimizing the parameters of the target parametrized quantum circuits using classical computation.\n",
        "\n",
        "QAOA has the potential to provide the optimal solutions to the target problems at a utility scale, thanks to the application of various error mitigation and suppression techniques. VQE has many applications (like quantum chemistry) in which it is less scalable. But there have been a number of eigenvalue-related approaches that have emerged to complement and augment VQE, including Krylov subspace diagonalization and sampling-based quantum diagonalization (SQD). Understanding VQE is an important first step in understanding the wide range of classical-quantum hybrid algorithms that have emerged.\n",
        "\n",
        "This module describes the fundamental concepts and implementation of VQE and QAOA. Further tutorials will explore advanced topics and techniques for scaling up these algorithms.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "d3f4500f",
      "metadata": {},
      "source": [
        "You require the following library in your environment to run this notebook.\n",
        "If you have not installed it yet, you can install it by un-commenting and running the following cell.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "2a8e650f",
      "metadata": {},
      "outputs": [],
      "source": [
        "# % pip install 'qiskit[visualization]' qiskit-ibm-runtime"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "be601762",
      "metadata": {},
      "source": [
        "## 2. Computing the minimum eigenvalue of a simple Hamiltonian\n",
        "\n",
        "We will start by applying VQE to a very simple case, to see how it works. We will compute the minimum eigenvalue of Pauli $Z$ matrix with VQE. We will start by importing some general packages.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "id": "e8f398b2",
      "metadata": {},
      "outputs": [],
      "source": [
        "import numpy as np\n",
        "from qiskit.circuit import ParameterVector, QuantumCircuit\n",
        "from qiskit.primitives import StatevectorEstimator, StatevectorSampler\n",
        "from qiskit.quantum_info import SparsePauliOp\n",
        "from scipy.optimize import minimize"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "d4395012",
      "metadata": {},
      "source": [
        "We now define the operator of interest and view it in matrix form.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "id": "a4f374a2",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[ 1.+0.j,  0.+0.j],\n",
              "       [ 0.+0.j, -1.+0.j]])"
            ]
          },
          "execution_count": 2,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "op = SparsePauliOp(\"Z\")\n",
        "op.to_matrix()"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "0a001e71",
      "metadata": {},
      "source": [
        "It is easy to obtain the eigenvalues classically, so we can check our work. This might become difficult as we scale toward utility. Here we use numpy.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "id": "c9188662",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Eigenvalues: [-1.  1.]\n"
          ]
        }
      ],
      "source": [
        "# compute eigenvalues with numpy\n",
        "result = np.linalg.eigh(op.to_matrix())\n",
        "print(\"Eigenvalues:\", result.eigenvalues)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "8f9af345",
      "metadata": {},
      "source": [
        "To obtain eigenvalues using a variational quantum algorithm, we construct a circuit with gates that take variational parameters:\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 4,
      "id": "99d5d36b",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Image src=\"/learning/images/courses/utility-scale-quantum-computing/variational-quantum-algorithms/extracted-outputs/99d5d36b-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "execution_count": 4,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "# define a variational form\n",
        "param = ParameterVector(\"a\", 3)\n",
        "qc = QuantumCircuit(1, 1)\n",
        "qc.u(param[0], param[1], param[2], 0)\n",
        "qc_estimator = qc.copy()\n",
        "qc.measure(0, 0)\n",
        "qc.draw(\"mpl\")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "9607ab51",
      "metadata": {},
      "source": [
        "If we want to estimate the expectation value of an operator (like $Z$), we should use Estimator. If we want to look at the states of the system, we use Sampler.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "id": "df8ed8d3",
      "metadata": {},
      "outputs": [],
      "source": [
        "sampler = StatevectorSampler()\n",
        "estimator = StatevectorEstimator()"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "880bccfb",
      "metadata": {},
      "source": [
        "We can compute counts of bitstrings 0 and 1 with random parameter values `[1, 2, 3]` using Sampler.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "id": "5301ee71",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "{'0': 783, '1': 241}"
            ]
          },
          "execution_count": 6,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "# compute counts of bitstrings with random parameter values by Sampler\n",
        "result = sampler.run([(qc, [1, 2, 3])]).result()\n",
        "counts = result[0].data.c.get_counts()\n",
        "counts"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "30ef9b8f",
      "metadata": {},
      "source": [
        "We know that we can compute the expectation value of Z by $\\langle Z \\rangle = p_0 - p_1$ with probabilities $\\{0: p_0, 1: p_1\\}$.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 7,
      "id": "0f220c8d",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "0.529296875"
            ]
          },
          "execution_count": 7,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "# compute the expectation value of Z based on the counts\n",
        "(counts.get(\"0\", 0) - counts.get(\"1\", 0)) / sum(counts.values())"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "695ee4bd",
      "metadata": {},
      "source": [
        "This circuit worked, but the parameter values chosen did not correspond to a very low-energy (or low-eigenvalue) state. The eigenvalue obtained is quite a bit higher than the minimum. The result is similar when using estimator.\n",
        "\n",
        "Note that Estimator takes quantum circuits without measurements.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 8,
      "id": "c9e05530",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array(0.54030231)"
            ]
          },
          "execution_count": 8,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "result = estimator.run([(qc_estimator, op, [1, 2, 3])]).result()\n",
        "result[0].data.evs"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "b538e8b5",
      "metadata": {},
      "source": [
        "We will need to search through parameters and find those that yield the lowest eigenvalue.\n",
        "We make a function to receive parameter values of the variational form and return the expectation value $\\langle Z \\rangle$.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 9,
      "id": "97fead46",
      "metadata": {},
      "outputs": [],
      "source": [
        "# define a cost function to look for the minimum eigenvalue of Z\n",
        "def cost(x):\n",
        "    result = sampler.run([(qc, x)]).result()\n",
        "    counts = result[0].data.c.get_counts()\n",
        "    expval = (counts.get(\"0\", 0) - counts.get(\"1\", 0)) / sum(counts.values())\n",
        "    # the following line shows the trajectory of the optimization\n",
        "    print(expval, counts)\n",
        "    return expval"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "a97eaeba",
      "metadata": {},
      "source": [
        "Let's apply SciPy's `minimize` function to find the minimum eigenvalue of Z.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 10,
      "id": "44f56300",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "1.0 {'0': 1024}\n",
            "0.494140625 {'0': 765, '1': 259}\n",
            "0.466796875 {'0': 751, '1': 273}\n",
            "0.564453125 {'0': 801, '1': 223}\n",
            "-0.4296875 {'1': 732, '0': 292}\n",
            "-0.984375 {'1': 1016, '0': 8}\n",
            "-0.8984375 {'1': 972, '0': 52}\n",
            "-0.990234375 {'1': 1019, '0': 5}\n",
            "-0.892578125 {'1': 969, '0': 55}\n",
            "-0.986328125 {'1': 1017, '0': 7}\n",
            "-0.861328125 {'1': 953, '0': 71}\n",
            "-1.0 {'1': 1024}\n",
            "-0.982421875 {'1': 1015, '0': 9}\n",
            "-0.99609375 {'1': 1022, '0': 2}\n",
            "-0.986328125 {'1': 1017, '0': 7}\n",
            "-1.0 {'1': 1024}\n",
            "-0.990234375 {'1': 1019, '0': 5}\n",
            "-0.998046875 {'1': 1023, '0': 1}\n",
            "-0.99609375 {'1': 1022, '0': 2}\n",
            "-1.0 {'1': 1024}\n",
            "-1.0 {'1': 1024}\n",
            "-1.0 {'1': 1024}\n",
            "-1.0 {'1': 1024}\n",
            "-0.998046875 {'1': 1023, '0': 1}\n",
            "-1.0 {'1': 1024}\n",
            "-1.0 {'1': 1024}\n",
            "-0.998046875 {'1': 1023, '0': 1}\n",
            "-0.998046875 {'1': 1023, '0': 1}\n",
            "-0.998046875 {'1': 1023, '0': 1}\n",
            "-1.0 {'1': 1024}\n",
            "-0.99609375 {'1': 1022, '0': 2}\n",
            "-1.0 {'1': 1024}\n",
            "-0.99609375 {'1': 1022, '0': 2}\n",
            "-0.998046875 {'1': 1023, '0': 1}\n",
            "-0.998046875 {'1': 1023, '0': 1}\n",
            "-0.99609375 {'1': 1022, '0': 2}\n",
            "-0.998046875 {'1': 1023, '0': 1}\n",
            "-1.0 {'1': 1024}\n",
            "-0.998046875 {'1': 1023, '0': 1}\n",
            "-0.998046875 {'1': 1023, '0': 1}\n",
            "-0.99609375 {'1': 1022, '0': 2}\n",
            "-1.0 {'1': 1024}\n",
            "-0.998046875 {'1': 1023, '0': 1}\n",
            "-1.0 {'1': 1024}\n",
            "-0.998046875 {'1': 1023, '0': 1}\n",
            "-0.998046875 {'1': 1023, '0': 1}\n",
            "-1.0 {'1': 1024}\n",
            "-0.998046875 {'1': 1023, '0': 1}\n",
            "-0.998046875 {'1': 1023, '0': 1}\n",
            "-1.0 {'1': 1024}\n",
            "-1.0 {'1': 1024}\n",
            "-1.0 {'1': 1024}\n",
            "-1.0 {'1': 1024}\n",
            "-1.0 {'1': 1024}\n",
            "-1.0 {'1': 1024}\n",
            "-0.998046875 {'1': 1023, '0': 1}\n",
            "-0.994140625 {'1': 1021, '0': 3}\n",
            "-1.0 {'1': 1024}\n",
            "-1.0 {'1': 1024}\n",
            "-1.0 {'1': 1024}\n",
            "-1.0 {'1': 1024}\n",
            "-1.0 {'1': 1024}\n",
            "-1.0 {'1': 1024}\n"
          ]
        },
        {
          "data": {
            "text/plain": [
              " message: Optimization terminated successfully.\n",
              " success: True\n",
              "  status: 1\n",
              "     fun: -1.0\n",
              "       x: [ 3.182e+00  1.338e+00  1.664e-01]\n",
              "    nfev: 63\n",
              "   maxcv: 0.0"
            ]
          },
          "execution_count": 10,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "# minimize the cost function with scipy's minimize\n",
        "min_result = minimize(cost, [0, 0, 0], method=\"COBYLA\", tol=1e-8)\n",
        "min_result"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 11,
      "id": "4198224e",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "{'0': 1, '1': 1023}"
            ]
          },
          "execution_count": 11,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "# check counts of bitstrings with the optimal parameters\n",
        "result = sampler.run([(qc, min_result.x)]).result()\n",
        "result[0].data.c.get_counts()"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "d4e77e29",
      "metadata": {},
      "source": [
        "### 2.1 Exercise\n",
        "\n",
        "Compute the minimum eigenvalue of $Z \\otimes Z$ with VQE.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 12,
      "id": "86a2d76d",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "SparsePauliOp(['ZZ'],\n",
            "              coeffs=[1.+0.j])\n",
            "[[ 1.+0.j  0.+0.j  0.+0.j  0.+0.j]\n",
            " [ 0.+0.j -1.+0.j  0.+0.j  0.+0.j]\n",
            " [ 0.+0.j  0.+0.j -1.+0.j  0.+0.j]\n",
            " [ 0.+0.j  0.+0.j  0.+0.j  1.+0.j]]\n"
          ]
        }
      ],
      "source": [
        "z2 = SparsePauliOp(\"ZZ\")\n",
        "print(z2)\n",
        "print(z2.to_matrix())"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 15,
      "id": "76092648-3b28-442e-8319-8a8416e0c9d5",
      "metadata": {},
      "outputs": [],
      "source": [
        "# compute eigenvalues with numpy"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 16,
      "id": "f92a23b5",
      "metadata": {},
      "outputs": [],
      "source": [
        "# define a variational form\n",
        "# qc = ..."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 17,
      "id": "496572fb",
      "metadata": {
        "scrolled": true
      },
      "outputs": [],
      "source": [
        "# compute counts of bitstrings with a random parameter values by Sampler\n",
        "# result = sampler.run(...)\n",
        "# result"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 18,
      "id": "b4cad5d0",
      "metadata": {},
      "outputs": [],
      "source": [
        "# compute the expectation value of ZZ based on the counts"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 19,
      "id": "5ce555fe",
      "metadata": {},
      "outputs": [],
      "source": [
        "# verify the expectation value of ZZ with Estimator"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 20,
      "id": "c210c435",
      "metadata": {},
      "outputs": [],
      "source": [
        "# define a cost function to look for the minimum eigenvalue of ZZ\n",
        "# def cost(x):\n",
        "#    expval = ...\n",
        "#    return expval"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 21,
      "id": "a75f3303",
      "metadata": {},
      "outputs": [],
      "source": [
        "# minimize the cost function with scipy's minimize\n",
        "# min_result = minimize(cost, [...], method=\"COBYLA\", tol=1e-8)\n",
        "# min_result"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 22,
      "id": "8c476943",
      "metadata": {},
      "outputs": [],
      "source": [
        "# check counts of bitstrings with the optimal parameter values\n",
        "# result = sampler.run(qc, min_result.x).result()\n",
        "# result"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "b82ba33a",
      "metadata": {},
      "source": [
        "#### Solutions of the exercise\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "7f2e0a4e",
      "metadata": {},
      "source": [
        "We define the operator of interest and view it in matrix form.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 13,
      "id": "a1f9d371",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "SparsePauliOp(['ZZ'],\n",
            "              coeffs=[1.+0.j])\n",
            "[[ 1.+0.j  0.+0.j  0.+0.j  0.+0.j]\n",
            " [ 0.+0.j -1.+0.j  0.+0.j  0.+0.j]\n",
            " [ 0.+0.j  0.+0.j -1.+0.j  0.+0.j]\n",
            " [ 0.+0.j  0.+0.j  0.+0.j  1.+0.j]]\n"
          ]
        }
      ],
      "source": [
        "z2 = SparsePauliOp(\"ZZ\")\n",
        "print(z2)\n",
        "print(z2.to_matrix())"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "b123e9d9",
      "metadata": {},
      "source": [
        "To obtain eigenvalues using a variational quantum algorithm, we construct a circuit with gates that take variational parameters:\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 14,
      "id": "7d5e894a",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Image src=\"/learning/images/courses/utility-scale-quantum-computing/variational-quantum-algorithms/extracted-outputs/7d5e894a-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "execution_count": 14,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "# define a variational form\n",
        "param = ParameterVector(\"a\", 6)\n",
        "qc = QuantumCircuit(2, 2)\n",
        "qc.u(param[0], param[1], param[2], 0)\n",
        "qc.u(param[3], param[4], param[5], 1)\n",
        "qc_estimator = qc.copy()\n",
        "qc.measure([0, 1], [0, 1])\n",
        "qc.draw(\"mpl\")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "f1d86563",
      "metadata": {},
      "source": [
        "If we want to estimate the expectation value of an operator (like $Z \\otimes Z$), we would use Estimator. If we want to look at the states of the system, we use Sampler.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 15,
      "id": "e3bda5ff",
      "metadata": {},
      "outputs": [],
      "source": [
        "sampler = StatevectorSampler()\n",
        "estimator = StatevectorEstimator()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 16,
      "id": "aeac09f7",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "{'10': 661, '11': 203, '01': 47, '00': 113}"
            ]
          },
          "execution_count": 16,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "# compute counts of bitstrings with random parameter values by Sampler\n",
        "result = sampler.run([(qc, [1, 2, 3, 4, 5, 6])]).result()\n",
        "counts = result[0].data.c.get_counts()\n",
        "counts"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 17,
      "id": "c96acbdf",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "-0.3828125"
            ]
          },
          "execution_count": 17,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "# compute the expectation value of ZZ based on the counts\n",
        "(\n",
        "    counts.get(\"00\", 0)\n",
        "    - counts.get(\"01\", 0)\n",
        "    - counts.get(\"10\", 0)\n",
        "    + counts.get(\"11\", 0)\n",
        ") / sum(counts.values())"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "d2cc2ed8",
      "metadata": {},
      "source": [
        "This circuit worked, but the parameter values chosen did not correspond to a very low-energy (or low-eigenvalue) state. The eigenvalue obtained is quite a bit higher than the minimum. The result is similar when using estimator.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 18,
      "id": "474b84c5",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array(-0.35316516)"
            ]
          },
          "execution_count": 18,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "# verify the expectation value of ZZ with Estimator\n",
        "result = estimator.run([(qc_estimator, z2, [1, 2, 3, 4, 5, 6])]).result()\n",
        "result[0].data.evs"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "25f4a655",
      "metadata": {},
      "source": [
        "We will need to search through parameters and find those that yield the lowest eigenvalue.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 19,
      "id": "f29322f3",
      "metadata": {},
      "outputs": [],
      "source": [
        "# define a cost function to look for the minimum eigenvalue of ZZ\n",
        "def cost(x):\n",
        "    result = sampler.run([(qc, x)]).result()\n",
        "    counts = result[0].data.c.get_counts()\n",
        "    expval = (\n",
        "        counts.get(\"00\", 0)\n",
        "        - counts.get(\"01\", 0)\n",
        "        - counts.get(\"10\", 0)\n",
        "        + counts.get(\"11\", 0)\n",
        "    ) / sum(counts.values())\n",
        "    print(expval, counts)\n",
        "    return expval"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 20,
      "id": "12d8ba03",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "1.0 {'00': 1024}\n",
            "0.578125 {'00': 808, '01': 216}\n",
            "0.5234375 {'00': 780, '01': 244}\n",
            "0.548828125 {'00': 793, '01': 231}\n",
            "0.3515625 {'00': 637, '10': 164, '11': 55, '01': 168}\n",
            "0.3359375 {'00': 638, '11': 46, '10': 174, '01': 166}\n",
            "0.283203125 {'00': 602, '10': 181, '01': 186, '11': 55}\n",
            "-0.087890625 {'01': 414, '00': 184, '10': 143, '11': 283}\n",
            "0.236328125 {'10': 27, '11': 623, '01': 364, '00': 10}\n",
            "-0.0625 {'11': 261, '01': 403, '00': 219, '10': 141}\n",
            "0.248046875 {'01': 366, '11': 628, '00': 11, '10': 19}\n",
            "-0.0625 {'10': 145, '11': 254, '01': 399, '00': 226}\n",
            "0.228515625 {'01': 373, '11': 609, '00': 20, '10': 22}\n",
            "0.0546875 {'11': 376, '10': 273, '01': 211, '00': 164}\n",
            "-0.447265625 {'01': 731, '10': 10, '11': 267, '00': 16}\n",
            "-0.71484375 {'01': 871, '11': 99, '00': 47, '10': 7}\n",
            "-0.46484375 {'01': 741, '00': 253, '10': 9, '11': 21}\n",
            "-0.87890625 {'01': 962, '00': 39, '11': 23}\n",
            "-0.640625 {'00': 176, '01': 837, '11': 8, '10': 3}\n",
            "-0.88671875 {'01': 966, '00': 41, '11': 17}\n",
            "-0.994140625 {'01': 1021, '11': 3}\n",
            "-0.91796875 {'01': 982, '11': 35, '00': 7}\n",
            "-0.994140625 {'01': 1021, '11': 2, '00': 1}\n",
            "-0.939453125 {'01': 993, '00': 31}\n",
            "-0.990234375 {'01': 1019, '11': 5}\n",
            "-0.90234375 {'01': 974, '00': 21, '11': 29}\n",
            "-0.98046875 {'01': 1014, '11': 10}\n",
            "-0.994140625 {'01': 1021, '00': 3}\n",
            "-0.990234375 {'01': 1019, '11': 4, '00': 1}\n",
            "-0.98828125 {'01': 1018, '11': 6}\n",
            "-0.990234375 {'01': 1019, '11': 4, '00': 1}\n",
            "-0.994140625 {'01': 1021, '11': 2, '00': 1}\n",
            "-0.99609375 {'01': 1022, '11': 2}\n",
            "-0.998046875 {'01': 1023, '00': 1}\n",
            "-0.99609375 {'01': 1022, '00': 2}\n",
            "-1.0 {'01': 1024}\n",
            "-1.0 {'01': 1024}\n",
            "-1.0 {'01': 1024}\n",
            "-0.998046875 {'01': 1023, '11': 1}\n",
            "-1.0 {'01': 1024}\n",
            "-1.0 {'01': 1024}\n",
            "-1.0 {'01': 1024}\n",
            "-1.0 {'01': 1024}\n",
            "-1.0 {'01': 1024}\n",
            "-1.0 {'01': 1024}\n",
            "-1.0 {'01': 1024}\n",
            "-0.998046875 {'01': 1023, '00': 1}\n",
            "-0.998046875 {'01': 1023, '11': 1}\n",
            "-0.998046875 {'01': 1023, '00': 1}\n",
            "-1.0 {'01': 1024}\n",
            "-1.0 {'01': 1024}\n",
            "-1.0 {'01': 1024}\n",
            "-1.0 {'01': 1024}\n",
            "-1.0 {'01': 1024}\n",
            "-1.0 {'01': 1024}\n",
            "-0.998046875 {'01': 1023, '11': 1}\n",
            "-0.998046875 {'01': 1023, '11': 1}\n",
            "-1.0 {'01': 1024}\n",
            "-1.0 {'01': 1024}\n",
            "-0.998046875 {'01': 1023, '11': 1}\n",
            "-0.998046875 {'01': 1023, '11': 1}\n",
            "-0.998046875 {'01': 1023, '00': 1}\n",
            "-1.0 {'01': 1024}\n",
            "-1.0 {'01': 1024}\n",
            "-0.998046875 {'01': 1023, '00': 1}\n",
            "-1.0 {'01': 1024}\n",
            "-1.0 {'01': 1024}\n",
            "-1.0 {'01': 1024}\n",
            "-1.0 {'01': 1024}\n",
            "-0.998046875 {'01': 1023, '11': 1}\n",
            "-0.998046875 {'01': 1023, '11': 1}\n",
            "-1.0 {'01': 1024}\n",
            "-0.998046875 {'01': 1023, '11': 1}\n",
            "-1.0 {'01': 1024}\n",
            "-1.0 {'01': 1024}\n",
            "-1.0 {'01': 1024}\n",
            "-0.998046875 {'01': 1023, '11': 1}\n",
            "-0.998046875 {'01': 1023, '11': 1}\n",
            "-1.0 {'01': 1024}\n",
            "-1.0 {'01': 1024}\n",
            "-0.998046875 {'01': 1023, '11': 1}\n",
            "-0.998046875 {'01': 1023, '11': 1}\n",
            "-0.998046875 {'01': 1023, '00': 1}\n",
            "-1.0 {'01': 1024}\n",
            "-1.0 {'01': 1024}\n",
            "-1.0 {'01': 1024}\n",
            "-0.998046875 {'01': 1023, '11': 1}\n",
            "-1.0 {'01': 1024}\n",
            "-0.99609375 {'01': 1022, '00': 1, '11': 1}\n",
            "-0.998046875 {'01': 1023, '11': 1}\n",
            "-0.998046875 {'01': 1023, '00': 1}\n",
            "-0.998046875 {'01': 1023, '11': 1}\n",
            "-1.0 {'01': 1024}\n",
            "-0.99609375 {'01': 1022, '11': 1, '00': 1}\n",
            "-1.0 {'01': 1024}\n",
            "-0.998046875 {'01': 1023, '00': 1}\n",
            "-0.994140625 {'01': 1021, '00': 3}\n",
            "-0.998046875 {'01': 1023, '00': 1}\n",
            "-0.99609375 {'01': 1022, '11': 2}\n",
            "-1.0 {'01': 1024}\n",
            "-1.0 {'01': 1024}\n",
            "-0.998046875 {'01': 1023, '11': 1}\n",
            "-1.0 {'01': 1024}\n",
            "-1.0 {'01': 1024}\n",
            "-1.0 {'01': 1024}\n",
            "-1.0 {'01': 1024}\n",
            "-1.0 {'01': 1024}\n",
            "-1.0 {'01': 1024}\n",
            "-0.998046875 {'01': 1023, '11': 1}\n",
            "-1.0 {'01': 1024}\n",
            "-1.0 {'01': 1024}\n",
            "-1.0 {'01': 1024}\n",
            "-0.998046875 {'01': 1023, '11': 1}\n",
            "-0.998046875 {'01': 1023, '11': 1}\n",
            "-1.0 {'01': 1024}\n",
            "-0.998046875 {'01': 1023, '00': 1}\n",
            "-1.0 {'01': 1024}\n",
            "-1.0 {'01': 1024}\n",
            "-1.0 {'01': 1024}\n",
            "-1.0 {'01': 1024}\n",
            "-1.0 {'01': 1024}\n",
            "-1.0 {'01': 1024}\n",
            "-0.998046875 {'01': 1023, '11': 1}\n",
            "-0.998046875 {'01': 1023, '11': 1}\n",
            "-0.998046875 {'01': 1023, '11': 1}\n",
            "-0.99609375 {'01': 1022, '11': 2}\n",
            "-1.0 {'01': 1024}\n",
            "-0.998046875 {'01': 1023, '11': 1}\n"
          ]
        },
        {
          "data": {
            "text/plain": [
              " message: Optimization terminated successfully.\n",
              " success: True\n",
              "  status: 1\n",
              "     fun: -0.998046875\n",
              "       x: [ 3.167e+00  6.940e-01  1.033e+00 -2.894e-02  8.933e-01\n",
              "            1.885e+00]\n",
              "    nfev: 128\n",
              "   maxcv: 0.0"
            ]
          },
          "execution_count": 20,
          "metadata": {},
          "output_type": "execute_result"
        },
        {
          "data": {
            "text/plain": [
              " message: Optimization terminated successfully.\n",
              " success: True\n",
              "  status: 1\n",
              "     fun: -0.99609375\n",
              "       x: [ 3.098e+00 -5.402e-01  1.091e+00 -1.004e-02  3.615e-01\n",
              "            6.913e-01]\n",
              "    nfev: 115\n",
              "   maxcv: 0.0"
            ]
          },
          "execution_count": 30,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "# minimize the cost function with scipy's minimize\n",
        "min_result = minimize(cost, [0, 0, 0, 0, 0, 0], method=\"COBYLA\", tol=1e-8)\n",
        "min_result"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "d3e2de8d",
      "metadata": {},
      "source": [
        "We obtained an eigenvalue extremely close to the minimum given to us from numpy.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 21,
      "id": "d34a8544",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "{'01': 1024}"
            ]
          },
          "execution_count": 21,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "# check counts of bitstrings with the optimal parameters\n",
        "result = sampler.run([(qc, min_result.x)]).result()\n",
        "result[0].data.c.get_counts()"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "bc51e7bf-e582-49ba-93f8-035624d56ccf",
      "metadata": {},
      "source": [
        "## 3. Quantum Optimization with Qiskit patterns\n",
        "\n",
        "In this how-to we will learn about Qiskit patterns and quantum approximate optimization. A Qiskit pattern is an intuitive, repeatable set of steps for implementing a quantum computing workflow:\n",
        "\n"
      ]
    },
    {
      "attachments": {},
      "cell_type": "markdown",
      "id": "636ed1de-fc34-4cdd-9398-6fbd7c7fc9c6",
      "metadata": {},
      "source": [
        "![\"Qiskit function\"](https://eu-de.quantum.cloud.ibm.com/learning/images/courses/utility-scale-quantum-computing/variational-quantum-algorithms/qiskit-function.avif)\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "8d3218ca-ce9e-40b4-a041-1b1d09bac8f5",
      "metadata": {},
      "source": [
        "We will apply the patterns to the context of **combinatorial optimization** and show how to solve the **Max-Cut** problem using the **Quantum Approximate Optimization Algorithm (QAOA)**, a hybrid (quantum-classical) iterative method.\n",
        "\n",
        "Note that this QAOA part is based on \"Part 1: Small-scale QAOA\" of the [Quantum approximate optimization algorithm](/docs/tutorials/quantum-approximate-optimization-algorithm) tutorial. See the tutorial to learn how to scale it up.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "1e943b1a-218a-468c-bb63-4269896ebebe",
      "metadata": {},
      "source": [
        "### 3.1 (Small-scale) Qiskit pattern for optimization\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "68fd0b4f-baa4-45dc-9f4c-d9cdff01a651",
      "metadata": {},
      "source": [
        "This part will use a small-scale Max-Cut problem to illustrate the steps required to solve an optimization problem using a quantum computer.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "74b92ba5-c48a-405c-9c4b-04e985a7afbc",
      "metadata": {},
      "source": [
        "The Max-Cut problem is an optimization problem that is hard to solve (more specifically, it is an NP-hard problem) with a number of different applications in clustering, network science, and statistical physics. This tutorial considers a graph of nodes connected by edges and aims to partition the nodes into two sets by \"cutting\" edges, such that the number of edges cut is maximized.\n",
        "\n"
      ]
    },
    {
      "attachments": {},
      "cell_type": "markdown",
      "id": "887eb6ea-f58e-482f-8965-953a08fceecf",
      "metadata": {},
      "source": [
        "![\"Maxcut\"](https://eu-de.quantum.cloud.ibm.com/learning/images/courses/utility-scale-quantum-computing/variational-quantum-algorithms/maxcut.avif)\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "893a25f2",
      "metadata": {},
      "source": [
        "To give some context before mapping this problem to a quantum algorithm, you can better understand how the Max-Cut problem becomes a classical combinatorial optimization problem by first considering the minimization of a function $f(x)$\n",
        "\n",
        "$$\n",
        "\\min_{x\\in \\{0, 1\\}^n}f(x),\n",
        "$$\n",
        "\n",
        "where the input $x$ is a vector whose components correspond to each node of a graph.  Then, constrain each of these components to be either $0$ or $1$ (which represent being included or not included in the cut). This small-scale example case uses a graph with $n=5$ nodes.\n",
        "\n",
        "You could write a function of a pair of nodes $i,j$ which indicates whether the corresponding edge $(i,j)$ is in the cut. For example, the function $x_i + x_j - 2 x_i x_j$ is 1 only if one of either $x_i$ or $x_j$ are 1 (which means that the edge is in the cut) and zero otherwise. The problem of maximizing the edges in the cut can be formulated as\n",
        "\n",
        "$$\n",
        "\\max_{x\\in \\{0, 1\\}^n} \\sum_{(i,j)} x_i + x_j - 2 x_i x_j,\n",
        "$$\n",
        "\n",
        "which can be rewritten as a minimization of the form\n",
        "\n",
        "$$\n",
        "\\min_{x\\in \\{0, 1\\}^n} \\sum_{(i,j)}  2 x_i x_j - x_i - x_j.\n",
        "$$\n",
        "\n",
        "The minimum of $f(x)$ in this case is when the number of edges traversed by the cut is maximal. As you can see, there is nothing relating to quantum computing yet. You need to reformulate this problem into something that a quantum computer can understand.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "e2105a90-027f-44d7-97d1-c2c99373d488",
      "metadata": {},
      "source": [
        "Initialize your problem by creating a graph with $n=5$ nodes.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 22,
      "id": "d3c0dfa7",
      "metadata": {},
      "outputs": [],
      "source": [
        "import matplotlib\n",
        "import matplotlib.pyplot as plt\n",
        "import numpy as np\n",
        "import rustworkx as rx\n",
        "from rustworkx.visualization import mpl_draw"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 23,
      "id": "99e763fe",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Image src=\"/learning/images/courses/utility-scale-quantum-computing/variational-quantum-algorithms/extracted-outputs/99e763fe-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "n = 5\n",
        "\n",
        "graph = rx.PyGraph()\n",
        "graph.add_nodes_from(range(1, n + 1))\n",
        "edge_list = [\n",
        "    (0, 1, 1.0),\n",
        "    (0, 2, 1.0),\n",
        "    (1, 2, 1.0),\n",
        "    (1, 3, 1.0),\n",
        "    (2, 4, 1.0),\n",
        "    (3, 4, 1.0),\n",
        "]\n",
        "graph.add_edges_from(edge_list)\n",
        "pos = rx.spring_layout(graph, seed=2)\n",
        "mpl_draw(graph, node_size=600, pos=pos, with_labels=True, labels=str)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "a06e4386-d7bd-4914-9baa-36a5cc60e3ab",
      "metadata": {},
      "source": [
        "### 3.2 Step 1. Map classical inputs to a quantum problem\n",
        "\n",
        "The first step of the pattern is to map the classical problem (graph) into quantum **circuits** and **operators**. To do this, there are three main steps to take:\n",
        "\n",
        "1. Utilize a series of mathematical reformulations, to represent this problem using the Quadratic Unconstrained Binary Optimization (QUBO) problems notation.\n",
        "2. Rewrite the optimization problem as a Hamiltonian for which the ground state corresponds to the solution which minimizes the cost function.\n",
        "3. Create a quantum circuit which will prepare the ground state of this Hamiltonian via a process similar to quantum annealing.\n",
        "\n",
        "**Note:** In the QAOA methodology, you ultimately want to have an operator (**Hamiltonian**) that represents the **cost function** of our hybrid algorithm, as well as a parametrized circuit (**Ansatz**) that represents quantum states with candidate solutions to the problem. You can sample from these candidate states and then evaluate them using the cost function.\n",
        "\n",
        "#### Graph → optimization problem\n",
        "\n",
        "The first step of the mapping is a notation change, The following expresses the problem in QUBO notation:\n",
        "\n",
        "$$\n",
        "\\min_{x\\in \\{0, 1\\}^n}x^T Q x,\n",
        "$$\n",
        "\n",
        "where $Q$ is a $n\\times n$ matrix of real numbers, $n$ corresponds to the number of nodes in your graph, $x$ is the vector of binary variables introduced above, and $x^T$ indicates the transpose of the vector $x$.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "c00a3493-abab-46d8-86a9-5b853979a575",
      "metadata": {},
      "source": [
        "```\n",
        "Problem name: maxcut\n",
        "\n",
        "Minimize\n",
        "  2*x_1*x_2 + 2*x_1*x_3 + 2*x_2*x_3 + 2*x_2*x_4 + 2*x_3*x_5 + 2*x_4*x_5 - 2*x_1\n",
        "  - 3*x_2 - 3*x_3 - 2*x_4 - 2*x_5\n",
        "\n",
        "Subject to\n",
        "  No constraints\n",
        "\n",
        "  Binary variables (5)\n",
        "    x_1 x_2 x_3 x_4 x_5\n",
        "```\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "a5b9e551-38a1-4543-b9f1-caaefb0ef3a9",
      "metadata": {},
      "source": [
        "#### Optimization problem → Hamiltonian\n",
        "\n",
        "You can then reformulate the QUBO problem as a **Hamiltonian** (here, a matrix that represents the energy of a system):\n",
        "\n",
        "$$\n",
        "H_C=\\sum_{ij}Q_{ij}Z_iZ_j + \\sum_i b_iZ_i.\n",
        "$$\n",
        "\n",
        "**Reformulation steps from the QAOA problem to the Hamiltonian**\n",
        "\n",
        "To demonstrate how the QAOA problem can be rewritten in this way, first replace the binary variables $x_i$ to a new set of variables $z_i\\in\\{-1, 1\\}$ via\n",
        "\n",
        "$$\n",
        "x_i = \\frac{1-z_i}{2}.\n",
        "$$\n",
        "\n",
        "Here you can see that if $x_i$ is $0$, then $z_i$ must be $1$. When the $x_i$'s are substituted for the $z_i$'s in the optimization problem ($x^TQx$), an equivalent formulation can be obtained.\n",
        "\n",
        "$$\n",
        "x^TQx=\\sum_{ij}Q_{ij}x_ix_j \\\\ =\\frac{1}{4}\\sum_{ij}Q_{ij}(1-z_i)(1-z_j) \\\\=\\frac{1}{4}\\sum_{ij}Q_{ij}z_iz_j-\\frac{1}{4}\\sum_{ij}(Q_{ij}+Q_{ji})z_i + \\frac{n^2}{4}.\n",
        "$$\n",
        "\n",
        "Now if we define $b_i=-\\sum_{j}(Q_{ij}+Q_{ji})$, remove the prefactor, and the constant $n^2$ term, we arrive at the two equivalent formulations of the same optimization problem.\n",
        "\n",
        "$$\n",
        "min_{x\\in\\{0,1\\}^n} x^TQx\\Longleftrightarrow \\min_{z\\in\\{-1,1\\}^n}z^TQz + b^Tz\n",
        "$$\n",
        "\n",
        "Here, $b$ depends on $Q$. Note that to obtain $z^TQz + b^Tz$ we dropped the factor of 1/4 and a constant offset of $n^2$ which do not play a role in the optimization.\n",
        "\n",
        "Now, to obtain a quantum formulation of the problem, promote the $z_i$ variables to a Pauli $Z$ matrix, such as a $2\\times 2$ matrix of the form\n",
        "\n",
        "$$\n",
        "Z_i = \\begin{pmatrix}1 & 0 \\\\ 0 & -1\\end{pmatrix}.\n",
        "$$\n",
        "\n",
        "When you substitute these matrices in the optimization problem above, you obtain the following Hamiltonian\n",
        "\n",
        "$$\n",
        "H_C=\\sum_{ij}Q_{ij}Z_iZ_j + \\sum_i b_iZ_i.\n",
        "$$\n",
        "\n",
        "*Also recall that the $Z$ matrices are embedded in the quantum computer's computational space, that is, a Hilbert space of size $2^n\\times 2^n$. Therefore, you should understand terms such as $Z_iZ_j$ as the tensor product $Z_i\\otimes Z_j$ embedded in the $2^n\\times 2^n$ Hilbert space. For example, in a problem with five decision variables the term $Z_1Z_3$ is understood to mean $I\\otimes Z_3\\otimes I\\otimes Z_1\\otimes I$ where $I$ is the $2\\times 2$ identity matrix.*\n",
        "\n",
        "This Hamiltonian is called the <b>cost function Hamiltonian</b>. It has the property that its ground state corresponds to the solution that <b>minimizes the cost function $f(x)$</b>.\n",
        "Therefore, to solve your optimization problem you now need to prepare the ground state of $H_C$ (or a state with a high overlap with it) on the quantum computer. Then, sampling from this state will, with a high probability, yield the solution to $\\min~f(x)$.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 24,
      "id": "47256f53-8e02-494e-864a-a6f45b10442a",
      "metadata": {},
      "outputs": [],
      "source": [
        "def build_max_cut_operator(graph: rx.PyGraph) -> tuple[SparsePauliOp, float]:\n",
        "    sp_list = []\n",
        "    constant = 0\n",
        "    for s, t in graph.edge_list():\n",
        "        w = graph.get_edge_data(s, t)\n",
        "        sp_list.append((\"ZZ\", [s, t], w / 2))\n",
        "        constant -= 1 / 2\n",
        "    return SparsePauliOp.from_sparse_list(\n",
        "        sp_list, num_qubits=graph.num_nodes()\n",
        "    ), constant"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 25,
      "id": "01a2d8eb-b63b-40bc-93c5-0b547c06b194",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Cost Function Hamiltonian: SparsePauliOp(['IIIZZ', 'IIZIZ', 'IIZZI', 'IZIZI', 'ZIZII', 'ZZIII'],\n",
            "              coeffs=[0.5+0.j, 0.5+0.j, 0.5+0.j, 0.5+0.j, 0.5+0.j, 0.5+0.j])\n",
            "Constant: -3.0\n"
          ]
        }
      ],
      "source": [
        "cost_hamiltonian, constant = build_max_cut_operator(graph)\n",
        "print(\"Cost Function Hamiltonian:\", cost_hamiltonian)\n",
        "print(\"Constant:\", constant)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "33f71b0d-4a2a-4082-8c1a-ce9d2b769048",
      "metadata": {},
      "source": [
        "#### Hamiltonian → quantum circuit\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "00431c46-30c2-40f9-99df-40baf8da98f6",
      "metadata": {},
      "source": [
        "The Hamiltonian $H_C$ contains the quantum definition of your problem. Now you can create a quantum circuit that will help *sample* good solutions from the quantum computer. The QAOA is inspired by quantum annealing and applies alternating layers of operators in the quantum circuit.\n",
        "\n",
        "The general idea is to start in the ground state of a known system, $H^{\\otimes n}|0\\rangle$ above, and then steer the system into the ground state of the cost operator that you are interested in. This is done by applying the operators $\\exp\\{-i\\gamma_k H_C\\}$ and $\\exp\\{-i\\beta_k H_m\\}$ with angles $\\gamma_1,...,\\gamma_p$ and $\\beta_1,...,\\beta_p~$.\n",
        "\n",
        "The quantum circuit that you generate is **parametrized** by $\\gamma_i$ and $\\beta_i$, so you can try out different values of $\\gamma_i$ and $\\beta_i$ and sample from the resulting state.\n",
        "\n"
      ]
    },
    {
      "attachments": {},
      "cell_type": "markdown",
      "id": "ca09d6bf-e421-4ada-9515-1f69687bb511",
      "metadata": {},
      "source": [
        "![\"QAOA circuit diagram\"](https://eu-de.quantum.cloud.ibm.com/learning/images/courses/utility-scale-quantum-computing/variational-quantum-algorithms/circuit-diagram.svg)\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "0d12bc53-7805-43e4-b4cb-754fd234b519",
      "metadata": {},
      "source": [
        "In this case we will try an example with 1 QAOA layer that contains two parameters: $\\gamma_1$ and $\\beta_1$.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 26,
      "id": "1f6215c0",
      "metadata": {},
      "outputs": [],
      "source": [
        "from qiskit.circuit.library import QAOAAnsatz"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 27,
      "id": "7bd8c6d4-f40f-4a11-a440-0b26d9021b53",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Image src=\"/learning/images/courses/utility-scale-quantum-computing/variational-quantum-algorithms/extracted-outputs/7bd8c6d4-f40f-4a11-a440-0b26d9021b53-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "execution_count": 27,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "circuit = QAOAAnsatz(cost_operator=cost_hamiltonian, reps=1)\n",
        "circuit.measure_all()\n",
        "circuit.draw(\"mpl\")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 28,
      "id": "148d2d62",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Image src=\"/learning/images/courses/utility-scale-quantum-computing/variational-quantum-algorithms/extracted-outputs/148d2d62-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "execution_count": 28,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "circuit.decompose(reps=3).draw(\"mpl\", fold=-1)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 29,
      "id": "315c495a",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "ParameterView([ParameterVectorElement(β[0]), ParameterVectorElement(γ[0])])"
            ]
          },
          "execution_count": 29,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "circuit.parameters"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "82f70daa-ff68-447a-8064-8b7df7a646cf",
      "metadata": {},
      "source": [
        "### 3.3 Step 2. Optimize circuits for quantum hardware execution\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "c08be444-e3ed-4178-a10b-414069b1b411",
      "metadata": {},
      "source": [
        "The circuit above contains a series of abstractions useful to think about quantum algorithms, but not possible to run on the hardware. To be able to run on a QPU, the circuit needs to undergo a series of operations that make up the **transpilation** or **circuit optimization** step of the pattern.\n",
        "\n",
        "The Qiskit library offers a series of **transpilation passes** that cater to a wide range of circuit transformations. You need to make sure that your circuit is **optimized** for your purpose.\n",
        "\n",
        "Transpilation might involves several steps, such as:\n",
        "\n",
        "* **Initial mapping** of the qubits in the circuit (such as decision variables) to physical qubits on the device.\n",
        "* **Unrolling** of the instructions in the quantum circuit to the hardware-native instructions that the backend understands.\n",
        "* **Routing** of any qubits in the circuit that interact to physical qubits that are adjacent with one another.\n",
        "* **Error suppression** by adding single-qubit gates to suppress noise with dynamical decoupling.\n",
        "\n",
        "More information about transpilation is available in our [documentation](/docs/guides/transpile).\n",
        "\n",
        "The following code transforms and optimizes the abstract circuit into a format that is ready for execution on one of devices accessible through the cloud using the **Qiskit IBM® Runtime service**.\n",
        "\n",
        "Note that you can test your programs locally by \"local testing mode\" before sending them to real quantum computers.\n",
        "More information about the local testing mode is available in the [documentation.](/docs/guides/local-testing-mode)\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "95cd3eed-0348-4373-b664-16a65d42f1e7",
      "metadata": {},
      "outputs": [
        {
          "name": "stderr",
          "output_type": "stream",
          "text": [
            "  service = QiskitRuntimeService(channel=\"ibm_quantum_platform\")\n"
          ]
        },
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "<IBMBackend('ibm_strasbourg')>\n"
          ]
        },
        {
          "data": {
            "text/plain": [
              "<Image src=\"/learning/images/courses/utility-scale-quantum-computing/variational-quantum-algorithms/extracted-outputs/95cd3eed-0348-4373-b664-16a65d42f1e7-2.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "execution_count": 31,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "from qiskit_ibm_runtime import QiskitRuntimeService\n",
        "from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager\n",
        "\n",
        "# Use a quantum device\n",
        "service = QiskitRuntimeService()\n",
        "backend = service.least_busy(min_num_qubits=127)\n",
        "# backend = service.backend(\"ibm_kingston\")\n",
        "\n",
        "# You can test your programs locally with a fake backend (local testing mode)\n",
        "# backend = FakeBrisbane()\n",
        "\n",
        "print(backend)\n",
        "\n",
        "# Create pass manager for transpilation\n",
        "pm = generate_preset_pass_manager(optimization_level=3, backend=backend)\n",
        "\n",
        "candidate_circuit = pm.run(circuit)\n",
        "candidate_circuit.draw(\"mpl\", fold=False, idle_wires=False)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "4e75cad7-f599-4937-b5fe-f4d01f53423c",
      "metadata": {},
      "source": [
        "### 3.4 Step 3. Execute using Qiskit primitives\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "9b99ce67-f121-4244-b62a-536be38fea86",
      "metadata": {},
      "source": [
        "In the QAOA workflow, the optimal QAOA parameters are found in an iterative optimization loop, which runs a series of circuit evaluations and uses a classical optimizer to find the optimal $\\beta_k$ and $\\gamma_k$ parameters. This execution loop is executed via the following steps:\n",
        "\n",
        "1. Define the initial parameters\n",
        "2. Instantiate a new `Session` containing the optimization loop and the primitive used to sample the circuit\n",
        "3. Once an optimal set of parameters is found, execute the circuit a final time to obtain a final distribution which will be used in the post-process step.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "00b2b0f1-9bad-4ad3-b93e-5cbf40395dbf",
      "metadata": {},
      "source": [
        "#### Define circuit with initial parameters\n",
        "\n",
        "We start with arbitrary chosen parameters.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 32,
      "id": "afa5747f-44dc-4e41-a875-7b6f896f13e2",
      "metadata": {},
      "outputs": [],
      "source": [
        "initial_gamma = np.pi\n",
        "initial_beta = np.pi / 2\n",
        "init_params = [initial_gamma, initial_beta]"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "b867f1b0-7196-4d34-9b28-e3fb1de8221c",
      "metadata": {},
      "source": [
        "#### Define backend and execution primitive\n",
        "\n",
        "Use the **Qiskit Runtime primitives** to interact with IBM® backends. The two primitives are Sampler and Estimator, and the choice of primitive depends on what type of measurement you want to run on the quantum computer. For the minimization of $H_C$, use the Estimator since the measurement of the cost function is simply the expectation value of $\\langle H_C \\rangle$.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "36c789f2",
      "metadata": {},
      "source": [
        "#### Run\n",
        "\n",
        "The primitives offer a variety of [execution modes](/docs/guides/execution-modes) to schedule workloads on quantum devices, and a QAOA workflow runs iteratively in a session.\n",
        "\n"
      ]
    },
    {
      "attachments": {},
      "cell_type": "markdown",
      "id": "122b7dc8-8d8a-45f2-813e-6199905d765b",
      "metadata": {},
      "source": [
        "![\"execution mode\"](https://eu-de.quantum.cloud.ibm.com/learning/images/courses/utility-scale-quantum-computing/variational-quantum-algorithms/execution-mode.avif)\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "fff58deb",
      "metadata": {},
      "source": [
        "You can plug the sampler-based cost function into the SciPy minimizing routine to find the optimal parameters.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 33,
      "id": "ff947109-cddc-4d3c-9119-2c729df73115",
      "metadata": {},
      "outputs": [],
      "source": [
        "def cost_func_estimator(params, ansatz, hamiltonian, estimator):\n",
        "    # transform the observable defined on virtual qubits to\n",
        "    # an observable defined on all physical qubits\n",
        "    isa_hamiltonian = hamiltonian.apply_layout(ansatz.layout)\n",
        "\n",
        "    pub = (ansatz, isa_hamiltonian, params)\n",
        "    job = estimator.run([pub])\n",
        "\n",
        "    results = job.result()[0]\n",
        "    cost = results.data.evs\n",
        "\n",
        "    objective_func_vals.append(cost)\n",
        "\n",
        "    return cost"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 34,
      "id": "f5f46775",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            " message: Optimization terminated successfully.\n",
            " success: True\n",
            "  status: 1\n",
            "     fun: -0.6557925874481715\n",
            "       x: [ 2.873e+00  9.414e-01]\n",
            "    nfev: 21\n",
            "   maxcv: 0.0\n"
          ]
        }
      ],
      "source": [
        "from qiskit_ibm_runtime import Session, EstimatorV2\n",
        "from scipy.optimize import minimize\n",
        "\n",
        "objective_func_vals = []  # Global variable\n",
        "with Session(backend=backend) as session:\n",
        "    # If using qiskit-ibm-runtime<0.24.0, change `mode=` to `session=`\n",
        "    estimator = EstimatorV2(mode=session)\n",
        "    estimator.options.default_shots = 1000\n",
        "\n",
        "    # Set simple error suppression/mitigation options\n",
        "    estimator.options.dynamical_decoupling.enable = True\n",
        "    estimator.options.dynamical_decoupling.sequence_type = \"XY4\"\n",
        "    estimator.options.twirling.enable_gates = True\n",
        "    estimator.options.twirling.num_randomizations = \"auto\"\n",
        "\n",
        "    result = minimize(\n",
        "        cost_func_estimator,\n",
        "        init_params,\n",
        "        args=(candidate_circuit, cost_hamiltonian, estimator),\n",
        "        method=\"COBYLA\",\n",
        "        tol=1e-2,\n",
        "    )\n",
        "    print(result)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "ad878b62",
      "metadata": {},
      "source": [
        "The optimizer was able to reduce the cost and find better parameters for the circuit.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 35,
      "id": "f923dd5d",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Image src=\"/learning/images/courses/utility-scale-quantum-computing/variational-quantum-algorithms/extracted-outputs/f923dd5d-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "plt.figure(figsize=(12, 6))\n",
        "plt.plot(objective_func_vals)\n",
        "plt.xlabel(\"Iteration\")\n",
        "plt.ylabel(\"Cost\")\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "e2ea6359",
      "metadata": {},
      "source": [
        "Once you have found the optimal parameters for the circuit, you can assign these parameters and sample the final distribution obtained with the optimized parameters. Here is where the *Sampler* primitive should be used since it is the probability distribution of bitstring measurements which correspond to the optimal cut of the graph.\n",
        "\n",
        "**Note:** This means preparing a quantum state $\\psi$ in the computer and then measuring it. A measurement will collapse the state into a single computational basis state - for example, `010101110000...` - which corresponds to a candidate solution $x$ to our initial optimization problem ($\\max f(x)$ or $\\min f(x)$ depending on the task).\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 36,
      "id": "f8dddf5a",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Image src=\"/learning/images/courses/utility-scale-quantum-computing/variational-quantum-algorithms/extracted-outputs/f8dddf5a-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "execution_count": 36,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "optimized_circuit = candidate_circuit.assign_parameters(result.x)\n",
        "optimized_circuit.draw(\"mpl\", fold=False, idle_wires=False)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 37,
      "id": "fd9669cf",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "{12: 0.0652, 31: 0.0089, 4: 0.0085, 13: 0.0731, 26: 0.0256, 28: 0.0246, 17: 0.0405, 25: 0.0591, 20: 0.031, 15: 0.0221, 8: 0.017, 21: 0.0371, 14: 0.0461, 16: 0.0229, 19: 0.0723, 23: 0.0199, 22: 0.0478, 18: 0.0708, 24: 0.0165, 6: 0.0525, 7: 0.0155, 5: 0.0245, 3: 0.0231, 29: 0.0121, 30: 0.0062, 10: 0.0363, 1: 0.0097, 9: 0.042, 27: 0.0094, 11: 0.0349, 0: 0.0129, 2: 0.0119}\n"
          ]
        }
      ],
      "source": [
        "from qiskit_ibm_runtime import SamplerV2\n",
        "\n",
        "# If using qiskit-ibm-runtime<0.24.0, change `mode=` to `backend=`\n",
        "sampler = SamplerV2(mode=backend)\n",
        "\n",
        "# Set simple error suppression/mitigation options\n",
        "sampler.options.dynamical_decoupling.enable = True\n",
        "sampler.options.dynamical_decoupling.sequence_type = \"XY4\"\n",
        "sampler.options.twirling.enable_gates = True\n",
        "sampler.options.twirling.num_randomizations = \"auto\"\n",
        "\n",
        "pub = (optimized_circuit,)\n",
        "job = sampler.run([pub], shots=int(1e4))\n",
        "counts_int = job.result()[0].data.meas.get_int_counts()\n",
        "counts_bin = job.result()[0].data.meas.get_counts()\n",
        "shots = sum(counts_int.values())\n",
        "final_distribution_int = {key: val / shots for key, val in counts_int.items()}\n",
        "final_distribution_bin = {key: val / shots for key, val in counts_bin.items()}\n",
        "print(final_distribution_int)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "2c89613f",
      "metadata": {},
      "source": [
        "### 3.5 Step 4. Post-process, return result in classical format\n",
        "\n",
        "The post-processing step interprets the sampling output to return a solution for your original problem. In this case, you are interested in the bitstring with the highest probability as this determines the optimal cut. The symmetries in the problem allow for four possible solutions, and the sampling process will return one of them with a slightly higher probability, but you can see in the plotted distribution below that four of the bitstrings are distinctively more likely than the rest.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 38,
      "id": "d4f7fc70-883f-4b6b-8e92-2fc4afbbea46",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Result bitstring: [1, 0, 1, 1, 0]\n"
          ]
        }
      ],
      "source": [
        "# auxiliary functions to sample most likely bitstring\n",
        "def to_bitstring(integer, num_bits):\n",
        "    result = np.binary_repr(integer, width=num_bits)\n",
        "    return [int(digit) for digit in result]\n",
        "\n",
        "\n",
        "keys = list(final_distribution_int.keys())\n",
        "values = list(final_distribution_int.values())\n",
        "most_likely = keys[np.argmax(np.abs(values))]\n",
        "most_likely_bitstring = to_bitstring(most_likely, len(graph))\n",
        "most_likely_bitstring.reverse()\n",
        "\n",
        "print(\"Result bitstring:\", most_likely_bitstring)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 39,
      "id": "32a3020e-c1ea-4aff-988c-d7910a690fa8",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Image src=\"/learning/images/courses/utility-scale-quantum-computing/variational-quantum-algorithms/extracted-outputs/32a3020e-c1ea-4aff-988c-d7910a690fa8-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "import matplotlib.pyplot as plt\n",
        "\n",
        "matplotlib.rcParams.update({\"font.size\": 10})\n",
        "final_bits = final_distribution_bin\n",
        "values = np.abs(list(final_bits.values()))\n",
        "top_4_values = sorted(values, reverse=True)[:4]\n",
        "positions = []\n",
        "for value in top_4_values:\n",
        "    positions.append(np.where(values == value)[0])\n",
        "fig = plt.figure(figsize=(11, 6))\n",
        "ax = fig.add_subplot(1, 1, 1)\n",
        "plt.xticks(rotation=45)\n",
        "plt.title(\"Result Distribution\")\n",
        "plt.xlabel(\"Bitstrings (reversed)\")\n",
        "plt.ylabel(\"Probability\")\n",
        "ax.bar(list(final_bits.keys()), list(final_bits.values()), color=\"tab:grey\")\n",
        "for p in positions:\n",
        "    ax.get_children()[p[0].item()].set_color(\"tab:purple\")\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "6cfcb278",
      "metadata": {},
      "source": [
        "#### Visualize best cut\n",
        "\n",
        "From the optimal bit string, you can then visualize this cut on the original graph.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 40,
      "id": "22a48124-e6b4-4144-bee1-f01fa4c7ccbb",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Image src=\"/learning/images/courses/utility-scale-quantum-computing/variational-quantum-algorithms/extracted-outputs/22a48124-e6b4-4144-bee1-f01fa4c7ccbb-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "colors = [\"tab:grey\" if i == 0 else \"tab:purple\" for i in most_likely_bitstring]\n",
        "mpl_draw(graph, node_size=600, pos=pos, with_labels=True, labels=str, node_color=colors)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "4c4803ab",
      "metadata": {},
      "source": [
        "And calculate the value of the cut. The solution is not optimal due to noise (the cut value of the optimal solution is 5).\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 41,
      "id": "7208ee7d",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "The value of the cut is: 5\n"
          ]
        }
      ],
      "source": [
        "from typing import Sequence\n",
        "\n",
        "\n",
        "def evaluate_sample(x: Sequence[int], graph: rx.PyGraph) -> float:\n",
        "    assert len(x) == len(\n",
        "        list(graph.nodes())\n",
        "    ), \"The length of x must coincide with the number of nodes in the graph.\"\n",
        "    return sum(\n",
        "        x[u] * (1 - x[v]) + x[v] * (1 - x[u]) for u, v in list(graph.edge_list())\n",
        "    )\n",
        "\n",
        "\n",
        "cut_value = evaluate_sample(most_likely_bitstring, graph)\n",
        "print(\"The value of the cut is:\", cut_value)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "1450be7e",
      "metadata": {},
      "source": [
        "This concludes the small-scale QAOA tutorial.\n",
        "You will learn how to adapt QAOA at a utility-scale in \"Part 2: scale it up!\" of the [Quantum approximate optimization algorithm](/docs/tutorials/quantum-approximate-optimization-algorithm) tutorial.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 42,
      "id": "2a4f85ab",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "'2.0.2'"
            ]
          },
          "execution_count": 42,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "# Check Qiskit version\n",
        "import qiskit\n",
        "\n",
        "qiskit.__version__"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "id": "a1b8767d",
      "source": "© IBM Corp., 2017-2026"
    }
  ],
  "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"
    },
    "widgets": {
      "application/vnd.jupyter.widget-state+json": {
        "state": {},
        "version_major": 2,
        "version_minor": 0
      }
    }
  },
  "nbformat": 4,
  "nbformat_minor": 5
}