2_customizing_workflows.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "ff380785",
   "metadata": {},
   "source": [
    "Customizing Workflows\n",
    "===================\n",
    "\n",
    "This notebook shows how to write workflows for BQSKit, using the Hubbard Model$^1$ as an exemplar. The Hubbard Model is only briefly introduced; full details on the Hubbard Model, the Variational Quantum Eigensolver, the effects of noise, and more, can be found as part of the scikit-quant tutorials.$^2$\n",
    "\n",
    "Prerequisites: OpenFermion, Qiskit, BQSKit, and SciPy (which will pull in NumPy), as well as optionally scikit-quant, e.g. install with pip:\n",
    "\n",
    "`$ python -m pip install openfermion qiskit bqskit scipy`\n",
    "\n",
    "If doing the noisy simulation part of the tutorial, then qiskit_algorithms, qiskit-aer, and scikit-quant are needed as well, e.g. install with pip:\n",
    "\n",
    "`$ python -m pip install qiskit_algorithms qiskit-aer scikit-quant`"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2c306def",
   "metadata": {},
   "source": [
    "Contents\n",
    "-----------\n",
    "\n",
    "This tutorial is setup to be run serially for the first several sections, as lessons learned (and code) are used in subsequent sections. However, afterwards it is possible to skip and jump ahead if a specific focus is of more interest.\n",
    "\n",
    "* Work through in series:\n",
    "   * [Introduction to the Hubbard Model](#Introduction)\n",
    "   * [Effects of circuit structure on synthesis](#Energy-Objective)\n",
    "   * [Running BQSKit synthesis in VQE](#Variational-Quantum-Eigensolver)\n",
    "* Work through in series, or pick and choose:\n",
    "   * [Customizing BQSKit workflows](#Customized-Compilation)\n",
    "   * [Seeded synthesis](#Seeded-Synthesis)\n",
    "   * [Implications of noise](#Implications-of-Noise)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "635652ff",
   "metadata": {},
   "source": [
    "Introduction\n",
    "--------------\n",
    "\n",
    "We will solve the Hubbard Model using the Variational Quantum Eigensolver (VQE) algorithm. The Hubbard Model (HM) has a range of uses, most notably in solid-state physics to explain phase transitions and other correlated behaviors. There are several variations of the HM, but in essence, the model describes interacting particles on a lattice. Despite its apparent simplicity, it is hard to solve on a classical computer because the known algorithms scale exponentially unless further approximations are made. In this tutorial, we will use the VQE algorithm to estimate the energy ground states of different configurations of an HM."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "276bca9d",
   "metadata": {},
   "source": [
    "Basic Fermi-Hubbard Model\n",
    "---------------------------------\n",
    "\n",
    "We will restrict ourselves to a Hubbard Model with fermions (such as electrons), i.e. the Fermi-Hubbard Model. A visualization of an example with a 2x2 grid with spinful electrons is shown below. Electrons being fermions means that two electrons can not occupy the same quantum state. Because the electrons have spin, there can thus be at most two electrons per site (one spin-up, one spin-down). The model is characterized by an interaction term 'U', which only contributes if a spin-up and a spin-down electron occupy the same site; and by a \"hopping\" term 't', representing the kinetic energy in the system. In this basic model, only hops to neighboring sites are considered.\n",
    "\n",
    "<table style=\"padding: 0px\"><tr></tr><tr>\n",
    "<td style=\"width: 65%; padding: 0px\">\n",
    "<div align=\"left\" style=\"text-align: left; font-size: 120%\">\n",
    "The energy in the model is described by the following Hamiltonian using fermion algebra:\n",
    "<br>\n",
    "<p>\\begin{equation*}\n",
    "H = -t \\sum_{ij,\\sigma=\\uparrow\\downarrow} {a_{i\\sigma}^\\dagger a_{j\\sigma}} + U \\sum_i {a_{i\\uparrow}^\\dagger a_{i\\uparrow} a_{i\\downarrow}^\\dagger a_{i\\downarrow}}\n",
    "\\end{equation*}</p>\n",
    "where:\n",
    "<ul>\n",
    "    <li>$a_{i\\sigma}^\\dagger$ = creation operator for site 'i' and spin $\\sigma$</li>\n",
    "    <li>$a_{i\\sigma}$ = annihilation operator for site 'i' and spin $\\sigma$</li>\n",
    "    <li>U = interaction term</li>\n",
    "    <li>t = hopping term</li>\n",
    "</ul>\n",
    "<p>The first sum quantifies moving electrons from each site 'j' to each neighboring site 'i' (first term, that is: $a_i^\\dagger$ \"creates\" an electron at site 'i' while $a_j$ simultaneously \"annihilates\" it from site 'j'), and vice versa (second term), considering both up and down spins ($\\sigma$) independently; the second sum quantifies the interaction contribution for having two electrons, one up and one down, at the same site 'i'.\n",
    "</p></div></td>\n",
    "    <td>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;</td>\n",
    "<td style=\"width: 35%\">\n",
    "    <img src='./images/basic_fermi_hubbard_model.png'>\n",
    "</td>\n",
    "</tr></table>\n",
    "\n",
    "<p>\n",
    "Although the pictorial shows the electrons as concrete particles, the model treats them as waves that extend over all sites. A solution to the model provides the probabilities with which any of the electrons can be measured, as a particle, at each one of the sites."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9feb427d",
   "metadata": {},
   "source": [
    "Mapping fermions to qubits\n",
    "---------------------------------\n",
    "\n",
    "Fermion algebra is only a subset of what a quantum computer can express, since fermions are interchangeable and many states are thus completely equivalent, whereas qubits are identifiable (each one has a well defined location on a chip or ion trap, for example). Thus, it is useful to perform a mapping from fermion to qubit algebra that preserves the fermionic relations, to ensure that the mapped operators only produce states inside that fermionic subset; meaning that the computation should not, noise aside, produce non-physical states. Doing so reduces the scope for errors.$^3$\n",
    "\n",
    "There are several possible mappings from fermion to qubit operators that preserve the fermion algebra. Here, we'll consider the Jordan-Wigner$^4$ and Bravyi-Kitaev$^5$ transforms, because these are commonly used and easy to understand and their differences$^6$ show an important point about how the optimizability of circuits.\n",
    "\n",
    "To get started, first load some relevant modules and helper functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4b7e427f",
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "50faf1cf",
   "metadata": {},
   "outputs": [],
   "source": [
    "# useful tools\n",
    "import logging\n",
    "import numpy as np\n",
    "import qiskit as qk\n",
    "import scipy\n",
    "from timeit import default_timer as timer\n",
    "\n",
    "# BQSKit tooling\n",
    "import bqskit\n",
    "from bqskit.ext import qiskit_to_bqskit\n",
    "\n",
    "# VQE objective function for Hamiltonians\n",
    "from hubbard import EnergyObjective, hamiltonian_qiskit, hubbard"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ee84a797",
   "metadata": {},
   "outputs": [],
   "source": [
    "def display_circuit(circuit, basis=['u', 'cx']):\n",
    "    \"\"\"perform basic, unoptimized decomposition for circuit display\"\"\"\n",
    "\n",
    "    uc =  qk.transpile(circuit, basis_gates=basis)\n",
    "    \n",
    "    for gate, count in uc.count_ops().items():\n",
    "        print(\"gate %s count: %d\" % (gate, count))\n",
    "\n",
    "    return uc.draw()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "402444d6",
   "metadata": {},
   "source": [
    "The Fermi-Hubbard model is defined by the dimensions of the 2D grid, plus the hopping and interaction terms. In addition, OpenFermion allows some variations on the team by making the electrons spinless and making the grid periodic. Below is a very small model, just 2 sites (2x1), which will result in 2 classical optimization parameters. You can try a larger model, but this small example is chosen for a tutorial session where it can be easily run on a laptop on battery power."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3fb9a577",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Setup for a 2x1 grid with a single electron.\n",
    "hm_lattice = {\n",
    "    'x_dimension' : 2,\n",
    "    'y_dimension' : 1,      # 2D x_dim * y_dim grid\n",
    "    'tunneling'   : 1.0,    # hpping term `t`\n",
    "    'coulomb'     : 2.0,    # interaction term `U`\n",
    "    'spinless'    : False,  # double negative, so do have spins\n",
    "    'periodic'    : False,  # do not allow wrap-around hopping\n",
    "}\n",
    "\n",
    "hm_filling = {\n",
    "    'n_electrons_up'   : 1,\n",
    "    'n_electrons_down' : 0,\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d4ce2011",
   "metadata": {},
   "source": [
    "Energy Objective\n",
    "--------------------\n",
    "\n",
    "The Variational Quantum Eigensolver is an iterative algorithm where a classical optimizer minimizes an energy objective (the expectation value of the Hamiltonian) as calculated on the quantum computer. Typically, an Hamiltonian has many components, which can be evaluated in parallel, on separate quantum chips/computers. The final result is obtained by summing the parts classically using their proper weights.\n",
    "<table style=\"padding: 0px\"><tr></tr><tr>\n",
    "<td style=\"width: 60%; padding: 0px\">\n",
    "<div align=\"left\" style=\"text-align: left; font-size: 120%\">\n",
    "The trial quantum states proposed by the classical optimizer are constructed using an <i>Ansatz</i> that is parametrized using classical values only.\n",
    "Here, we wil use <i>Unitary Coupled Cluster</i> Ansatz for single and double excitations (UCCSD). Briefly, this Ansatz explores the search space by moving single electrons and by swapping pairs.\n",
    "This Ansatz is easier on the classical optimizer than many alternatives, but produces rather deep circuits, which is why the use of a tool like BQSKit becomes a necessity.\n",
    "</div></td>\n",
    "<td><img src='./images/vqe_structure.png'></td>\n",
    "</tr></table>\n",
    "To complete a circuit, we need to fill it with electrons and then instantiate the circuit for a set of optimization paremeters. Note that the meaning of these parameters depends on the chosen mapping."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a3223436",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Hubbard Model Hamiltonian using Jordan-Wigner mapping\n",
    "hm_jw = hamiltonian_qiskit(fermion_transform='jordan-wigner', **hm_lattice)\n",
    "\n",
    "# create an energy objective, which can be minimized classically, from the mapped hamiltonian\n",
    "jw_objective = EnergyObjective(hm_jw, **hm_filling)\n",
    "\n",
    "# instantiate a circuit from a set of classical parameters\n",
    "jw_circuit_instance = jw_objective.generate_circuit([0.1]*jw_objective.npar())\n",
    "\n",
    "# visualize the circuit, using basic gates, and collect some statistics\n",
    "c_jw = qk.transpile(jw_circuit_instance, basis_gates=[\"u3\", \"cx\"])\n",
    "\n",
    "for gate, count in c_jw.count_ops().items():\n",
    "    print(\"gate %s count: %d\" % (gate, count))\n",
    "    \n",
    "c_jw.draw()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c678faab",
   "metadata": {},
   "source": [
    "**Exercise:** Instantiate an energy objective like before, but now for the Bravyi-Kitaev mapping. Draw and compare the generated circuits, find the main differences in structure, and try to guess how that might affect the synthesis, in particular based on what you have learned about partitioning.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c16f51be",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Hubbard Model Hamiltonian using Bravyi-Kitaev mapping\n",
    "hm_bk = hamiltonian_qiskit(fermion_transform='bravyi-kitaev', **hm_lattice)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dd8975b8",
   "metadata": {},
   "source": [
    "Effect of Fermion to Qubit Mapping\n",
    "-----------------------------------------\n",
    "\n",
    "The mapping has to encode the positions of the electrons and do so with a guarantee of anti-symmetry of the fermionic wave functions, to respect fermion algebra. There are thus two pieces of information: the electron occupancy on each site and the overall parity. The Jordan-Wigner mapping fully localizes the occupation number (each qubit represents the occupation on a site), which then causes the parity to be completely non-localized: whenever an electron \"hops\" sites, the parity update requires an update to all qubits, propagated through a chain of controlled gates. The Bravyi-Kitaev transform encodes both occupancy and parity information non-locally, leading to a more balanced result.\n",
    "\n",
    "**Exercise:** Check the circuits again and find that the Jordan-Wigner mapping results in many more series of CNOTs involving all qubits than does Bravyi-Kitaev. These series of CNOTs across all qubits are a worst-case scenario for synthesis, as they prevent further partitioning and are already optimal.\n",
    "\n",
    "<div class=\"alert alert-block alert-info\">\n",
    "<b>Note:</b> the circuits above were generated using the same number of Trotter steps (a method of approximating an exponential) and although the Bravyi-Kitaev circuit is significantly shorter in that case, the approximation is actually slightly worse. After adjusting the number of steps and the method to achieve the same quality of approximation (use <tt>trotter_steps=2, trotter_mode='trotter'</tt> for Jordan-Wigner and <tt>trotter_steps=3, trotter_mode='trotter'</tt>), both methods are much closer in the number of CNOTs in the unoptimized circuit. Since Trotter steps are just repetitions, this does not change any conclusions about the relative effectiveness of BQSKit's circuit synthesis on the respective mappings.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d54235b4",
   "metadata": {},
   "source": [
    "**Exercise:** Use what you learned in the __[BQSKit introductory notebook](1_comping_with_bqskit.ipynb)__ to compile both circuit instances. Try for `optimization_level=2` and `max_synthesis_size=3` first, adjust and play with the parameters if time allows (do note that larger block size and/or higher optimization level can result in significantly longer compilation times). Compare and explain the results.\n",
    "\n",
    "**Optional exercise:** The examples so far have presumed all-to-all. Experiment with one of the realistic QPU models, or your own. Are the results and conclusions the same?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "edad8c68",
   "metadata": {},
   "outputs": [],
   "source": [
    "opt_circuit = ...\n",
    "\n",
    "# Print new circuit statistics\n",
    "print(\"Compiled Circuit Statistics:\")\n",
    "print(\"  Gate Counts:\", opt_circuit.gate_counts)\n",
    "print(\"  Connectivity:\", opt_circuit.coupling_graph)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "db3f8861",
   "metadata": {},
   "source": [
    "Variational Quantum Eigensolver\n",
    "--------------------------------------\n",
    "\n",
    "In the VQE algorithm, circuits are synthesized based on parameters provided by a classical minimizer. All we need to do, is to pass the energy objective function to a classical optimizer. Since we are not considering noise here (see the [scikit-quant tutorials](https://github.com/scikit-quant/scikit-quant/tree/master/tutorials), if interested), most any optimizer will do. Here, we will use [BFGS from SciPy](https://docs.scipy.org/doc/scipy/reference/optimize.html#module-scipy.optimize)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f9850094",
   "metadata": {},
   "outputs": [],
   "source": [
    "min_energy = scipy.optimize.minimize(jw_objective, [0.1]*jw_objective.npar(), method='BFGS')\n",
    "print(\"minimum energy: %.5f (exact: %.5f)\" % (min_energy.fun, hubbard.exact(hm_jw, **hm_filling)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dda91c89",
   "metadata": {},
   "source": [
    "**Exercise:** Verify the results for the objective using the Bravyi-Kitaev mapping. Optionally, try different fillings and grid sizes (note that larger grids will have local minima where the optimizer can get stuck if not seeded close to the global minimum)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c68cbe8b",
   "metadata": {},
   "outputs": [],
   "source": [
    "min_energy = ...\n",
    "print(\"minimum energy: %.5f (exact: %.5f)\" % (min_energy.fun, hubbard.exact(hm_bk, **hm_filling)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "42ac5a2d",
   "metadata": {},
   "source": [
    "The Ansatz forms the bulk of any VQE circuit: both the input state and the Hamiltonian components only require single qubit rotations, at most, at the beginning and at the very end to measure the desired projection components for each qubit.\n",
    "It thus makes sense to optimize the Ansatz part of the circuit independently and re-use it (with at most a fast peep-hole pass to merge single qubit gates at start and end where possible).\n",
    "\n",
    "**Exercise:** The `EnergyObjective` is written such that a derived class can implement the `bqskit` optimization. Fill in the remainder below, and verify the results. Use the added pass to collect statistics on the `bqskit` results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7c6fbc11",
   "metadata": {},
   "outputs": [],
   "source": [
    "from hubbard import EnergyObjective\n",
    "\n",
    "# to show BQSKit is running\n",
    "bqskit.enable_logging()\n",
    "\n",
    "class MyObjective(EnergyObjective):\n",
    "    def optimize_evolution(self, bqskit_circuit):\n",
    "        opt_circuit = ... # use bqskit to improve the circuit\n",
    "        return opt_circuit\n",
    "    \n",
    "bk_objective = MyObjective(hm_bk, **hm_filling)\n",
    "start = timer()\n",
    "min_energy = scipy.optimize.minimize(bk_objective, [0.1]*bk_objective.npar(), method='BFGS')\n",
    "print(\"converged in %.0f seconds\" % (timer() - start))\n",
    "print(\"minimum energy: %.5f (exact: %.5f)\" % (min_energy.fun, hubbard.exact(hm_bk, **hm_filling)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b6a9f7e6",
   "metadata": {},
   "source": [
    "## Customizing Compilation for Better and Faster Workflows\n",
    "-----------------------------\n",
    "\n",
    "It is often useful to customize the compilation workflow that is set up by the `compile()` function, either because there are specific BQSKit tools you want to use, because you want to fine-tune settings, or simply because you want to add diagnostics to understand the synthesis performance better.\n",
    "In BQSKit, a [`Workflow`](https://bqskit.readthedocs.io/en/latest/source/autogen/bqskit.compiler.Workflow.html#bqskit.compiler.Workflow) object consists of a list of compilation passes executed in order. You can compile a circuit with a `Workflow` using a [`Compiler`](https://bqskit.readthedocs.io/en/latest/source/autogen/bqskit.compiler.Compiler.html#bqskit.compiler.Compiler) object, which handles the parallel execution environment. The off-the-shelf `compile()` function builds and executes a default workflow transparently. Let's take a look at this default workflow:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6dab19b0",
   "metadata": {},
   "outputs": [],
   "source": [
    "from bqskit.compiler.compile import build_workflow\n",
    "\n",
    "default_opt1_workflow = build_workflow(qiskit_to_bqskit(c_jw), optimization_level=1)\n",
    "print(default_opt1_workflow)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0bd1eff0",
   "metadata": {},
   "source": [
    "**Exercise:** Higher optimization workflows will, by definition, include more optimizations. Print out the workflows for higher optimization levels and identify the differences in the space below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "de459628",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(build_workflow(qiskit_to_bqskit(c_jw), ...)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8e082030",
   "metadata": {},
   "source": [
    "These default workflows are multi-faceted because they must deal with arbitrary input and work well in any situation.\n",
    "This necessity is not there when designing custom workflows as we often know great detail about our problem that we can exploit,\n",
    "leading to more straightforward, faster, and often higher-quality compilation workflows.\n",
    "\n",
    "The most common paradigm in BQSKit is to:\n",
    "\n",
    "1. Partition a circuit\n",
    "2. Do something on each partition\n",
    "3. Unpartition (unfold) the circuit\n",
    "\n",
    "For example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "20485ecf",
   "metadata": {},
   "outputs": [],
   "source": [
    "from bqskit.compiler import Workflow\n",
    "from bqskit.passes import QuickPartitioner, ForEachBlockPass, ScanningGateRemovalPass, UnfoldPass\n",
    "\n",
    "basic_gate_deletion_workflow = Workflow([\n",
    "    QuickPartitioner(3),  # Partition into 3-qubit blocks\n",
    "    ForEachBlockPass(ScanningGateRemovalPass()),  # Apply gate deletion to each block (in parallel)\n",
    "    UnfoldPass(),  # Unfold the blocks back into the original circuit\n",
    "])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "df832baa",
   "metadata": {},
   "source": [
    "The available passes are listed in the [documentation](https://bqskit.readthedocs.io/en/latest/source/passes.html). There are passes for partitioning, synthesis, conversions, iterative, conditionals, checkpointing, etc. etc. Of course, you can also write your own custom passes.\n",
    "\n",
    "Once you have a `Workflow`, you can use the BQSKIt compiler to execute it on a circuit:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9121c092",
   "metadata": {},
   "outputs": [],
   "source": [
    "from bqskit.compiler import Compiler\n",
    "\n",
    "with Compiler() as compiler:\n",
    "    opt_circuit = compiler.compile(qiskit_to_bqskit(c_jw), workflow=basic_gate_deletion_workflow)\n",
    "\n",
    "print(opt_circuit.gate_counts)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "00d4a6d3",
   "metadata": {},
   "source": [
    "**Exercise:** The `compile()` function uses the LEAP algorithm by default. Below are some example workflows to get started trying different partitioners and algorithms. Make sure to look up each pass in the [documentation](https://bqskit.readthedocs.io/en/latest/source/passes.html), and don't hesitate to try your own mix. Try them either with one of the earlier single circuits, or implement them in the `MyObjective` function from the previous exercise to use with a complete VQE run.\n",
    "\n",
    "```\n",
    "from bqskit import MachineModel\n",
    "from bqskit.passes.mapping.setmodel import SetModelPass\n",
    "from bqskit.passes import *\n",
    "\n",
    "workflow = [\n",
    "    SetModelPass(model),\n",
    "    QuickPartitioner(),\n",
    "    ForEachBlockPass([\n",
    "        QSearchSynthesisPass(),\n",
    "        ScanningGateRemovalPass()\n",
    "    ]),\n",
    "    UnfoldPass()\n",
    "]\n",
    "```\n",
    "\n",
    "**Optional Exercise:** To improve your understanding, here are two possible diagnostic passes that can be inserted into the workflow to show/debug the effect of individual passes. Add them into your workflow as needed and/or create your own:\n",
    "\n",
    "```\n",
    "from bqskit.compiler.basepass import BasePass\n",
    "from bqskit.ir.gates import CircuitGate, CNOTGate\n",
    "\n",
    "class PrintCNOTsPass(BasePass):\n",
    "    async def run(self, circuit, data) -> None:\n",
    "        print(\"BQSKit step, current CNOT count:\", circuit.count(CNOTGate()))\n",
    "        \n",
    "class PrintPartitionInfoPass(BasePass):\n",
    "    async def run(self, circuit, data):\n",
    "        acm = 0\n",
    "        for op in circuit:            \n",
    "            if not isinstance(op.gate, CircuitGate):\n",
    "                print(\"Circuit not yet partitioned\")\n",
    "                return\n",
    "\n",
    "            acm += op.gate._circuit.num_operations\n",
    "\n",
    "        print(\"Number of partitions:\", circuit.num_operations)\n",
    "        print(\"Average block size:\", acm / circuit.num_operations)\n",
    "```\n",
    "<div class=\"alert alert-warning\">\n",
    "<b>Important:</b> Any custom pass or object passed into a pass must be pickleable. Therefore, any custom pass must be defined in a separate importable module that isn't `__main__`.\n",
    "</div>\n",
    "\n",
    "<div class=\"alert alert-block alert-info\">\n",
    "<b>Note:</b> The printing passes will be run in a runtime worker process, and jupyter will therefore not pick up the output. You will find the printed output in the terminal where jupyter was started. However, you can use the Python logging module normally, and the runtime will ship any logs with `bqskit` as a parent logger to the original endpoint. Granted, logging has to be enabled via `bqskit.enable_logging`.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6f5ffc29",
   "metadata": {},
   "outputs": [],
   "source": [
    "from bqskit.compiler.compiler import Compiler\n",
    "from bqskit import MachineModel\n",
    "from bqskit.passes import *\n",
    "from bqskit import enable_logging\n",
    "\n",
    "enable_logging()\n",
    "\n",
    "coupling_graph = [\n",
    "    (0, 1), (0, 2), (0, 3), (1, 2), (1, 3), (2, 3)\n",
    "]\n",
    "model = ...\n",
    "\n",
    "workflow = [\n",
    "    ...\n",
    "]\n",
    "\n",
    "with Compiler() as compiler:\n",
    "    opt_circuit = compiler.compile(qiskit_to_bqskit(c_bk), workflow)\n",
    "    \n",
    "print(\"Gate Counts:\", opt_circuit.gate_counts)\n",
    "print(\"Connectivity:\", opt_circuit.coupling_graph)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "24a16001",
   "metadata": {},
   "source": [
    "Seeded Synthesis for Improving Variational Workflows\n",
    "----------------------------------------------------\n",
    "\n",
    "This section discusses how we can leverage knowledge about our specific problem, a variational Hubbard model simulation, to build a compilation workflow that fits its needs much better than BQSKit off-the-shelf.\n",
    "In a variational setting, we run many circuits sequentially, modifying them slightly between each iteration. The time cost between each iteration is significant in an experimental environment as QPU access time is not cheap for two main reasons.\n",
    "First, QPU time may cost actual currency and be very expensive, so we want to minimize the time idling waiting to compile a new circuit. Second, while idling the QPU may be drifting causing a type of noise that may be hard to handle in the algorithm.\n",
    "Therefore, a long-running compilation workflow may not be advantageous in a variational setting.\n",
    "\n",
    "This being said, there are two positive observations that we can make from this environment that we can use to guide the creation of a custom workflow:\n",
    "\n",
    "1. We have practically unlimited time before accessing the QPU to perform analysis and compilation on the first circuit in sequence.\n",
    "2. The circuits are not likely to change much between each iteration.\n",
    "\n",
    "The first observation allows us to perform deep compilation and analysis routines on the circuit before starting the algorithm.\n",
    "In contrast, the second point enables us to \"seed\" a compilation workflow from the results of a previous iteration.\n",
    "\n",
    "Over the rest of this section, we will first learn how to do compilation seed analysis, then combine that with first-compilation detection to build a variational compilation workflow that greatly improves compile time without hurting compile quality.\n",
    "\n",
    "But before we start, let's get a baseline:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eeeb0e30",
   "metadata": {},
   "outputs": [],
   "source": [
    "from hubbard import EnergyObjective\n",
    "from bqskit.ir.gates import CNOTGate\n",
    "from bqskit.passes.partitioning import QuickPartitioner\n",
    "from bqskit.passes.synthesis.qsearch import QSearchSynthesisPass\n",
    "from bqskit.passes.processing import ScanningGateRemovalPass \n",
    "from bqskit.passes.control.foreach import ForEachBlockPass\n",
    "from bqskit.compiler import Compiler\n",
    "from bqskit.compiler.compile import get_instantiate_options\n",
    "\n",
    "# To record how well we are doing\n",
    "compiled_circuit_cnot_counts = []\n",
    "instantiate_options = get_instantiate_options(optimization_level=1)\n",
    "\n",
    "workflow = [\n",
    "    QuickPartitioner(3), \n",
    "    ForEachBlockPass(\n",
    "        QSearchSynthesisPass(instantiate_options=instantiate_options),\n",
    "        ScanningGateRemovalPass(instantiate_options=instantiate_options)\n",
    "    ),\n",
    "    UnfoldPass(),\n",
    "]\n",
    "\n",
    "# Start a compiler\n",
    "compiler = Compiler()\n",
    "\n",
    "class BaselineObjective(EnergyObjective):\n",
    "    def optimize_evolution(self, bqskit_circuit):\n",
    "        compiled_circuit = compiler.compile(bqskit_circuit, workflow)\n",
    "        compiled_circuit_cnot_counts.append(compiled_circuit.count(CNOTGate()))\n",
    "        return compiled_circuit\n",
    "    \n",
    "bk_objective = BaselineObjective(hm_bk, **hm_filling)\n",
    "start = timer()\n",
    "min_energy = scipy.optimize.minimize(bk_objective, [0.1]*bk_objective.npar(), method='BFGS')\n",
    "print(\"converged in %.0f seconds\" % (timer() - start))\n",
    "print(\"minimum energy: %.5f (exact: %.5f)\" % (min_energy.fun, hubbard.exact(hm_bk, **hm_filling)))\n",
    "print(\"Average CNOT gate count: %.2f\" % np.mean(compiled_circuit_cnot_counts))\n",
    "\n",
    "# Close the compiler\n",
    "compiler.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1728b944",
   "metadata": {},
   "source": [
    "Now, the goal of this section is to improve upon that runtime while not harming the CNOT quality.\n",
    "\n",
    "We are going to do this using a form of seeded synthesis. For a detailed walkthrough of this concept, refer to the [QSeed project](https://github.com/mtweiden/qseed).\n",
    "The general idea is that compiling similar circuits will be discover similar structures of our block-level optimizations, and by recognizing these, we can skip the block compilations altogether.\n",
    "\n",
    "Below are some helper functions that plot histograms of found structures and draw specified circuits.\n",
    "\n",
    "Structures from runs of synthesis are stored in the `AnalysisObjective.structures` list. These structures can be drawn using the `draw_structure` function defined below. `plot_structure_histograms` should take as input a list of structure circuits. It will show you which structures are most common in an optimized circuit. The structure associated with a label in the plotted histogram can be retrieved using `label_to_structure`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b87edac5",
   "metadata": {},
   "outputs": [],
   "source": [
    "from bqskit import Circuit\n",
    "from bqskit.ir.gates import CNOTGate, U3Gate\n",
    "import matplotlib.pyplot as plt\n",
    "from bqskit.ir.lang.qasm2.qasm2 import OPENQASM2Language\n",
    "from qiskit import QuantumCircuit\n",
    "\n",
    "\n",
    "def plot_structure_histograms(structures_list):\n",
    "    unique_structures = sorted(list(set(structures_list.copy())), key=len)\n",
    "    counts = [structures_list.count(s) for s in unique_structures]\n",
    "    labels = [l for (l, _) in enumerate(unique_structures)]\n",
    "    plt.figure(figsize=(12, 4))\n",
    "    plt.bar(labels, counts)\n",
    "    plt.ylabel('Structure Frequency')\n",
    "    plt.xlabel('Structure Labels')\n",
    "    plt.xticks(labels)\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "def draw_structure(structure):\n",
    "    circuit = structure_to_circuit(structure, 3)\n",
    "    lang = OPENQASM2Language()\n",
    "    qasm = lang.encode(circuit)\n",
    "    qk_circ = QuantumCircuit.from_qasm_str(qasm)\n",
    "    print(qk_circ)\n",
    "\n",
    "\n",
    "def label_to_structure(structures_list, label):\n",
    "    unique_structures = sorted(list(set(structures_list.copy())), key=len)\n",
    "    return unique_structures[label]\n",
    "\n",
    "\n",
    "def structure_to_circuit(structure, num_qudits = 3):\n",
    "    circuit = Circuit(num_qudits)\n",
    "    for q in range(num_qudits):\n",
    "        circuit.append_gate(U3Gate(), [q])\n",
    "    for (q1, q2) in structure:\n",
    "        circuit.append_gate(CNOTGate(), [q1, q2])\n",
    "        circuit.append_gate(U3Gate(), [q1])\n",
    "        circuit.append_gate(U3Gate(), [q2])\n",
    "    return circuit\n",
    "\n",
    "\n",
    "def count_cnots(partitioned_circuit: Circuit):\n",
    "    circuit_copy = partitioned_circuit.copy()\n",
    "    circuit_copy.unfold_all()\n",
    "    return circuit_copy.gate_counts[CNOTGate()]\n",
    "\n",
    "\n",
    "def get_structures(partitioned_circuit):\n",
    "    circuit_copy = partitioned_circuit.copy()\n",
    "    structures = []\n",
    "    for block in circuit_copy:\n",
    "        subcirc = Circuit.from_operation(block)\n",
    "        subcirc.unfold_all()\n",
    "        structures.append(tuple([\n",
    "            tuple(op.location) for op in subcirc \n",
    "            if op.num_qudits == 2\n",
    "        ]))\n",
    "    return structures"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b7ecc1b6",
   "metadata": {},
   "source": [
    "Let's analyze a few of the iterations during the objective function analysis for good seeds:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0a82d531",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Since we are \"off-the-QPU\" during analysis, we can take our time and\n",
    "# perform high-quality compilation. We will use a workflow that does\n",
    "# QSearch synthesis with high-quality settings followed by gate deletion.\n",
    "instantiate_options = get_instantiate_options(optimization_level=3)\n",
    "workflow = [\n",
    "    QuickPartitioner(3), \n",
    "    ForEachBlockPass(\n",
    "        QSearchSynthesisPass(instantiate_options=instantiate_options),\n",
    "        ScanningGateRemovalPass(instantiate_options=instantiate_options)\n",
    "    ),\n",
    "]\n",
    "\n",
    "# To save time during analysis, we only want to do compilation and\n",
    "# analysis once in a while during the algorithm. As you decrease this,\n",
    "# you will get better seeds, but it will take longer.\n",
    "analysis_period = 4\n",
    "\n",
    "class AnalysisObjective(EnergyObjective):\n",
    "    def __init__(self, *args, **kwargs):\n",
    "        super().__init__(*args, **kwargs)\n",
    "        self.structures = []\n",
    "        self.block_counter = 0\n",
    "\n",
    "    def optimize_evolution(self, bqskit_circuit):\n",
    "        if self.block_counter < analysis_period:\n",
    "            self.block_counter += 1\n",
    "            return bqskit_circuit\n",
    "        \n",
    "        self.block_counter = 0\n",
    "        with Compiler() as compiler:\n",
    "            opt_circuit = compiler.compile(bqskit_circuit, workflow)\n",
    "            self.structures.append(get_structures(opt_circuit))\n",
    "            opt_circuit.unfold_all()\n",
    "        return opt_circuit\n",
    "\n",
    "# We are only running this to get the structures, so we aren't printing anything here\n",
    "bottomup_objective = AnalysisObjective(hm_bk, **hm_filling)\n",
    "min_energy = scipy.optimize.minimize(bottomup_objective, [0.1]*bottomup_objective.npar(), method='BFGS')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e71d38fe",
   "metadata": {},
   "source": [
    "Now try profiling the circuits using `plot_structure_histogram`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "61d86042",
   "metadata": {},
   "outputs": [],
   "source": [
    "structures = []\n",
    "for structure_batch in bottomup_objective.structures:\n",
    "    structures.extend(structure_batch)\n",
    "plot_structure_histograms(structures)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "07f6013e",
   "metadata": {},
   "source": [
    "**Exercise:** Visualize seeds of interest by putting their labels in the `seed_labels` list below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f3846d0c",
   "metadata": {},
   "outputs": [],
   "source": [
    "seed_labels = [...]  # Put seeds of interest here!\n",
    "for label in seed_labels:\n",
    "    print(f'structure: {label}')\n",
    "    draw_structure(label_to_structure(structures, label))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0141130a",
   "metadata": {},
   "source": [
    "Let's put this together now that we have done our analysis. Ideally, we have selected our seeds to start synthesis with, but this might require some experimenting.\n",
    "\n",
    "**Exercise:** After profiling, select one or more seeds by entering the associated integer label into the `seed_labels` list below.\n",
    "Note that some less common seeds may still be worth investigating if they are similar to more frequent ones."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c11d70b8",
   "metadata": {},
   "outputs": [],
   "source": [
    "seed_labels = [...]\n",
    "\n",
    "assert len(seed_labels) > 0, \"Make sure to put some seeds in the list above!\"\n",
    "\n",
    "# Build the seed structures\n",
    "seeds = [structure_to_circuit(label_to_structure(structures, s)) for s in seed_labels]\n",
    "\n",
    "# Build the workflow\n",
    "layer_gen = SeedLayerGenerator(seeds)\n",
    "instantiate_options = get_instantiate_options(optimization_level=3)\n",
    "workflow = [\n",
    "    QuickPartitioner(3), \n",
    "    ForEachBlockPass(\n",
    "        QSearchSynthesisPass(\n",
    "            layer_generator=layer_gen,\n",
    "            instantiate_options=instantiate_options,\n",
    "        ),\n",
    "        ScanningGateRemovalPass(instantiate_options=instantiate_options)\n",
    "    ),\n",
    "    UnfoldPass()\n",
    "]\n",
    "\n",
    "# To record how well we are doing\n",
    "compiled_circuit_cnot_counts = []\n",
    "\n",
    "# Start a compiler\n",
    "compiler = Compiler()\n",
    "\n",
    "class SeededObjective(EnergyObjective):\n",
    "    def optimize_evolution(self, bqskit_circuit):\n",
    "        compiled_circuit = compiler.compile(bqskit_circuit, workflow)\n",
    "        compiled_circuit_cnot_counts.append(compiled_circuit.count(CNOTGate()))\n",
    "        return compiled_circuit\n",
    "\n",
    "bk_objective = SeededObjective(hm_bk, **hm_filling)\n",
    "start = timer()\n",
    "min_energy = scipy.optimize.minimize(bk_objective, [0.1]*bk_objective.npar(), method='BFGS')\n",
    "print(\"converged in %.0f seconds\" % (timer() - start))\n",
    "print(\"minimum energy: %.5f (exact: %.5f)\" % (min_energy.fun, hubbard.exact(hm_bk, **hm_filling)))\n",
    "print(\"Average CNOT gate count: %.2f\" % np.mean(compiled_circuit_cnot_counts))\n",
    "\n",
    "# Close the compiler\n",
    "compiler.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "07dbbb18",
   "metadata": {},
   "source": [
    "How much speed up did you observe by using seeds? How much was the final CNOT count impacted compared to the baseline workflow? How many seeds do you use to get these results?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5fdb64a2",
   "metadata": {},
   "source": [
    "More information about a machine learning approach to seed recommendation can be found in the [QSeed paper](https://arxiv.org/abs/2306.05622) [7]. Circuit seeding in bqskit is done with the use of the [SeedLayerGenerator](https://bqskit.readthedocs.io/en/latest/source/autogen/bqskit.passes.SeedLayerGenerator.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d2c2de68",
   "metadata": {},
   "source": [
    "Implications of Noise\n",
    "-------------------------\n",
    "\n",
    "The use of BQSKit truly comes into its own when there is noise present in the system, as is the case in any current, real, QPU: the circuit optimization will show real improvements in the final result. As a note of caution, however: when noise is present, several things change. For example, a noise-resilient optimizer is needed and even then, the final energy result obtained is likely to be higher than the global minimum, even when the final parameters are spot-on. Additionally, simulation with noise takes considerably longer than without, and multiple runs are needed for convincing statistics.\n",
    "\n",
    "Install scikit-quant, to make its optimizers available as a Qiskit optimizer:\n",
    "\n",
    "    `$ python -m pip install scikit-quant`\n",
    "\n",
    "Construct a simple noise model (the circuits as-is will, even after optimization, be too long to be run on current hardware; you can try that by setting `noise_model='realistic'`) and pass it through the constructor of the `EnergyObjective`. Run with the base `EnergyObjective` and with your derived implementation using BQSKit. Compare the results with and without `bqskit` optimization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2acae37e",
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiskit_algorithms.optimizers import IMFIL\n",
    "from hubbard import noise_model\n",
    "\n",
    "logging.getLogger('hubbard').setLevel(logging.INFO)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7097f04f",
   "metadata": {},
   "outputs": [],
   "source": [
    "noise = noise_model.create(p_gate_1q=0.01, p_gate_2q=0.02)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6c06f5eb",
   "metadata": {},
   "outputs": [],
   "source": [
    "objective = EnergyObjective(hm_bk, noise_model=noise, **hm_filling)\n",
    "npar = objective.npar()\n",
    "\n",
    "result = IMFIL(maxiter=100).minimize(objective, np.array([0.1]*npar), bounds=np.array([(-1,1)]*npar))\n",
    "print(\"minimum energy: %.5f (at solution %s: %.5f; exact: %.5f)\" % (\n",
    "    result.fun, result.x, EnergyObjective(hm_bk, **hm_filling)(result.x), hubbard.exact(hm_bk, **hm_filling)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "428eb4ae",
   "metadata": {},
   "outputs": [],
   "source": [
    "objective = MyObjective(hm_bk, noise_model=noise, **hm_filling)\n",
    "npar = objective.npar()\n",
    "\n",
    "result = IMFIL(maxiter=100).minimize(objective, np.array([0.1]*npar), bounds=np.array([(-1,1)]*npar))\n",
    "print(\"minimum energy: %.5f (at solution %s: %.5f; exact: %.5f)\" % (\n",
    "    result.fun, result.x, EnergyObjective(hm_bk, **hm_filling)(result.x), hubbard.exact(hm_bk, **hm_filling)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5f07163a",
   "metadata": {},
   "source": [
    "## References\n",
    "-------------\n",
    "\n",
    "[1] J. Hubbard, Proc. Roy. Soc. (Lond.) 276 (1963) 238\n",
    "<br>[2] https://github.com/scikit-quant/scikit-quant/tree/master/tutorials\n",
    "<br>[3] https://arxiv.org/abs/cond-mat/0012334, Ortiz et al., \"Quantum Algorithms for Fermionic Simulations\"\n",
    "<br>[4] P. Jordan and E. Wigner, \"Über das Paulische Äquivalenzverbot\", Z. Phys. 47, 631 (1928)\n",
    "<br>[5] https://arxiv.org/abs/quant-ph/0003137 S. Bravyi, A. Kitaev, \"Fermionic quantum computation\"\n",
    "<br>[6] https://arxiv.org/abs/1812.02233, A. Tranter, et al., \"A comparison of the Bravyi-Kitaev and Jordan-Wigner transformations for the quantum simulation of quantum chemistry\"\n",
    "<br>[7] https://arxiv.org/abs/2306.05622, M. Weiden, et al., \"Improving Quantum Circuit Synthesis with Machine Learning\""
   ]
  }
 ],
 "metadata": {
  "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.10.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}