{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Quantum Approximate Optimization Algorithm\n",
        "\n",
        "This notebook introduces the API of the QAOA implementation available in the `qat.opt` namespace.\n",
        "\n",
        "We assume that the user is familiar with the algorithm itself and only detail the API of the library and the few algorithmical tools used for circuit synthesis.\n",
        "\n",
        "\n",
        "## The CombinatorialProblem class\n",
        "\n",
        "The interface of the library is concentrated into a single class `CombinatorialProblem`.\n",
        "\n",
        "This class allows to:\n",
        "* declare boolean variables\n",
        "* add new clauses (i.e boolean formulae) to the final cost function"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "from qat.opt import CombinatorialProblem\n",
        "\n",
        "# Declaring a fresh problem\n",
        "my_problem = CombinatorialProblem()\n",
        "\n",
        "# Declaring a new variable\n",
        "v0 = my_problem.new_var()\n",
        "v1 = my_problem.new_var()\n",
        "\n",
        "# Or several variables\n",
        "v_array = my_problem.new_vars(4)\n",
        "\n",
        "# Variable are indexed starting from 0\n",
        "print(v0, v1)\n",
        "print(\", \".join(str(v) for v in v_array))\n",
        "# Clauses are built using boolean operators (|, &, ^, ~) and variables\n",
        "print(v0 | v1)\n",
        "print(v_array[0] & v_array[2])\n",
        "print(v0 ^ v_array[0])\n",
        "print(~v0)\n",
        "print(~(v0 ^v_array[3] | v1))\n",
        "\n",
        "# Clauses are added to a problem using the `add_clause` method\n",
        "my_problem.add_clause(v0 ^ v1)\n",
        "# Clauses can be weighted\n",
        "my_problem.add_clause(v0 | v1, weight=2.)\n",
        "for clause, weight in my_problem.clauses:\n",
        "    print(clause, weight)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "By default, the class assumes that the described problem is a minimization problem.\n",
        "It is possible to specify maximization problems by adding an argument in the constructor.\n",
        "\n",
        "In practice, this will simply flip the sign of the cost function (or more precisely, its Hamiltonian encoding)."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "my_maximization_problem = CombinatorialProblem(maximization=True)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## From problem to variational Ansätze\n",
        "\n",
        "Once a problem is declared, it is straightforward to construct a QAOA variational Ansatz from it: "
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "# We just need to specify a number of layers\n",
        "my_problem = CombinatorialProblem()\n",
        "variables = my_problem.new_vars(5)\n",
        "for i in range(4):\n",
        "    my_problem.add_clause(variables[i]^variables[i+1])\n",
        "\n",
        "depth = 3\n",
        "ansatz = my_problem.to_job(\"qaoa\", depth).circuit\n",
        "ansatz.display()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The variational ansatz is parametrized by abstract variables $\\gamma_0,...,\\gamma_{l-1}$ and $\\beta_0,...,\\beta_{l-1}$.\n",
        "\n",
        "Variables can be listed as follows:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "print(\"Variables:\", ansatz.get_variables())"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "You can see that their name is latex compliant, just for a nice display.\n",
        "\n",
        "It is possible to bind these variables using their names:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "import numpy as np\n",
        "ansatz_gamma_0_pi = ansatz.bind_variables({\"\\\\gamma_{0}\": np.pi})\n",
        "# or equivalently\n",
        "ansatz_gamma_0_pi = ansatz(**{\"\\\\gamma_{0}\": np.pi})\n",
        "ansatz_gamma_0_pi.display()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Observable synthesis\n",
        "\n",
        "In order to be able to generate a QAOA Anstaz the Problem class first encodes each clause into a small Hamiltonian using the following inductive definition:\n",
        "If boolean clauses are represented using the following grammar:\n",
        "\n",
        "$exp := exp \\lor exp | exp \\land exp | exp \\oplus exp | \\neg exp | V$\n",
        "\n",
        "Then the Hamiltonian encoding proceeds as follow:\n",
        "\n",
        "$H(e_1\\lor e_2) = H(e_1) + H(e_2) - H(e_1)H(e_2)$\n",
        "\n",
        "$H(e_1 \\land e_2) = H(e_1) * H(e_2)$\n",
        "\n",
        "$H(e1 \\oplus e2) = H(e1) + H(e2) - 2H(e1)H(e2)$\n",
        "\n",
        "$H(\\neg e) = 1 - H(e)$\n",
        "\n",
        "$H(V(i)) = \\frac{1 - \\sigma_i^z}{2}$\n",
        "\n",
        "The complete encoding is then obtained by summing these smaller Hamiltonian (with some eventual coefficients to account for the weights). \n",
        "\n",
        "Finally, if the problem is a maximization problem, the sign of the Hamiltonian is flipped, so that the problem becomes a minimization problem.\n",
        "\n",
        "The Hamiltonian, i.e. the *terms*-type of Observable for the problem, can be obtained using the `.get_observable()` method with the `\"terms\"` argument:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "my_problem = CombinatorialProblem()\n",
        "variables = my_problem.new_vars(5)\n",
        "for i in range(4):\n",
        "    my_problem.add_clause(variables[i]^variables[i+1])\n",
        "print(\"Minimization:\\n\", my_problem.get_observable(\"terms\"))\n",
        "\n",
        "my_problem = CombinatorialProblem(maximization=True)\n",
        "variables = my_problem.new_vars(5)\n",
        "for i in range(4):\n",
        "    my_problem.add_clause(variables[i]^variables[i+1])\n",
        "print(\"Maximization:\\n\", my_problem.get_observable(\"terms\"))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Circuit synthesis\n",
        "\n",
        "Once the observable is generated, there are two distinct circuit synthesis algorithm that can be used to extract an Ansatz from the cost Hamiltonian:\n",
        "\n",
        "* The \"default\" algorithm naively produces a subcircuit per term in the Hamiltonian for each layer of the Ansatz. For most applications, this algorithm is enough and will provide a relatively efficient Ansatz.\n",
        "\n",
        "* The \"coloring\" heuristics does pretty much the same but optimizes the ordering of the terms in order to minimize circuit depth.\n",
        "\n",
        "* The \"gray_synth\" heuristics uses Amy et al phase polynomial synthesis algorithm to implement the entangling portion of the Ansatz. This can help reduce the CNOT count of the resulting circuit."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "my_problem = CombinatorialProblem()\n",
        "n = 7\n",
        "variables = my_problem.new_vars(n)\n",
        "for i in range(n - 2):\n",
        "    my_problem.add_clause(variables[i] ^ variables[i+1] ^ variables[i+2])\n",
        "print(\"Cost Hamiltonian:\\n\", my_problem.get_observable(\"terms\"))\n",
        "\n",
        "circuit1 = my_problem.to_job(\"qaoa\", 1, strategy=\"default\").circuit # '1' is the depth\n",
        "circuit2 = my_problem.to_job(\"qaoa\", 1, strategy=\"coloring\").circuit\n",
        "circuit3 = my_problem.to_job(\"qaoa\", 1, strategy=\"gray_synth\").circuit\n",
        "circuit1.display()\n",
        "circuit2.display()\n",
        "circuit3.display()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Running a full algorithm - MaxCut\n",
        "\n",
        "The `qat.opt` namespace also contains a very simple wrapper to produce problems describing a MAXCUT instance.\n",
        "\n",
        "The class can be instantiated using a networkx graph:\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "import networkx as nx\n",
        "\n",
        "graph = nx.generators.random_graphs.erdos_renyi_graph(10, 0.5)\n",
        "nx.draw(graph)\n",
        "from qat.opt import MaxCut\n",
        "problem = MaxCut(graph)\n",
        "print(problem)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can now make use of a variational plugin to optimize a QAOA Ansatz:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "from qat.qpus import get_default_qpu\n",
        "from qat.plugins import ScipyMinimizePlugin\n",
        "qpu = get_default_qpu()\n",
        "stack = ScipyMinimizePlugin(method=\"COBYLA\",\n",
        "                            tol=1e-5, \n",
        "                            options={\"maxiter\": 200}) | qpu\n",
        "# We can directly call the to_job method of the Problem class to pack an Ansatz and \n",
        "# the cost observable in a single abstract Job\n",
        "job = problem.to_job(\"qaoa\", 3) # Here 3 is the depth of the Ansatz\n",
        "result = stack.submit(job)\n",
        "print(\"Final energy:\", result.value)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt\n",
        "plt.plot(eval(result.meta_data[\"optimization_trace\"]))\n",
        "plt.xlabel(\"steps\")\n",
        "plt.ylabel(\"energy\")\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Printing the most probable states\n",
        "\n",
        "We now print the most probable states of the distribution corresponding to the optimized QAOA parameters:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {},
      "outputs": [],
      "source": [
        "#Binding the variables:\n",
        "sol_job = job(**eval(result.meta_data[\"parameter_map\"]))\n",
        "\n",
        "#Rerunning in 'SAMPLE' mode to get the most probable states:\n",
        "sampling_job = sol_job.circuit.to_job()\n",
        "\n",
        "sol_res = qpu.submit(sampling_job)\n",
        "print(\"Most probable states are:\")\n",
        "for sample in sol_res:\n",
        "    if sample.probability > 0.05:\n",
        "        print(sample.state, \"{:.2f}%\".format(100 * sample.probability))\n",
        "        \n",
        "# We can also directly cast states into bitstrings for practical use:\n",
        "print(\"And as bitstrings:\")\n",
        "for sample in sol_res:\n",
        "    if sample.probability > 0.05:\n",
        "        print(sample.state.value, \"{:.2f}%\".format(100 * sample.probability))\n",
        "        \n",
        "# Emulating a reasonnable setup:\n",
        "# Drawing 1024 cuts\n",
        "sampling_job = sol_job.circuit.to_job(nbshots=1024)\n",
        "sol_res = qpu.submit(sampling_job)\n",
        "# Picking the most probable cut\n",
        "best_cut = max([(s.state.value[0], s.probability) for s in sol_res], key=lambda s: s[1])[0]\n",
        "print(\n",
        "    \"Most probable cut:\", \n",
        "    set(i for i in graph.nodes() if best_cut[i] == '1'), \n",
        "    set(i for i in graph.nodes() if best_cut[i] == '0')\n",
        ")\n",
        "# Plotting the cut nicely\n",
        "nx.draw(\n",
        "    graph, labels={i: str(i) for i in graph.nodes()},\n",
        "    node_color=[\"cyan\" if b == '1' else 'lightgreen' for b in best_cut],\n",
        "    edge_color=['red' if best_cut[a] != best_cut[b] else 'black' for a, b in graph.edges()]\n",
        ")"
      ]
    }
  ],
  "metadata": {
    "authors": [
      "Simon Martiel"
    ],
    "kernelspec": {
      "display_name": "Python 3 (ipykernel)",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.12.11"
    },
    "tags": [
      "variational"
    ]
  },
  "nbformat": 4,
  "nbformat_minor": 2
}