Lab_UCLA_blank.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "67d8f6e3-705f-4ea8-aff7-ca7dfe263579",
   "metadata": {},
   "source": [
    "# Lab: Time Evolution and Variational Algorithm\n",
    "\n",
    "Wenyang Qian, April 26, 1 PM - 3 PM\n",
    "\n",
    "Qiskit Lab Section (two hours) for Quantum Computing Lectures at UCLA"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12cac13e-e4a8-4f0a-b738-d0f15de3d79a",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">\n",
    "\n",
    "**Readme**  \n",
    "    \n",
    "There are 17 exercises (marked in Green) in total in this Jupyter notebook. These questions will teach you on\n",
    "\n",
    "Part 0. Basics in Qiskit to perform time evolution and expectation measurement (Exercises 1-6)\n",
    "\n",
    "Part 1. Real-time evolution with the Ising model (Exercises 7-12)\n",
    "\n",
    "Part 2. Variational Quantum Eigensolvers from scratch to solve the Maxcut Problem (Exercises 13-17)\n",
    "    \n",
    "Most of these exercises are straightforward, with solutions in one or two lines. Basically, you need to fill the space between:\n",
    "\n",
    "\"\n",
    "##### ==================================\n",
    "\\# Write your solution in here.  \n",
    "\n",
    "\n",
    "\n",
    "##### ==================================\n",
    "\"\n",
    "\n",
    "This is not a test so be feel free to ask your mentors and discss with your classmates if you have questions.\n",
    "\n",
    "We really hope you enjoy learning quantum computing with Qiskit. If you have any question after the class, be free to email qian dot wenyang at usc dot es.\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f6723436-85ae-41ae-94de-ac54f683f983",
   "metadata": {},
   "outputs": [],
   "source": [
    "# make sure these libraries are install, if not, use pip to install them\n",
    "# took me 2-3 mins (Apr 22, 2024)\n",
    "\n",
    "# %%capture\n",
    "\n",
    "# !pip install qiskit==0.44\n",
    "# !pip install qiskit-aer==0.12\n",
    "# !pip install qiskit_algorithms==0.2\n",
    "# !pip install matplotlib\n",
    "# !pip install pylatexenc\n",
    "# !pip install numpy\n",
    "# !pip install networkx\n",
    "# !pip install scipy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a14c1347-f68f-416e-a885-1a8ec604e907",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import networkx as nx\n",
    "from time import time as current_time\n",
    "from scipy import linalg as LA\n",
    "\n",
    "from qiskit import execute, Aer, transpile\n",
    "from qiskit.circuit import QuantumCircuit\n",
    "from qiskit.circuit.library import PauliEvolutionGate, TwoLocal\n",
    "from qiskit.quantum_info import Pauli, SparsePauliOp, Statevector\n",
    "from qiskit.tools.visualization import plot_histogram, plot_distribution\n",
    "from qiskit.primitives import Estimator, Sampler\n",
    "from qiskit.synthesis import LieTrotter, SuzukiTrotter\n",
    "\n",
    "import qiskit_aer\n",
    "import qiskit\n",
    "import qiskit_algorithms\n",
    "import qiskit.tools.jupyter\n",
    "\n",
    "%qiskit_version_table"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "029d5082-2d4f-452b-931d-04cb927c4228",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Part 0. Preparation\n",
    "\n",
    "In this part, you will learn the basic tools to construct operators $\\hat{H}$, simulate quantum circuit $\\exp(-i\\hat{H}t)$ in real time, and evaluate observable such as $\\langle \\hat{O} \\rangle$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "387284fd-742c-441e-91f2-7c7578dcc70f",
   "metadata": {},
   "source": [
    "### 0.1. Pauli operations in Qiskit"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "155ef393-cd9d-4324-b426-351b10bc902b",
   "metadata": {},
   "source": [
    "In Qiskit, building the Pauli strings to represent the Hamiltonian is usually the first step. Typically, we can use the *SparsePauliOp* class. You can read more about *SparsePauliOp* class [here](https://qiskit.org/documentation/stubs/qiskit.quantum_info.SparsePauliOp.html)\n",
    "\n",
    "We start by writing the standard **Pauli matrices** in Qiskit. Make sure you check these single-qubit Pauli operators in matrix form and it should agree with your expectation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d4287265-50d7-4156-84bb-50812b7b0315",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "X = SparsePauliOp('X')\n",
    "Y = SparsePauliOp('Y')\n",
    "Z = SparsePauliOp('Z')\n",
    "I = SparsePauliOp('I')\n",
    "\n",
    "print(\"Pauli-X operator:\")\n",
    "print(X.to_matrix())\n",
    "print(\"\\nPauli-Y operator:\")\n",
    "print(Y.to_matrix())\n",
    "print(\"\\nPauli-Z operator:\")\n",
    "print(Z.to_matrix())\n",
    "print(\"\\nPauli-I operator:\")\n",
    "print(I.to_matrix())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e818d838-735e-4198-bd47-000699c3c912",
   "metadata": {},
   "source": [
    "These are single-qubit Pauli matrices, which are the building blocks of more complicated Pauli objects (typically over multiple qubits).  \n",
    "\n",
    "To build more interesting Pauli operators, we need the following the Pauli arithmetic operations.\n",
    "\n",
    "(a) **Multiple-qubit Pauli Operators**. You need represent tensor product of Paulis with $^$ symbol. See the examples below."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9b0c6ad3-0d13-421e-a728-5e3a08685bfe",
   "metadata": {},
   "source": [
    "Example 1: 4-qubit tensor product, $I \\otimes X\\otimes Y \\otimes Z$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b7677475-f858-4d8b-9342-08293a5066dd",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "I^X^Y^Z"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c01c475a-b51a-4016-8ad1-042ec98bd302",
   "metadata": {},
   "source": [
    "Note, Qiskit used a reversed notation, so it is actually $I_3 \\otimes X_2\\otimes Y_1 \\otimes Z_0$. (PS: Typically if you use a consistently, you can decide to whether to invert your measured outcomes at the end of simulation or not.)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "858f6cf8-976a-4471-b92f-3cdeeb166f67",
   "metadata": {},
   "source": [
    "(b) **Scalar Product**. You will need to use the $*$ symbol between the coefficient and the Pauli."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3a01325a-374e-4120-9b7b-8ead99d79082",
   "metadata": {},
   "source": [
    "Example 2: 4-qubit tensor product, $4 I \\otimes X\\otimes Y \\otimes Z$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fb249131-f189-4370-a585-57e514be065b",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "4*I^X^Y^Z"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3f18baa1-d384-47a8-aeeb-7d5bc32365dd",
   "metadata": {},
   "source": [
    "(c) **Linear Combinations of Pauli Operators**. You will use the $+$ symbol to write a linear combination of the Pauli operators. Only Pauli operators of the same number of qubits can be joined. Different components must bed grouped with parenthesis.\n",
    "\n",
    "Example 3: 4-qubit linear combinations, $X\\otimes X\\otimes I \\otimes I + 2I \\otimes Y\\otimes  Y\\otimes I + 3I \\otimes I\\otimes  Z\\otimes Z$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8546469d-7ef3-4de1-be53-1eadf6eb6216",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "(X^X^I^I) + 2*(I^Y^Y^I) + 3*(I^I^Z^Z)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "da78dcab-938b-4ebc-b3af-f0684c334b39",
   "metadata": {},
   "source": [
    "(d) **Product of Pauli Operators**. *SparsePauliOp* class really shines when you want to compute qubitwise product of pauli operators. You will use the $@$ symbol. \n",
    "\n",
    "In the example below, you can see it supports algebraic expression and knows Pauli matrix properties. The final result can be further simplified with *.simplify()* method. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "21ed15ce-d1ec-47d5-b306-2be218cd5f4d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# contructing an interesting Pauli operator\n",
    "\n",
    "op1 = X^Y\n",
    "op2 = Y^X\n",
    "\n",
    "res = op1@op1 + op2@op2 + op1@op2 + op2@op1\n",
    "\n",
    "print(res)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5a13a4cc-7e7c-47c8-8963-83daed8a75a0",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "print('\\nSum of product after .simplify()')\n",
    "print(res.simplify())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9d540fae-05e3-446f-a6b8-9300a079f1e5",
   "metadata": {},
   "source": [
    "### 0.2. The Swap Hamiltonian"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7b98726e-68ca-43ee-8a00-a6168bb77846",
   "metadata": {},
   "source": [
    "Now, in the exercise below, let us onstruct a swap Hamiltonian with what you learned about Pauli operations in Qiskit"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ed8f6ea1-e5ef-4feb-9b34-6f2b9f4b6309",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">\n",
    "\n",
    "**Exercise 1. Building Multi-qubit Pauli Operator**  \n",
    "    \n",
    "Build Pauli Operator $$H_{sw} = \\frac{1}{2}(I_0I_1+X_0X_1 + Y_0Y_1+ Z_0Z_1)$$ using your X, Y, Z, I variables defined above. Or you can use other options of *SparsePauliOp* directly its documentation.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4a233f69-182d-43ba-b59d-a3b22a10f471",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "##### ==================================\n",
    "# Write your solution in here.\n",
    "\n",
    "##### ==================================\n",
    "\n",
    "print(H_sw)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "29b63554-0493-4771-ba46-104ffb864ce5",
   "metadata": {},
   "source": [
    "Great, now you have build an important quantum operator! (You will see why in a few more steps.)\n",
    "\n",
    "Note that manual construction of the Pauli operators can be tedious, especially when involving many terms and multiple qubits. For example, what if we want to build the same operator for 10 qubits, instead of 2 qubits?\n",
    "\n",
    "Luckily, it can done easily with *SparsePauliOp* class using *from_sparse_list* method. In each tuple, one specify the Pauli types, qubit locations, and the Pauli coefficient."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "79a8b23a-19a6-450d-86e7-72fc05096903",
   "metadata": {},
   "outputs": [],
   "source": [
    "II_tuples = [(\"II\", [0, 1], 1/2)] # \n",
    "XX_tuples = [(\"XX\", [0, 1], 1/2)]\n",
    "YY_tuples = [(\"YY\", [0, 1], 1/2)]\n",
    "ZZ_tuples = [(\"ZZ\", [0, 1], 1/2)]\n",
    "L = 10\n",
    "\n",
    "SparsePauliOp.from_sparse_list([*II_tuples, *XX_tuples, *YY_tuples, *ZZ_tuples], num_qubits=L)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c0558b10-f222-4909-8d38-d6f51ef051dc",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">\n",
    "\n",
    "**Exercise 2. Building a Function for Multi-qubit Pauli Operator**  \n",
    "    \n",
    "Build a function that constructs the Pauli Operator\n",
    "$$H_{sw} = \\frac{1}{2}(I_0I_{N-1}+X_0X_{N-1} + Y_0Y_{N-1}+ Z_0Z_{N-1})$$ \n",
    "over $N$ qubits.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f14e97fe-3d63-41ea-9bbd-d2fce49f28fe",
   "metadata": {},
   "outputs": [],
   "source": [
    "def H_sw_operator(N):\n",
    "    ##### ==================================\n",
    "    # Write your solution in here.\n",
    "\n",
    "    \n",
    "    ##### ==================================\n",
    "    res = SparsePauliOp.from_sparse_list([*II_tuples, *XX_tuples, *YY_tuples, *ZZ_tuples], num_qubits=N)\n",
    "    return res.simplify()\n",
    "\n",
    "N = 2\n",
    "print(H_sw_operator(N))\n",
    "\n",
    "N = 6\n",
    "print(H_sw_operator(N))\n",
    "\n",
    "N = 10\n",
    "print(H_sw_operator(N))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cb4738fa-e0ab-4c08-b6fd-bd2f0629a36c",
   "metadata": {},
   "source": [
    "Please make sure what you obtain match your expectation. \n",
    "\n",
    "Now you should be confident to build any Pauli operator in Qiskit. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ce9c85e6-5b41-468c-8700-a321f7035b12",
   "metadata": {},
   "source": [
    "### 0.3. Quantum Circuit Time Evolution (Single Trotter Step)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e7b3155d-51ed-411c-b1d7-bf03139dbab6",
   "metadata": {},
   "source": [
    "In class, we learned how to simulate any Pauli-based exponential $\\exp(-iHt)$. Let us briefly recall,\n",
    "\n",
    "* Simplify $H$ into disjoint, commuting Pauli terms.\n",
    "* Ultilize different Trotter formulas, such as Lie-Trotter to linear order, to approximate the expression.\n",
    "* Build towers of Rz (Rotation Z gate) and CX (Control-Not gate) for purely Pauli-Z based operator\n",
    "* Prepend and Append H, or SH gates for Pauli-X and Pauli-Y based operators.\n",
    "* The coefficients naturally appear in the Rz rotation angle.\n",
    "\n",
    "Though these steps are straightforward, manual construction of a generic circuit is tedius and complicated. Luckily, with Qiskit, you do not worry about it. All these can be easily constructed using *PauliEvolutionGate* class. (there are also many more classes)\n",
    "\n",
    "See an example below, to construct a two-qubit $\\exp\\big[-i (XY) t\\big]$ for t = 1.5 seconds. Note, we used **First-order** Lie Trotter formula here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "97a4aca0-3001-43e8-8da8-2a68c76f5ebf",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "peg = PauliEvolutionGate(operator=(X^Y), time=1.5, synthesis=LieTrotter())\n",
    "peg"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "87ff6c6f-0f7e-4f9d-9c65-fbe2543ff514",
   "metadata": {},
   "source": [
    "Note *PauliEvolutionGate* is essentially quantum circuit instructions. To build it into a Qiskit quantum circuit, you will need to append these instruction to some quantum circuits, for example below"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f50c7e48-5780-47b8-97d7-9b8bfe7f7c98",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "qc = QuantumCircuit(2) # our main circuit\n",
    "qc.append(peg, [0, 1]) \n",
    "# [0, 1] indicates the location of the appending operation. \n",
    "# Here, it is qubit-0 to qubit-0 and qubit-1 to qubit-1\n",
    "# as you can see in the circuit, 0 is connected to 0, 1 is connected to 1\n",
    "# usually for arbitrary length N, just do qc.append(peg, range(N))\n",
    "qc.draw(output='mpl')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "79f55654-486d-4597-9f97-84e2b7b9f96b",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">\n",
    "\n",
    "**Exercise 3. Building a Function to Construct a Quantum Circuit to Evolve a Pauli Operator**  \n",
    "    \n",
    "Build a generic function to construct a quantum circuit for\n",
    "\n",
    "$$\\exp(-i H_{operator} t)$$\n",
    "\n",
    "using *PauliEvolutionGate* class. \n",
    "\n",
    "Note, it must return a quantum circuit object with an input Pauli operator and an input time.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9f10f293-f3e7-459e-b733-3b236bfd07d2",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def time_evolve(H_operator, time, synthesis=LieTrotter()):\n",
    "    num_qubit = H_operator.num_qubits\n",
    "    ##### ==================================\n",
    "    # Write your solution in here.    \n",
    "\n",
    "\n",
    "    \n",
    "    ##### ==================================\n",
    "    return qc\n",
    "\n",
    "# Test your circuit\n",
    "N = 2\n",
    "H_operator = H_sw_operator(N)\n",
    "qc_test = time_evolve(H_operator, 1.5)\n",
    "qc_test.decompose().draw(output='mpl')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b520fe65-d3b9-4a7d-92b7-403d6b74b3a1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Test your circuit\n",
    "N = 4\n",
    "H_operator = H_sw_operator(N)\n",
    "qc_test = time_evolve(H_operator, 1.5)\n",
    "qc_test.decompose().draw(output='mpl')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "03a365a2-d71f-48ac-be8f-091294664f27",
   "metadata": {},
   "source": [
    "You should see 4 Pauli exponential blocks here. Why? Because all these Pauli terms are commutative. So, we can simulate them blockwise and without error!"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e832bd97-f6f8-4957-b7a6-efd7214f8c1a",
   "metadata": {},
   "source": [
    "### 0.4. Quantum Circuit Transpilation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "415308a5-9955-46fe-9896-5cce73fef277",
   "metadata": {},
   "source": [
    "You may ask, how do we verify them individually and see the tower structure of CX gates and RZ gates we saw in class?\n",
    "\n",
    "There are more than one ways to decompose the quantum circuit. To see what we saw in class, we need to transpile the quantum circuits to the basis gates that we are assigning them.\n",
    "\n",
    "Introducing the *transpile* function, which takes a quantum circuit and a list of basis gates. Run the command below and verify with your expectation."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ebc47ef8-7883-4c65-bebe-2cda5c70ecab",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">\n",
    "\n",
    "**Exercise 4. Transpiling Quantum Circuit to Specified Gate Sets**  \n",
    "    \n",
    "Transpile your quantum circuit to basis gates \"h\", \"s\", \"sdg\", \"rz\", \"cx\". Does it match your expectation?\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9c7e9b1e-c61f-4cd4-a9b1-5bd0de79733b",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "N = 2\n",
    "H_operator = H_sw_operator(N)\n",
    "qc_test = time_evolve(H_operator, 1.5).decompose()\n",
    "\n",
    "transpiled_qc = transpile(qc_test, basis_gates=['h', 's', 'sdg', 'rz', 'cx'], optimization_level=0) # no optimization = 0\n",
    "transpiled_qc.draw(output='mpl')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d718a1ac-fa12-4ced-b0b8-6f67ecf12f2e",
   "metadata": {},
   "source": [
    "Note, you should try to identify **three** towers of Rz gates with CXs. \n",
    "\n",
    "The $\\exp(-iIIt)=\\cos(t) II - i \\sin(t) II = \\Big(cos(t)-i\\sin(t)\\Big) II$ becomes part of the global phase, because  so Qiskit automatically skipped its simulation."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c0b81b4-1958-4597-858a-dc68f7bd0bdd",
   "metadata": {},
   "source": [
    "### 0.5. Quantum Circuit Initialization and Measurement"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c04bfb7f-5dc7-4853-a787-7f6f940cee53",
   "metadata": {},
   "source": [
    "Before simulating your quantum circuit, we need to add an initialization block and a measurement block.\n",
    "    \n",
    "To make things interesting, try initialize qubit-0 as $|1\\rangle$, qubit-1 as $|0\\rangle$ state, and then time evolve your two-qubit $H_{sw}$ for $\\pi/2$ period, and finally add a measurement block.\n",
    "\n",
    "Now we have a complete circuit!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "87946a94-855d-4f08-8df1-7f56519320cc",
   "metadata": {},
   "outputs": [],
   "source": [
    "qc_test = QuantumCircuit(2)\n",
    "\n",
    "# 1. Initialize qubit_0 as |1>, qubit_1 as |0> state, simply apply X gate on qubit 0\n",
    "qc_test.x(0)\n",
    "\n",
    "# 2. Time evolve for t = pi/2, using your function *time_evolve* \n",
    "H_operator = H_sw_operator(N=2)\n",
    "qc_test.append(time_evolve(H_operator, np.pi/2), range(2))\n",
    "\n",
    "# 3. Measure all\n",
    "qc_test.measure_all()\n",
    "\n",
    "qc_test.decompose().draw(output='mpl')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cb7ca788-91f5-4db7-81a6-46bd10e71dfa",
   "metadata": {},
   "source": [
    "To measure the outcome of any quantum circuit, we can use *Sampler* class in Qiskit. \n",
    "\n",
    "* With *shots=0*, we mimic an exact simulation using the state vector. \n",
    "\n",
    "* With *shots* = a positive integer, we mimic a practical simulation with however many shots on an **ideal** quantum computer. Obviously, the more shots one have, the closer to the exact simulation.\n",
    "\n",
    "For a realistic simulation on a current **Noisy** quantum computer, one could either introduce a noise model or use IBMQ cloud quantum device (for free!). I will leave this topic to another time given the limited time. [If you are interested, we can discuss later.]\n",
    "\n",
    "Qiskit has gone through several implementations back and forth, starting with the *execute* function to *Sampler* class. Here, I will simplify write an auxiliary function *sample_with_shots* for our class to use."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3f3034ea-6b1c-4369-b35d-81422f4053d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "def pad_int(k, n):\n",
    "    return \"0\"*(n-len(bin(k)[2:]))+bin(k)[2:]\n",
    "    \n",
    "def sample_with_shots(circuit, num_qubits, shots=None):\n",
    "    ''' Using Sampler class, based on Qiskit (v0.44)\n",
    "        Exact simulation when shots = None\n",
    "        Shot simulation when shots = integer\n",
    "    '''\n",
    "    shot_sampler = Sampler(options={\"shots\": shots})\n",
    "    dist = shot_sampler.run(circuit).result().quasi_dists[0]\n",
    "    dist = {pad_int(k, num_qubits): v for k, v in dist.items()}  # this is added because of Qiskit newly version does not like bitstring anymore\n",
    "    return dist"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d20f628a-db7a-47a8-b7b0-5ee34c7e9d21",
   "metadata": {},
   "source": [
    "It is important to compare exact simulation with shot simulations with our Hamiltonian"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2df49d09-9116-4aab-9c42-3dab299d0131",
   "metadata": {},
   "outputs": [],
   "source": [
    "dist_exact  = sample_with_shots(qc_test, 2, None)\n",
    "dist_10     = sample_with_shots(qc_test, 2, 10)\n",
    "dist_100    = sample_with_shots(qc_test, 2, 100)\n",
    "dist_1000   = sample_with_shots(qc_test, 2, 1000)\n",
    "dist_10000  = sample_with_shots(qc_test, 2, 10000)\n",
    "\n",
    "dist_list = [dist_10, dist_100, dist_1000, dist_10000, dist_exact]\n",
    "\n",
    "plot_distribution(dist_list, legend=[\"10 shots\", \"100 shots\", \"1000 shots\", \"10000 shots\", \"Exact\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d0cd6253-39fd-42b0-8d74-9d5a4ba64006",
   "metadata": {},
   "source": [
    "Since we started with $|q_1q_0\\rangle = |01\\rangle$ state, you will see the result in **any shot** is now a $|10\\rangle = |2\\rangle$ state. \n",
    "\n",
    "Yes, that is right. You have just implemented a SWAP gate. There are many implementations of SWAP gates. Here, you simulate a SWAP gate by $\\exp(-i\\frac{(II+XX+YY+ZZ)}{2}\\theta)$ with an angle $\\theta=\\pi/2$, which your simulation time here.\n",
    "\n",
    "Now, what will happen if you try to rotate only half the angle? What do you think will happen.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "644209fe-d781-4293-9217-3fe633a5c5ff",
   "metadata": {
    "tags": []
   },
   "source": [
    "<div class=\"alert alert-block alert-success\">\n",
    "\n",
    "**Exercise 5. Simulate with Ideal Shot-based Simulator**  \n",
    "    \n",
    "Now let us try on our own. What if we time evolve the H_sw_operator on $N=4$ qubits for an evolution time to $\\pi/4$. \n",
    "    \n",
    "Try to run different number of shots, 10, 100, 1000, 10000, and exact simulation. What do you observe? You may also try it several times.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a2c1c00b-310a-4274-a883-000aac694334",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "\n",
    "##### ==================================\n",
    "# Write your solution in here.\n",
    "qc_test = QuantumCircuit(4)\n",
    "\n",
    "# 1. Initialize qubit_0 as |1>, qubit_1 as |0> state\n",
    "\n",
    "\n",
    "# 2. Time evolve for t = pi/4\n",
    "\n",
    "\n",
    "# 3. Measure all\n",
    "\n",
    "\n",
    "##### ==================================\n",
    "\n",
    "dist_exact  = sample_with_shots(qc_test, 4, None)\n",
    "dist_10     = sample_with_shots(qc_test, 4, 10)\n",
    "dist_100    = sample_with_shots(qc_test, 4, 100)\n",
    "dist_1000   = sample_with_shots(qc_test, 4, 1000)\n",
    "dist_10000  = sample_with_shots(qc_test, 4, 10000)\n",
    "\n",
    "dist_list = [dist_10, dist_100, dist_1000, dist_10000, dist_exact]\n",
    "\n",
    "plot_distribution(dist_list, legend=[ \"10 shots\", \"100 shots\", \"1000 shots\", \"10000 shots\", \"Exact\",])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "75a85efc-1eb3-4705-bafe-e9d1f90137af",
   "metadata": {},
   "source": [
    "It should be an even superposition of $|0001\\rangle=|1\\rangle$ and $|1000\\rangle=|8\\rangle$ states right? With the number of shots going larger and larger, the probability distributions should be 0.5 for each.\n",
    "\n",
    "In practice, this is very important property, with different angles, you can generate desire weighted sum of these one-hot encoded bitstrings, which has many applications especially in quantum machine learning and combinatorial problems."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a587ac95-cfc7-4be3-b529-21cc35f27283",
   "metadata": {},
   "source": [
    "Now that, we have talked about simulating any Hamiltonian. We will switch to discussing extracting observables from the quantum circuit. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6dbd7025-1cd8-4fa6-b0f7-e9eec816aba3",
   "metadata": {},
   "source": [
    "### 0.6. Quantum Circuit Expectation Value of Observable"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6efe7447-da8f-4e90-87af-2ffc89c04f69",
   "metadata": {
    "tags": []
   },
   "source": [
    "More often or not, we are interesting in obtaining the expectation values of certain observable for our evolve state. Those information is useful and can be used to characterize the state of our quantum system, such as the energy.\n",
    "\n",
    "For example, in our previous simulation, we can also introduce an operator \n",
    "\n",
    "$$O=\\frac{1}{2}((I-Z)\\otimes I\\otimes I \\otimes I - I\\otimes I\\otimes I \\otimes (I-Z))$$\n",
    "\n",
    "to describe the asymmetry of the sampled outcome. Remember $\\frac{(I-Z)}{2}|0\\rangle = 0|0\\rangle=0$ and $\\frac{(I-Z)}{2}|1\\rangle = 1|1\\rangle$. Predict what you get if measuring this observable for the exact simulation distribution above?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "892f4e6f-ae38-47f8-9d46-ffdb172ec85e",
   "metadata": {
    "tags": []
   },
   "source": [
    "<div class=\"alert alert-block alert-success\">\n",
    "\n",
    "**Exercise 6. Building Pauli Observables**  \n",
    "    \n",
    "Building Pauli observables is the same as constructing Pauli operators. Let us build\n",
    "$$O=\\frac{1}{2}((I-Z)\\otimes I\\otimes I \\otimes I - I\\otimes I\\otimes I \\otimes (I-Z))$$\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "088fc3d9-702c-42ac-91bc-f897112ed130",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "##### ==================================\n",
    "# Write your solution in here.\n",
    "\n",
    "\n",
    "##### ==================================\n",
    "O_test"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6dab04bd-d785-49bc-8d3e-f2f32dcff246",
   "metadata": {},
   "source": [
    "To run the estimation, there are multiple methods in Qiskit. For convenience, I have written a method for you to estimate your observable using the *Estimator* class.\n",
    "\n",
    "Play with the simulation with different shots below, and convince yourself everything makes sense."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a0da1138-7e22-4e6f-bcc6-c5d8fb700c88",
   "metadata": {},
   "outputs": [],
   "source": [
    "def estimate_with_shots(circuit, observable, shots=None):\n",
    "    ''' Using Estimator class, based on Qiskit (v0.44)\n",
    "        Exact estimation when shots = None\n",
    "        Shot estimation when shots = integer\n",
    "    '''    \n",
    "    _circuit = circuit.copy()\n",
    "    _circuit.remove_final_measurements()\n",
    "    shot_estimator = Estimator(options={\"shots\": shots})\n",
    "    expectation = shot_estimator.run(_circuit, observable).result().values[0]\n",
    "    return expectation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "655c12f9-196d-4e96-aafc-ddd480fb34d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "shot_list = np.round(np.geomspace(10, 10000, 30))\n",
    "expectation_list = []\n",
    "\n",
    "for shots in shot_list:\n",
    "    expectation = estimate_with_shots(qc_test, O_test, shots)\n",
    "    expectation_list.append(expectation)\n",
    "    \n",
    "exact_est = estimate_with_shots(qc_test, O_test, None)\n",
    "plt.hlines(exact_est, 0, 10000, ls='solid', color='red', label='exact simulation')\n",
    "plt.scatter(shot_list, expectation_list, label='shot simulation')\n",
    "plt.ylabel(r'$\\langle O\\rangle$', fontsize=20)\n",
    "plt.xlabel('shots', fontsize=20)\n",
    "plt.xticks(fontsize=15); plt.yticks(fontsize=15)\n",
    "plt.ylim(-0.4, 0.4)\n",
    "plt.legend(fontsize=14)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "451b781a-16a1-4e3d-8074-099b03e2816b",
   "metadata": {},
   "source": [
    "You should observe that the expectation value approaches 0 as the number of shots increase. It is because, our observable is basically a number operator (or counting operator if you like). \n",
    "\n",
    "It counts the number of 1's in highest and the lowest qubit and calculate their difference. From the distribution earlier, you should expect zero, since we have an even distribution of $|0001\\rangle$ and $|1000\\rangle$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aef34d9d-f3b0-4dad-9309-bfeeb1275ae1",
   "metadata": {},
   "source": [
    "You should also see that having enough shots is crucial to simulation results. On the other hand, having too many shots will not help the simulation after certain point. It is essentially closely related the state space of the many-qubit systems. In quantum information science, the theoretical number of shots $\\sim 1/\\mathcal{O}(\\epsilon^2)$ where $\\epsilon$ is the accuracy."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2b2b1953-a790-405b-b447-832c3ba21785",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Part 1. Ising Model\n",
    "\n",
    "In this part, you will simulate the 1-dimensional Ising model Hamiltonian in real time using quantum circuit. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "77a1f657-8043-4c55-aff6-d69d4b9db10c",
   "metadata": {},
   "source": [
    "### 1.1. The Ising Hamiltonian and Observables"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "b59d0179-8ca9-418a-a0e6-81f834c82eed",
   "metadata": {},
   "source": [
    "The Ising model is one of the simplest yet very powerful model in statistical mechanics. It can be used to explain many property of matter such as phase transitions. Here, we want to use quantum simulation to study the property of the Ising model in one dimension. To make things more interesting, we will look at the Ising model in an arbitrary background field."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aed5696b-5232-4447-9305-72be8f05bb4f",
   "metadata": {},
   "source": [
    "<div>\n",
    "    <img src=\"ising_model.png\" width=600>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "536907b6-14d8-4a02-b279-e934d61290d1",
   "metadata": {},
   "source": [
    "In one dimension, the spins in the Ising model is placed along the horizontal X direction. Their spins are pointing in the direction of either $Z$ and $-Z$ direction. The Hamiltonian can be written as follows\n",
    "\n",
    "$$\n",
    "H = -J \\sum_{i=0}^{L-1} \\sigma_i \\sigma_{i+1} - h \\sum_{i=0}^{L-1} \\sigma_i\\cdot \\hat{n}\n",
    "$$\n",
    "\n",
    "where the $J$ is the spin coupling strength and $h$ is magnetic field strength. The unit vector $\\hat{n}$ points in an arbitrary direction in the X-Z plane. For simplicity, we use uniform strength for all spins. The spins $\\sigma_i$ takes values of +1 and -1, and they satisfy the periodic boundary condition, $\\sigma_L = \\sigma_0$."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "fa8db331-4d0f-4642-b010-fa6cf9caa6e2",
   "metadata": {},
   "source": [
    "On the quantum circuit, one can represent the spin up (down) with the quantum state $|0\\rangle$ ($|1\\rangle$). Then, the Ising Hamiltonian can be naturally mapped to Pauli bases by using the mapping of $\\sigma_i \\mapsto Z_i$. The \"qubitized\" Ising Hamiltonian becomes,\n",
    "\n",
    "$$\n",
    "H = -J \\sum_{i=0}^{L-1} Z_i Z_{i+1} - h \\sum_{i=0}^{L-1} \\big(\\sin(\\alpha)Z_i+\\cos(\\alpha)X_i\\big)\n",
    "$$\n",
    "\n",
    "where angle $\\alpha$ is used to parameterize the field. $\\alpha=0$ indicates a purely transverse field. $\\alpha=\\pi/2$ gives a purely vertical field."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f08d668f-7a6c-454c-af70-c5cdc52551b1",
   "metadata": {
    "tags": []
   },
   "source": [
    "<div class=\"alert alert-block alert-success\">\n",
    "\n",
    "**Exercise 7. Building Ising Hamiltonian with an External Field**  \n",
    "\n",
    "Let us first write a function to generate the Ising Hamiltonian with an arbitrary external field. The function you write should work for any lattice of $L$ spins with their spin coupling $J$ and an external field with its strength $h$ and angle $\\alpha$.\n",
    "\n",
    "Tip: since the Ising Hamiltonian is generic, you may want to look at Exercise 2 and use For-loops to generate it. To make it simple, we already wrote part of the code in the function.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "78c3899b-4525-415e-a770-305ba191f534",
   "metadata": {},
   "outputs": [],
   "source": [
    "def H_Ising_operator(L, J, h, alpha=0):\n",
    "\n",
    "    # List of Hamiltonian terms as 3-tuples containing\n",
    "    # (1) the Pauli string,\n",
    "    # (2) the qubit indices corresponding to the Pauli string,\n",
    "    # (3) the coefficient.\n",
    "    ##### ==================================\n",
    "    # Write your solution in here.\n",
    "\n",
    "\n",
    "    \n",
    "    ##### ==================================\n",
    "    \n",
    "    # We create the Hamiltonian as a SparsePauliOp, via the method\n",
    "    # `from_sparse_list`, and multiply by the interaction term.\n",
    "    hamiltonian = SparsePauliOp.from_sparse_list([*ZZ_tuples, *Z_tuples, *X_tuples], num_qubits=L)\n",
    "    return hamiltonian.simplify()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "85db61c6-d66a-44ed-a5f9-a7c1bd4fe856",
   "metadata": {},
   "source": [
    "Test to see if your Hamiltonian make sense at different angles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "73b02511-de96-4818-b4e5-e67e4a87cde4",
   "metadata": {},
   "outputs": [],
   "source": [
    "L = 4\n",
    "J = 0.2\n",
    "h = 1.2\n",
    "alpha = 0\n",
    "\n",
    "H_operator = H_Ising_operator(L=L, J=J, h=h, alpha=alpha)\n",
    "H_operator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6ae8be74-cfd1-4b3e-8f56-e995c4e0524e",
   "metadata": {},
   "outputs": [],
   "source": [
    "L = 4\n",
    "J = 0.2\n",
    "h = 1.2\n",
    "alpha = np.pi/2\n",
    "\n",
    "H_operator = H_Ising_operator(L=L, J=J, h=h, alpha=alpha)\n",
    "H_operator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eac8bebd-542b-438a-bc57-7daa5d3673d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "L = 4\n",
    "J = 0.2\n",
    "h = 1.2\n",
    "alpha = np.pi/4\n",
    "\n",
    "H_operator = H_Ising_operator(L=L, J=J, h=h, alpha=alpha)\n",
    "H_operator"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5a313566-3685-4f6c-a0f7-05b0c4118de7",
   "metadata": {},
   "source": [
    "Now that your Hamiltonian is working. Let us look at how to several useful observables for studying the Ising model.\n",
    "\n",
    "The magnetization operator $M$ is a sum of all the spins in the system, characterizing the total magnetization,\n",
    "$$M = \\sum_{i}^{L-1} \\sigma_i.$$\n",
    "The correlation operator $C$ is a sum of all neighboring spin correlations,\n",
    "$$C = \\frac{1}{L}\\sum_{i}^{L-1} \\sigma_i\\sigma_{i+1}.$$\n",
    "Note these are global observables. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "00e1d407-244d-460d-a437-9b1607ba25f1",
   "metadata": {
    "tags": []
   },
   "source": [
    "<div class=\"alert alert-block alert-success\">\n",
    "\n",
    "**Exercise 8. Building Magnetization and Correlation Operator**  \n",
    "\n",
    "Let us build these global operators $M$ and $C$ that will be used for the simulation. Mapping to Pauli matrices is straightforward and we leave it for you to decide.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f0fdbd63-790a-4c3d-8811-a76d204f7058",
   "metadata": {},
   "outputs": [],
   "source": [
    "def magnetization_op(L):\n",
    "    ##### ==================================\n",
    "    # Write your solution in here.\n",
    "\n",
    "    \n",
    "    ##### ==================================\n",
    "    return magnetization_op\n",
    "\n",
    "def correlation_op(L):\n",
    "    ##### ==================================\n",
    "    # Write your solution in here.\n",
    "\n",
    "    \n",
    "    ##### ==================================\n",
    "    return correlation_op"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "35618ea7-b2e2-4096-8d2e-21efb8820910",
   "metadata": {},
   "outputs": [],
   "source": [
    "L = 6\n",
    "print(magnetization_op(L))\n",
    "print(correlation_op(L))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ca30c3b0-91cc-45ad-b00e-16a314cd4e95",
   "metadata": {},
   "source": [
    "### 1.2. Quantum Circuit Time Evolution (full steps)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "efd6f661-f3fc-4a48-b650-5455601295b7",
   "metadata": {},
   "source": [
    "Now, we have all the ingrediences for the Ising model except for the evolution part. Notice that previously we wrote a *time_evolve* function that evolves a given Hamiltonian $H$ for time duration $t$ based on some trotterization strategy. From the Trotter Formula, we know that in general this will indice an error ~ $\\mathcal{O}(t^2)$ if $H$ contains non-commutating parts. When $t$ is small, this is okay. However, in quantum information science, we are typically interested in the long time simulation in a real-time evolution. \n",
    "\n",
    "Therefore, naturally, to reduce the error, we would want to divide the full time evolution $t$ into many tiny steps $\\Delta t$, and use trotterization scheme within each time slice. Namely, with $N_s = t / \\Delta t$ representing the total number of steps, we write\n",
    "$$ e^{iHt} = \\prod_{i=0}^{N_s}e^{iH\\Delta t} = \\underbrace{e^{iH\\Delta t}e^{iH\\Delta t}... e^{iH\\Delta t}}_{N_s \\text{ many time steps}}. $$\n",
    "In each little $e^{iH\\Delta t}$ evolution, we apply the trotter formula on the internal parts of the $H$.\n",
    "\n",
    "The typical trotter formulas are:\n",
    "* Lie trotter formula (1st order): $e^{i(A+B)\\Delta t} = e^{iA\\Delta t}e^{iB\\Delta t} + \\mathcal{O}((\\Delta t)^2)$   \n",
    "* Suzuki trotter formula (2nd order): $e^{i(A+B)\\Delta t} = e^{iB/2\\Delta t}e^{iA\\Delta t}e^{iB/2\\Delta t} + \\mathcal{O}((\\Delta t)^3)$   \n",
    "* ...\n",
    "\n",
    "We can see that higher order trotterization will reduce the total evolution error. Lie trotterization on each step will give a total trotterization error up to $\\mathcal{O}((\\Delta t)^2)$ and Suzuki trotterization with higher order $p$ gives an error up to $\\mathcal{O}((\\Delta t)^{p+1})$.\n",
    "\n",
    "Nonetheless, higher order trotterizations will introduce many more quantum gates than lower order ones. Having smaller step size will probably be fine with first order trotterization. In principle, there is a trade off between step size $\\Delta t$ and trotterization order $p$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bf4cbd13-b2de-4e08-9640-108d4b7a3b46",
   "metadata": {},
   "source": [
    "Your first task to expand the previously written *time_evolve* function for an arbitrary time period and an arbitrary number of total steps."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "30502ad8-3127-4fdf-85bb-2a523ac225c6",
   "metadata": {
    "tags": []
   },
   "source": [
    "<div class=\"alert alert-block alert-success\">\n",
    "\n",
    "**Exercise 9. Writing a Function to Create a Time Evolution Circuit with Many Trotterization Steps** \n",
    "\n",
    "Your function should takes in an Hamiltonian operator, *H_operator*, and an initial quantum circuit, *initial_circuit*. It will return a quantum circuit that evolves the Hamiltonian for a total number of steps, *total_steps*, where each step size is *dt*.\n",
    "\n",
    "Trotterization strategy should also be passed in as an argument using the *synthesis* variable. By default, we use first-order Lie Trotter formula.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0a080176-88aa-4359-904a-14d2b07ce9a6",
   "metadata": {},
   "outputs": [],
   "source": [
    "def time_evolve_many_steps(H_operator, dt, total_steps, initial_circuit, synthesis=LieTrotter()):\n",
    "    N = H_operator.num_qubits\n",
    "    qc = initial_circuit.copy()\n",
    "    ##### ==================================\n",
    "    # Write your solution in here.\n",
    "\n",
    "\n",
    "    \n",
    "    ##### ==================================\n",
    "    qc.measure_all()\n",
    "    return qc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3528bf8f-564c-4064-9662-d68df03ade14",
   "metadata": {},
   "outputs": [],
   "source": [
    "# test to see if it is working\n",
    "L = 4\n",
    "J = 0.2\n",
    "h = 1.2\n",
    "alpha = np.pi/4\n",
    "\n",
    "H_operator = H_Ising_operator(L=L, J=J, h=h, alpha=alpha)\n",
    "initial_circuit = QuantumCircuit(L)\n",
    "initial_circuit.x(1)\n",
    "initial_circuit.x(2)\n",
    "\n",
    "dt = 0.1\n",
    "total_steps = 2\n",
    "synthesis = LieTrotter()\n",
    "\n",
    "qc = time_evolve_many_steps(H_operator, dt, total_steps, initial_circuit, synthesis)\n",
    "qc.decompose().draw(output=\"mpl\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "947d69d9-402d-4ef6-9673-1843c482f881",
   "metadata": {},
   "source": [
    "### 1.3. Real Time Evolution of the Ising Model - Global Observable"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9fe9508b-fb98-49f9-825f-c34f933874f5",
   "metadata": {
    "tags": []
   },
   "source": [
    "<div class=\"alert alert-block alert-success\">\n",
    "\n",
    "**Exercise 10. Real-time Evolution of the Ising Model**  \n",
    "\n",
    "Now, let us evolve the Ising Hamiltonian in real time. You can use our following code block. Here, we evolve the system and measure the energy, total magnetization, and total correlation, to $t=12$ seconds.\n",
    "\n",
    "Specifically, we start with $L = 6$, and a coupling $J =0.2 \\ll h =1.2$ and $\\alpha=0$. What does it mean? Does your result agree with your expectation?\n",
    "\n",
    "What if we use a smaller time step?\n",
    "\n",
    "You can also play with the parameters to see different effects on the system.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "47abb79f-ec04-4bc7-a53d-dd112eb7c016",
   "metadata": {},
   "outputs": [],
   "source": [
    "###### Ising Model Setup ######\n",
    "L = 6\n",
    "J = 0.2\n",
    "h = 1.2\n",
    "alpha = 0\n",
    "\n",
    "H_operator = H_Ising_operator(L=L, J=J, h=h, alpha=alpha)\n",
    "magnet_op = magnetization_op(L)\n",
    "corr_op = correlation_op(L)\n",
    "\n",
    "initial_spin_downs = [0, L-1]    # spin downs\n",
    "init_qc = QuantumCircuit(L)\n",
    "init_qc.x(initial_spin_downs)    # create initial spin downs at specified locations\n",
    "\n",
    "###### Real-time Evolution Setup ######\n",
    "dt = 0.5\n",
    "total_steps = 24   # t/dt = 12/0.5 = 24\n",
    "total_time = total_steps*dt\n",
    "shots = None\n",
    "synthesis = LieTrotter()\n",
    "\n",
    "###### Quantum Simulation ######\n",
    "lie_exp_vals = []; ts_vals = []\n",
    "t0 = current_time()\n",
    "for steps in range(0, total_steps+1):\n",
    "    time = steps*dt\n",
    "    qc = time_evolve_many_steps(H_operator, dt, steps, initial_circuit=init_qc, synthesis=synthesis)\n",
    "    E_val = estimate_with_shots(qc, H_operator, shots)\n",
    "    M_val = estimate_with_shots(qc, magnet_op, shots)\n",
    "    C_val = estimate_with_shots(qc, corr_op, shots)\n",
    "    print(f'steps = {steps:>4} | time = {time:>5.2f}s | E = {E_val:>5.2f} | M = {M_val:>5.2f} | C = {C_val:>5.2f}')\n",
    "    lie_exp_vals.append([E_val, M_val, C_val])\n",
    "    ts_vals.append(time)\n",
    "print(f\"time elapsed {current_time()-t0:.2f}s\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c74e4ef2-004e-40c0-a988-f0732ff70fa8",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(3, figsize=(8, 6)); fs = 12\n",
    "# simulation\n",
    "lie_exp_vals = np.array(lie_exp_vals)\n",
    "axes[0].plot(ts_vals, lie_exp_vals[:,0], label=\"Lie (1st order)\", marker=\"x\", c=\"darkmagenta\", ls=\"dashed\", lw=0.8)\n",
    "axes[1].plot(ts_vals, lie_exp_vals[:,1], label=\"Lie (1st order)\", marker=\"x\", c=\"darkmagenta\", ls=\"dashed\", lw=0.8)\n",
    "axes[2].plot(ts_vals, lie_exp_vals[:,2], label=\"Lie (1st order)\", marker=\"x\", c=\"darkmagenta\", ls=\"dashed\", lw=0.8)\n",
    "axes[0].set_ylabel(r\"$\\langle H \\rangle$\", fontsize=fs)\n",
    "axes[1].set_ylabel(r\"$\\langle M \\rangle$\", fontsize=fs)\n",
    "axes[2].set_ylabel(r\"$\\langle C \\rangle$\", fontsize=fs)\n",
    "axes[2].set_xlabel(\"Time (s)\", fontsize=fs)\n",
    "fig.suptitle(\"Observable evolution\")\n",
    "legend = fig.legend(*axes[0].get_legend_handles_labels(), bbox_to_anchor=(1.0, 0.5), loc=\"center left\", framealpha=0.5,)\n",
    "fig.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fdaf28a8-2d82-41d5-a6af-bd75d036c9bc",
   "metadata": {},
   "source": [
    "### 1.4. Classical Time Evolution using Matrix Exponential\n",
    "\n",
    "To verify whether our \"quantum\" simulation is correct or not, we directly perform exact matrix exponential on an initial statevector. Here we use *scipy.linalg.expm* to evaluate the time evolution operator $e^{iH\\Delta t}$ approximately at each step. Such approximation is very good for low number of qubits ~ 10 qubits, so they are basically exact. For higher number of qubits, one would need to obtain the eigenvalue and eigenvectors of the Hamiltonian that can be very costy very quickly (Not a subject in our lab section today).\n",
    "\n",
    "Here, we simply gave the code below for you to test. It may take sometime to simulate ~ 40 seconds. In particular, we used double the amount of total steps to make the curve smooth.\n",
    "\n",
    "Does your exact result agree with first order simulation?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9fefb167-a40a-429d-9fc0-69c854d27d36",
   "metadata": {},
   "outputs": [],
   "source": [
    "###### Ising Model Setup ######\n",
    "L = 6\n",
    "J = 0.2\n",
    "h = 1.2\n",
    "alpha = 0\n",
    "\n",
    "H_operator = H_Ising_operator(L=L, J=J, h=h, alpha=alpha)\n",
    "magnet_op = magnetization_op(L)\n",
    "corr_op = correlation_op(L)\n",
    "\n",
    "H_matrix = H_operator.to_matrix()\n",
    "sv_init = Statevector.from_label(\"100001\")\n",
    "\n",
    "###### Real-time Evolution Setup ######\n",
    "dt = 0.5\n",
    "total_steps = 24   # t/dt = 12/0.5 = 24\n",
    "total_time = total_steps*dt\n",
    "exact_ts_val = np.linspace(0, total_time, total_steps*2)\n",
    "\n",
    "###### Exact Classical Simulation ######\n",
    "t0 = current_time()\n",
    "exact_evolution = [sv_init.evolve(LA.expm(-1j * time * H_matrix)) for time in exact_ts_val]\n",
    "print(f\"time elapsed {current_time()-t0:.2f}s\")\n",
    "\n",
    "exact_energy = np.real([sv.expectation_value(H_operator) for sv in exact_evolution])\n",
    "exact_magnetization = np.real([sv.expectation_value(magnet_op) for sv in exact_evolution])\n",
    "exact_correlation = np.real([sv.expectation_value(corr_op) for sv in exact_evolution])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1de1e292-35a7-47a8-92a0-084390d45099",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(3, figsize=(8, 6)); fs = 12\n",
    "# simulation\n",
    "lie_exp_vals = np.array(lie_exp_vals)\n",
    "axes[0].plot(ts_vals, lie_exp_vals[:,0], label=\"Lie (1st order)\", marker=\"x\", c=\"darkmagenta\", ls=\"dashed\", lw=0.8)\n",
    "axes[1].plot(ts_vals, lie_exp_vals[:,1], label=\"Lie (1st order)\", marker=\"x\", c=\"darkmagenta\", ls=\"dashed\", lw=0.8)\n",
    "axes[2].plot(ts_vals, lie_exp_vals[:,2], label=\"Lie (1st order)\", marker=\"x\", c=\"darkmagenta\", ls=\"dashed\", lw=0.8)\n",
    "# exact\n",
    "axes[0].plot(exact_ts_val, exact_energy, label=\"Exact\", c=\"k\", ls=\"-\", lw=1)\n",
    "axes[1].plot(exact_ts_val, exact_magnetization, label=\"Exact\", c=\"k\", ls=\"-\", lw=1)\n",
    "axes[2].plot(exact_ts_val, exact_correlation, label=\"Exact\", c=\"k\", ls=\"-\", lw=1)\n",
    "axes[0].set_ylabel(r\"$\\langle H \\rangle$\", fontsize=fs)\n",
    "axes[1].set_ylabel(r\"$\\langle M \\rangle$\", fontsize=fs)\n",
    "axes[2].set_ylabel(r\"$\\langle C \\rangle$\", fontsize=fs)\n",
    "axes[2].set_xlabel(\"Time (s)\", fontsize=fs)\n",
    "fig.suptitle(\"Observable evolution\")\n",
    "legend = fig.legend(*axes[0].get_legend_handles_labels(), bbox_to_anchor=(1.0, 0.5), loc=\"center left\", framealpha=0.5,)\n",
    "fig.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d9805cba-3613-414b-af2c-7c93d67dd0b6",
   "metadata": {
    "tags": []
   },
   "source": [
    "<div class=\"alert alert-block alert-success\">\n",
    "\n",
    "**Exercise 11. Real-time Evolution of the Ising Model, part 2**  \n",
    "\n",
    "Does your real time simulation agree with the exact simulation? What is going on if they are in agreement? How should you fix them? \n",
    "\n",
    "Let us see if you can produce a real time simulation that matches the exact solution up to 20 seconds.\n",
    "\n",
    "Hint: trotterization step, or higher order trotter formula may be helpful?\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b8012200-77e7-4eef-8bf6-b3b14267cbf4",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "44d01960-ed0f-4894-8845-4ea31097ef8f",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fe7cc5d7-d7f4-40fa-8527-537257750c77",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "67366a40-26ad-4abc-b3b0-d676ad9a1401",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "bffad13f-fa87-4108-8ce6-fef66a51a212",
   "metadata": {},
   "source": [
    "### 1.5. Real Time Evolution of the Ising Model - Local Observable"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "28385062-42eb-49e4-b309-142b70714304",
   "metadata": {},
   "source": [
    "It may also be fun to visualize the magnetization at each site in the lattice and observe how they evolve over time. These local magnetization operator at each spin site $i$ is basically $M_i = \\sigma_i$. With exact simulation result, we can see the magnetization damps as time evolves."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "01e1970e-2314-4da1-89f4-66392504eb2d",
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import cm\n",
    "\n",
    "# An inner list comprehension loops over the terms of the SparsePauliOp magnetization_op, which corresponds to the magnetization of each one of the sites\n",
    "local_magnet = np.real([[sv.expectation_value(term) for term in magnet_op] for sv in exact_evolution])\n",
    "\n",
    "# The shape of magnetizations is (101, 6), containing <Z>(t) for each site 0, 1, ..., 5\n",
    "plt.figure(figsize=(14, 2))\n",
    "\n",
    "# Create the 2-dim xx and yy arrays tiling the grid with the x and y values\n",
    "xx, yy = np.meshgrid(exact_ts_val, np.arange(L))\n",
    "plt.pcolor(xx, yy, local_magnet.T, vmin=-1, vmax=+1, cmap=\"RdBu\")\n",
    "\n",
    "# Force the figure to have all y ticks from 0 to 5\n",
    "plt.yticks(np.arange(L))\n",
    "plt.ylabel(\"Site $i$\")\n",
    "plt.xlabel(\"Time\")\n",
    "plt.colorbar(label=\"$\\\\langle Z_i \\\\rangle$\", aspect=1.8)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "973c5573-c1c2-409a-ba7d-e170c56a8311",
   "metadata": {
    "tags": []
   },
   "source": [
    "<div class=\"alert alert-block alert-success\">\n",
    "\n",
    "**Exercise 12. Real-time Evolution of the Ising Model, part 3**  \n",
    "\n",
    "Let us use your best simulation and compare to see with the classical exact results.\n",
    "\n",
    "Hint: you may need to re-run to compute *local_magnet* results and also adjust the *exact_ts_val* to the *ts_val* used in your quantum simulation code.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1b578688-7fee-4af4-a978-4b5df6f0ec54",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "85e73eeb-45dc-48c0-af31-6d4ef2a9b7c3",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4f68492f-89eb-4401-8ae0-548157c8f0fd",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "eac9bb97-cadf-4235-bc41-88c33c4b8288",
   "metadata": {},
   "source": [
    "### Open Questions"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "68bc3a34-52d5-4e56-85ca-47dcab8cd285",
   "metadata": {
    "tags": []
   },
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "\n",
    "**Further Questions**  \n",
    "* Faster simulation with optimal trotter step size and trotter order, for example using *TrotterQRTE* class from *quantum_algorithm* module. See tutorial [here](https://qiskit-community.github.io/qiskit-algorithms/tutorials/13_trotterQRTE.html). Important note: on real quantum device, due to no-cloning theorem, one needs to restart the full evolution for any observable.\n",
    "* Comparing with Classical Monte Carlo simulation using Metropolis Algorithm to observe phase transition at high and low temperature limits.\n",
    "* Generalization to 2D Ising Model\n",
    "* Generalization to finite temperature Ising Model evolution\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3c99d775-a6d2-4a7f-bf6a-622ae4b34919",
   "metadata": {},
   "source": [
    "## Part 2. Maximum Cut Problem"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fbdd5ea7-72e9-4687-90fc-c5de6c49c1ca",
   "metadata": {},
   "source": [
    "The well-known Maximum cut problem (or Maxcut problem) is stated as follows:\n",
    "\n",
    "\"Given a graph of n nodes, how do we cut the graph into two partitions such that the number of edges between the two partitions are maximized?\"\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4268937e-2889-42a7-a4d9-0a9fe5d488f7",
   "metadata": {},
   "source": [
    "<div>\n",
    "    <img src=\"maxcut_problem.png\" width=300>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "34c1ce50-52c8-4a2b-9b6a-518593c9c486",
   "metadata": {},
   "source": [
    "This so-called Maxcut problem, is an NP-hard problem in general. It has excited many algorithms to calcuate its solution. In general, a bruteforce that enumerate all possible partitions will break down at large graph limit. People have dived into approximation algorithms to solve this combinatorial problem. It has become an important subfield in Discrete Mathematics.\n",
    "\n",
    "With the advent of quantum computing, the Maxcut problem has been applied extensively using variational quantum algorithms and for some specific subfield of this interesting problem, the quantum approach is able to defeat the classical approach theoretically. \n",
    "\n",
    "Today, we will look at the formulation of a simplified Maxcut problem (i.e., no constraints, no edge weight, etc), and we will use one of the popular quantum simulation techniques, the variational quantum eigensolver, to solve it."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fd6f2977-c8be-4b20-9e1c-c4b35467a475",
   "metadata": {},
   "source": [
    "### 2.1. Play with the Maximum Cut Problem"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "499155a3-ebd0-4f26-92df-a2a86f4f318d",
   "metadata": {},
   "source": [
    "We start by loading the networkx module and some helper functions..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "605a1d0f-03b6-43ea-9c09-d5cafc46a7a7",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "\n",
    "# load helpful functions\n",
    "def prepare_maxcut(n, edges):\n",
    "    '''\n",
    "    n: number of nodes in G\n",
    "    edges: list of connected edge pairs in G, make sure edges are indexed from 0, 1, to (n-1)\n",
    "    '''\n",
    "    for ei, ej in edges:\n",
    "        assert 0 <= ei < n and 0 <= ej < n, f\"Invalid edges! Edge ({ei}, {ej}) is Not allowed in a {n}-node Graph.\"\n",
    "    G = nx.Graph()\n",
    "    G.add_nodes_from(np.arange(0, n, 1))\n",
    "    G.add_edges_from(edges)\n",
    "    pos = nx.shell_layout(G)  # fix plot position, more options see \n",
    "    return G, pos\n",
    "\n",
    "def plot_maxcut(G, pos, colors=None):\n",
    "    '''\n",
    "    Plot the Maxcut problem instance\n",
    "    '''\n",
    "    if colors == None:\n",
    "        colors = [\"lightgray\" for node in G.nodes()]\n",
    "    default_axes = plt.axes(frameon=True)\n",
    "    nx.draw_networkx(G, node_color=colors, node_size=1000, alpha=1.0, ax=default_axes, pos=pos)\n",
    "    nx.draw_networkx_edges(G, pos=pos)\n",
    "    \n",
    "def get_adj_matrix(G):\n",
    "    '''\n",
    "    Return an adjacency matrix of G, representing w_ij\n",
    "    '''\n",
    "    w = np.zeros([n, n])\n",
    "    for i in range(n):\n",
    "        for j in range(n):\n",
    "            temp = G.get_edge_data(i, j, default=0)\n",
    "            if temp != 0:\n",
    "                w[i, j] = 1\n",
    "    return w\n",
    "\n",
    "def get_bruteforce_sol(G, verbose=False):\n",
    "    '''\n",
    "    Return the bruteforce solution for the adjacency matrix of a Maxcut graph.\n",
    "    '''\n",
    "    G_adj = get_adj_matrix(G)\n",
    "    best_cost_brute = 0\n",
    "    for b in range(2**n):\n",
    "        # b = 0, x = [0, 0, 0, 0]\n",
    "        # b = 1, x = [1, 0, 0, 0]\n",
    "        # b = 2, x = [0, 1, 0, 0] ...\n",
    "        x = [int(t) for t in reversed(list(bin(b)[2:].zfill(n)))]\n",
    "        cost = 0\n",
    "        for i in range(n):\n",
    "            for j in range(n):\n",
    "                cost = cost + G_adj[i, j] * x[i] * (1 - x[j])\n",
    "        if best_cost_brute < cost:\n",
    "            best_cost_brute = cost\n",
    "            xbest_brute = x\n",
    "        if verbose: print(\"case = \" + str(x) + \" cost = \" + str(cost))\n",
    "    return xbest_brute, best_cost_brute"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6d73f9e6-4511-4af1-8697-eb5ce5d54075",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">\n",
    "\n",
    "**Exercise 13. Get Familiar with Maxcut Problem**  \n",
    "    \n",
    "Play with these functions above \n",
    "- *prepare_maxcut*\n",
    "- *plot_maxcut*\n",
    "- *get_adj_matrix*\n",
    "- *get_bruteforce_sol*\n",
    "    \n",
    "See if you can come up with an interesting Maxcut problem. \n",
    "    \n",
    "Warning, be careful with the number of nodes, n, in the Graph. Your computer might crash if n is too big.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "df7061e7-359c-48e4-b6f1-eae3b8dcaf8b",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Example 1.\n",
    "n = 4\n",
    "edges = [[0, 1], [1, 2], [2, 3], [3, 0]]\n",
    "G, pos = prepare_maxcut(n, edges)\n",
    "print(\"Adjacency matrix W = \\n\", get_adj_matrix(G))\n",
    "plot_maxcut(G, pos)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f7356f6e-b5e4-4e30-8d7d-eafcc2d70405",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "xbest_brute, best_cost_brute = get_bruteforce_sol(G, verbose=True)\n",
    "print(\"\\nBest solution = \" + str(xbest_brute) + \"  Maxcut cost = \" + str(best_cost_brute))\n",
    "colors = [\"salmon\" if xbest_brute[i] == 0 else \"cornflowerblue\" for i in range(n)]\n",
    "plot_maxcut(G, pos, colors)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b488bd9f-ec41-4844-8826-4c2c6f486cee",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Example 2.\n",
    "n = 6\n",
    "edges = [[0, 1], [1, 2], [2, 3], [3, 4], [4, 5], [5, 0], [3, 0], [1, 3], [3, 5]]\n",
    "G, pos = prepare_maxcut(n, edges)\n",
    "print(\"Adjacency matrix W = \\n\", get_adj_matrix(G))\n",
    "plot_maxcut(G, pos)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6ef88f9c-2e68-498a-9794-f29f84dc742f",
   "metadata": {},
   "outputs": [],
   "source": [
    "xbest_brute, best_cost_brute = get_bruteforce_sol(G, verbose=True)\n",
    "print(\"\\nBest solution = \" + str(xbest_brute) + \"  Maxcut cost = \" + str(best_cost_brute))\n",
    "colors = [\"salmon\" if xbest_brute[i] == 0 else \"cornflowerblue\" for i in range(n)]\n",
    "plot_maxcut(G, pos, colors)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4b5b60f8-8ac9-4161-a759-9ed7268e6627",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Example 3.\n",
    "n = 8\n",
    "edges = [[i, i+1] for i in range(0, n-1)] + [[n-1, 0]]\n",
    "edges += [[0, 2], [1, 3], [2, 4], [3, 7], [2, 6], [3,4], [6,7], [5,7], [0, 4], [0, 6]]\n",
    "G, pos = prepare_maxcut(n, edges)\n",
    "xbest_brute, best_cost_brute = get_bruteforce_sol(G)\n",
    "print(\"Best solution = \" + str(xbest_brute) + \"  Maxcut cost = \" + str(best_cost_brute))\n",
    "colors = [\"salmon\" if xbest_brute[i] == 0 else \"cornflowerblue\" for i in range(n)]\n",
    "plot_maxcut(G, pos, colors)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6bc88246-92f7-400f-86b0-13ea172d4476",
   "metadata": {},
   "source": [
    "Did you come up with something more interesting?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "96cc08b2-9652-4e37-87d1-e95380bfd3f7",
   "metadata": {},
   "source": [
    "### 2.2. The Maxcut Hamiltonian"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "607ea762-3182-4324-bf09-7f7cf568b91d",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">\n",
    "\n",
    "**Exercise 14. Building the Maxcut Hamiltonian**  \n",
    "    \n",
    "From the lecture, we know the Maxcut Hamiltonian looks like,\n",
    "    \n",
    "$$    H = \\sum_{ij} w_{ij}x_i(1-x_j) = \\sum_{ij} w_{ij}\\frac{I_i-Z_i}{2}\\Big(I- \\frac{I_j-Z_j}{2}\\Big) = \\sum_{i} c_{i} Z_i + \\sum_{ij} d_{ij} Z_iZ_j$$\n",
    "    \n",
    "Use the helper function below (or come up with your own) to build the Hamiltonian for Maxcut.\n",
    "    \n",
    "Important to verify, your final, simplied Hamiltonian must only contain Pauli Z expressions.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2a0fd62e-d80f-4772-b665-a841dea5fc63",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# helpful functions\n",
    "def get_x_op_at_i(i, n):\n",
    "    ''' Return an operator x_i = (I_i - Z_i)/2 at specified qubit index i, with n qubits in total '''\n",
    "    candidates = [I]*n\n",
    "    candidates[i] = (I - Z)/2\n",
    "    res = candidates[0]\n",
    "    for i in range(1, n):\n",
    "        res = res^candidates[i]\n",
    "    return res\n",
    "\n",
    "def get_I_op(n):\n",
    "    ''' Return an identity operator over n qubits'''\n",
    "    return SparsePauliOp([\"I\"*n], [1])\n",
    "\n",
    "print(get_x_op_at_i(1, 4))\n",
    "print(get_I_op(4))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "64e24f21-c5f3-49e8-b624-e7d16376be0d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def get_Hamiltonian_of_graph(G):\n",
    "    '''\n",
    "    Get Pauli-ZZ Hamiltonian for the adjacency matrix of a graph\n",
    "    '''\n",
    "    G_adj = get_adj_matrix(G)\n",
    "    n = len(G_adj)\n",
    "    Hamiltonian = 0\n",
    "\n",
    "    for i in range(n):\n",
    "        for j in range(n):\n",
    "            if G_adj[i][j] != 0: # to speed up, avoid 0\n",
    "                ##### ==================================\n",
    "                # Write your solution in here.  \n",
    "                # basically, you try to add all nonzero pauli terms to the Hamiltonian variable\n",
    "                # each pauli term is a product of two pauli operator\n",
    "                \n",
    "                \n",
    "                ##### ==================================\n",
    "                \n",
    "    Hamiltonian = Hamiltonian.simplify()  # simplify algebraic expression\n",
    "    return Hamiltonian"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9397fb20-2eda-4a6b-9a0f-c5d79b50fc98",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Verification, Hamiltonian must contain only Pauli-Z\n",
    "# Try start with a simple graph first, then move to a larger graph.\n",
    "n = 4\n",
    "edges = [[0, 1], [1, 2], [2, 3], [3, 0]]\n",
    "G, pos = prepare_maxcut(n, edges)\n",
    "H = get_Hamiltonian_of_graph(G)\n",
    "print(H)\n",
    "plot_maxcut(G, pos)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "50d34cae-0cac-47f1-aabe-095c9c04ff14",
   "metadata": {},
   "source": [
    "### 2.3. The Variational Quantum Eigensolver - Ansatz"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5df42207-d3a7-462a-9b21-eeaf4c4efa47",
   "metadata": {},
   "source": [
    "Let us first quickly review what is the Variational Quantum Eigensolver (VQE). The VQE optimization loop consists of three essential parts\n",
    "\n",
    "- Variational ansatz\n",
    "\n",
    "- Quantum evolution and expectation evaluation, i.e. $\\langle H \\rangle = \\langle \\psi | H | \\psi \\rangle $\n",
    "    \n",
    "- Classical optimizer to update the next parameter, or exit the loop.\n",
    "\n",
    "Instead of just calling some function from Qiskit, it is very important to know how to write the Vartional Quantum Eigensolver (VQE) from scratch. Here, we will help you to build your own VQE algorithm piece by piece.\n",
    "\n",
    "The first task is to write a variational ansatz. This is probably the most diversified part of a VQE algorithm. One could implement this part like before using a time evolution of some Hamiltonian ansatz. Instead, we show you how to implement this using a heuristic variational circuit that is a good representation of the Hilbert space. \n",
    "\n",
    "For this part, we use the so-called Two-Local ansatz (available as *TwoLocal* class from Qiskit), where only adjacent qubit are connected (entangled). It is a commonly-used hardware-efficient ansatz, commonly used in quantum optimization problems."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0e268a2b-be5f-494f-b87c-171245d0eecd",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from qiskit.circuit.library import TwoLocal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "545b2e82-3f96-439c-bfaa-63a538bd39ad",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# An example on how to build 4-qubit TwoLocal ansatz with \"ry\" rotation gates, \"cz\" entanglement gates, and a depth of 5 repeated layers. \n",
    "# which is good starting place for us.\n",
    "ansatz = TwoLocal(4, rotation_blocks=[\"ry\"], entanglement_blocks=[\"cz\"], reps=5, entanglement=\"linear\")\n",
    "ansatz.decompose().draw(output='mpl')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a3aaf63c-990f-4a94-ae4e-0884f1800046",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">\n",
    "\n",
    "**Exercise 15. Vartional Quantum Eigensolver from Scratch, part 1**  \n",
    "    \n",
    "Looking at the example above and let's write a function to build our ansatz for a number of qubits $n$ and an arbitrary number of layers. You can keep the other function arguments same as the previous example.\n",
    "    \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e225431e-db43-4d37-a5b0-5c4ff92e6bb1",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# build our own function using the example above \n",
    "def get_ansatz(n, layers):\n",
    "    ##### ==================================\n",
    "    # Write your solution in here.  \n",
    "\n",
    "    \n",
    "    ##### ==================================\n",
    "    return ansatz\n",
    "\n",
    "# test\n",
    "get_ansatz(6, 3).decompose().draw(output='mpl')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e6aae0b1-c783-46e2-8e15-872695e19e11",
   "metadata": {},
   "source": [
    "### 2.4. The Variational Quantum Eigensolver - Cost Function"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "45f0ed85-4e9c-4448-ae8c-701b2a455644",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">\n",
    "\n",
    "**Exercise 14. Vartional Quantum Eigensolver, part 2**  \n",
    "    \n",
    "We now look at the measuring the observable.\n",
    "    \n",
    "We write an *get_expectation* function to take in a variantional ansatz, a list of parameters, an observable, and a shot number. It returns the expectation value. \n",
    "    \n",
    "Try take advantage of our own function *estimate_with_shots* defined earlier.\n",
    "    \n",
    "Note, the ansatz qubits and the observable qubits must be the identical.\n",
    "        \n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8d6f93bb-033b-45c1-86a9-346abc0ccb52",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def get_expectation(ansatz, params, observable, shots):\n",
    "    assert ansatz.num_qubits == observable.num_qubits, f\"ansatz qubits = {ansatz.num_qubits}, observable qubits = {observable.num_qubits}\"\n",
    "    ##### ==================================\n",
    "    # Write your solution in here, you can use the method .bind_parameters() to inject parameters into the ansatz\n",
    "\n",
    "    \n",
    "    ##### ==================================\n",
    "    return expectation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "904ed3c3-0fd1-435e-9404-f01e5c183dc5",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# test\n",
    "# maxcut problem instance\n",
    "n = 4\n",
    "edges = [[0, 1], [1, 2], [2, 3], [3, 0]]\n",
    "G, pos = prepare_maxcut(n, edges)\n",
    "\n",
    "# building ansatz and observable\n",
    "ansatz_ex = get_ansatz(n, 3)\n",
    "parameters = [0.1]*ansatz_ex.num_parameters\n",
    "observable = get_Hamiltonian_of_graph(G)\n",
    "\n",
    "# expectation\n",
    "get_expectation(ansatz_ex, parameters, observable, 1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "59c933d7-4f5a-45bd-8e54-d0fe7640c13e",
   "metadata": {},
   "source": [
    "Note, your *get_expectation* function should return a value, which is the expectation value of quantum state generated by the ansatz. By tuning the parameters, you will see that this expectation value is different.\n",
    "\n",
    "Essentially, the VQE algorithm finds the optimal parameter of the ansatz such that your expectation is minimizes. By the Variational Principle, when the expectation minimizes, our guess wavefunction approaches the ideal wavefunction.\n",
    "\n",
    "In short, your ansatz with optimal parameters becomes your ideal wavefunction."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c86ea341-69ce-4045-b93c-1a958f57cddf",
   "metadata": {},
   "source": [
    "### 2.5. The Variational Quantum Eigensolver - Optimizer"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "139ee1f2-9850-4c72-865e-f4f4c27bd484",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">\n",
    "\n",
    "**Exercise 16. Vartional Quantum Eigensolver, part 3**  \n",
    "    \n",
    "We now look at the last step, which is the optimization using classical optimizer, to complete our VQE function.\n",
    "    \n",
    "Our *vqe_from_scratch* takes a Maxcut problem instance G, a variational ansatz, and number of shots.\n",
    "    \n",
    "Specifically, we use Qiskit's *COBYLA* optimizers. Since *COBYLA* only has a *minimize* method, and we try to maximize the expectation observable (because it is called Maxcut!), we should evaluate $-\\langle H \\rangle$ in the *COBYLA.minimize* method. Maximum iterations, i.e., *maxiter*, should also be supplied.\n",
    "        \n",
    "To make results somewhat comparable, we add *random_seed* argument for initializing the parameter. However, simulation itself still involves randomness (a topic for another day).\n",
    "\n",
    "There are a lot of optimizers you can use if you are interested, available in [qiskit.algorithms.optimizers](https://qiskit.org/documentation/stubs/qiskit.algorithms.optimizers.html). If you are using other classes, call_back functions need modify, see comments within.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8c9eb22c-309b-44a4-980c-b5c03cfa5efa",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from qiskit_algorithms.optimizers import COBYLA\n",
    "# from qiskit.algorithms.optimizers import COBYLA # if last line does not work, use this one (but I recommend first line)\n",
    "\n",
    "def maxcut_vqe_from_scratch(G, ansatz, shots, maxiter=300, random_seed=0, verbose=True):\n",
    "    '''\n",
    "    Our VQE from scratch for Maxcut problem G\n",
    "    \n",
    "    Return the VQE optimization result and its history\n",
    "    \n",
    "    '''\n",
    "    \n",
    "    # Helper function\n",
    "    exp_val_list = []\n",
    "    def callback(params):\n",
    "        '''\n",
    "        Callback to record expectation history and more!\n",
    "           It is always optimizer dependent, so very tricky in general, for other optimizers, \n",
    "           you will need modify the argument to callback(num_eval, param, f_eval, step, acc) instead\n",
    "        '''\n",
    "        exp_val = get_expectation(ansatz, params, observable, shots)\n",
    "        exp_val_list.append(exp_val)\n",
    "        if verbose == True: print(f\"step {len(exp_val_list):>4}:  <H> = {exp_val:>6.3f}\")\n",
    "    \n",
    "    initial_parameters = [0.1]*ansatz_ex.num_parameters\n",
    "    observable = get_Hamiltonian_of_graph(G)\n",
    "    \n",
    "    def minimize_func(params):\n",
    "        ##### ==================================\n",
    "        # Write your solution in here.     \n",
    "\n",
    "        \n",
    "        ##### ==================================\n",
    "        return value_to_maximize\n",
    "\n",
    "    np.random.seed(random_seed)\n",
    "    initial_params = 0.1* np.random.random(ansatz.num_parameters)\n",
    "    opt = COBYLA(maxiter=maxiter, callback=callback)\n",
    "    \n",
    "    res = opt.minimize(minimize_func, initial_params)\n",
    "    \n",
    "    return res, exp_val_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "58d5a0b5-89eb-44b7-8532-c571e993799b",
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "# test\n",
    "# maxcut problem instance\n",
    "n = 4\n",
    "edges = [[0, 1], [1, 2], [2, 3], [3, 0]]\n",
    "G, pos = prepare_maxcut(n, edges)\n",
    "\n",
    "# building ansatz and observable\n",
    "ansatz = get_ansatz(n, 1)\n",
    "\n",
    "# vqe\n",
    "opt_res, exp_val_list = maxcut_vqe_from_scratch(G, ansatz, shots=1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6605d741-e985-4d8a-90b2-8d7aaa236674",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def plot_vqe_exp_val_list(exp_val_list, G=None):\n",
    "    N = len(exp_val_list)\n",
    "    xs = np.arange(1, N+1)\n",
    "    plt.scatter(xs, exp_val_list, label='quantum simulation')\n",
    "    plt.ylabel(r'Maxcut', fontsize=20)\n",
    "    plt.xlabel('Iteration step', fontsize=20)\n",
    "    plt.xticks(fontsize=15)\n",
    "    plt.yticks(fontsize=15)\n",
    "    if G != None:\n",
    "        xBest, costBest = get_bruteforce_sol(G)\n",
    "        plt.hlines(costBest, 0, N, ls='dashed', color='red', label='bruteforce')\n",
    "    plt.legend(fontsize=15)\n",
    "    plt.show()\n",
    "    \n",
    "# visualize result   \n",
    "plot_vqe_exp_val_list(exp_val_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "21b41717-3f0c-452c-8492-be2497c9bba1",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# visualize with bruteforce\n",
    "plot_vqe_exp_val_list(exp_val_list, G)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e66075a6-763b-48b7-916b-8705704ace11",
   "metadata": {},
   "source": [
    "If you see your simulation result agrees with the bruteforce as step increases, that indicates your algorithm is correct!\n",
    "\n",
    "Big congratulations! You have implemented VQE from scratch, essentially."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b7731712-764a-4a4a-a8de-a46a2d92f46a",
   "metadata": {},
   "source": [
    "### 2.6. Extract Maxcut Solution from Quantum Simulation "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5114bdfb-16ee-4a8c-8f82-e68d56fbde78",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-success\">\n",
    "\n",
    "**Exercise 17. Vartional Quantum Eigensolver, part 4**  \n",
    "    \n",
    "Hopefully your VQE optimization is working and you obtain the expected result with bruteforce calculation. Congratulations! It is a big step.\n",
    "    \n",
    "Now, you may wonder how do I know my optimal solution, i.e., the binary variable $x$? Afterall, we also need to know how to partition the graph, where the maximum cut is reached? \n",
    "    \n",
    "The answer is you just need to sample the quantum state with your final, optimal parameters of the ansatz! Use your *sample_with_shots* or *sample_with_shots_v2* class defined before.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "45fc13f1-e2ea-4ae4-ab34-d3be129f871b",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# your optimal parameters is in the x method of your opt_res object\n",
    "opt_res.x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f5c141ef-0ee5-488a-9b72-eb3674b4bbbc",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def sample_probability_dist(ansatz, optimization_result, shots):\n",
    "    ##### ==================================\n",
    "    # Write your solution in here. \n",
    "    \n",
    "\n",
    "    \n",
    "    ##### ==================================\n",
    "    return dist"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0e3a98e1-14cd-487f-b65b-c19051e9bc9c",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def plot_probability_dist(dist):\n",
    "    '''\n",
    "    Convenience class to plot distribution\n",
    "    '''\n",
    "    if type(list(dist.keys())[0]) != str: # depend on your sample function\n",
    "        new_dist = {}\n",
    "        for k, v in dist.items():\n",
    "            k = \"\".join([t for t in reversed(list(bin(k)[2:].zfill(n)))])\n",
    "            new_dist[k] = v\n",
    "        dist = new_dist\n",
    "    \n",
    "    return plot_histogram(dist)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "560e1847-dce7-4774-b21a-63af279bc365",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "dist = sample_probability_dist(ansatz, opt_res, 1024)\n",
    "print(dist)\n",
    "plot_distribution(dist)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f9b5be6d-3baf-4ca2-a47b-9933fe0c804f",
   "metadata": {},
   "source": [
    "Note both 0101 and 1010 are correct solutions, because it does not matter if you call which one in partition which. It is just a label. So we literally achieved nearly 100 % accuracy!\n",
    "\n",
    "Now, checking with bruteforce solution : )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0c5270a5-3c32-42b4-b01a-1623c5ab9db6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# cross check with bruteforce\n",
    "xbest_brute, best_cost_brute = get_bruteforce_sol(G, verbose=False)\n",
    "print(\"Best solution = \" + str(xbest_brute) + \"  Maxcut cost = \" + str(best_cost_brute))\n",
    "colors = [\"salmon\" if xbest_brute[i] == 0 else \"cornflowerblue\" for i in range(n)]\n",
    "plot_maxcut(G, pos, colors)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3980e1f8-8d30-4e42-b3e9-2ffe0f216d2e",
   "metadata": {
    "jp-MarkdownHeadingCollapsed": true
   },
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "5ecf18a8-f090-4042-839f-6137eea617b1",
   "metadata": {},
   "source": [
    "### 2.7. Experiment with the Maxcut Problem"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "41385d07-111a-4400-ab57-ccd20115c354",
   "metadata": {},
   "source": [
    "Now it is time to test your VQE algorithm on some randomly generated maxcut problems.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8359d2d5-1437-4d70-b8d7-03e96771c3ce",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# You can use my function to generate a randomly hard maxcut graph\n",
    "def prepare_random_maxcut(n, connectivity=0.3):\n",
    "    '''\n",
    "    Some hard maxcut graph generator, not an ideal one\n",
    "    '''\n",
    "    assert 0 <= connectivity < 1, \"Connectivity is between 0 and 1\"\n",
    "    edges = []\n",
    "    for i in range(n):\n",
    "        for j in range(n):\n",
    "            if i == j: continue\n",
    "            e, re = [i, j], [j, i]\n",
    "            if re in edges: continue # can be optimized with cache\n",
    "            if np.random.random() < connectivity:\n",
    "                edges.append(e)\n",
    "    G, pos = prepare_maxcut(n, edges)\n",
    "    return G, pos"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b38d9c31-c480-409c-8aa5-ff3926332bcd",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "n = 10\n",
    "G, pos = prepare_random_maxcut(n, 0.3)\n",
    "xbest_brute, best_cost_brute = get_bruteforce_sol(G)\n",
    "print(\"Best solution = \" + str(xbest_brute) + \"  Maxcut cost = \" + str(best_cost_brute))\n",
    "colors = [\"salmon\" if xbest_brute[i] == 0 else \"cornflowerblue\" for i in range(n)]\n",
    "plot_maxcut(G, pos, colors)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c09e0de6-226e-4f11-8dc1-5617b714f32b",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# building ansatz\n",
    "ansatz = get_ansatz(n, 1)\n",
    "ansatz.decompose().draw(scale=0.8, output='mpl')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a4bdf1ba-21b8-44d1-bdcd-e174c56343aa",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# vqe optimization\n",
    "opt_res, exp_val_list = maxcut_vqe_from_scratch(G, ansatz, shots=3000, maxiter=500, random_seed=42, verbose=False)\n",
    "plot_vqe_exp_val_list(exp_val_list, G)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c97b78fc-b285-40c2-afee-44f491f9a194",
   "metadata": {},
   "source": [
    "### Open Questions"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11be1c9a-8f7c-4fab-bc55-ff205cf941f0",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-warning\">\n",
    "\n",
    "**Further Questions**  \n",
    "    \n",
    "* Improve your VQE function to solve the Maxcut Problem\n",
    "* Play with your *maxcut_vqe_from_scratch* on a more complicated Maxcut problem to see how it behaves! At some point, for some graph, the algorithm becomes insufficient, indicating the room for improvement such as the optimizer, expectation evaluations, etc. \n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ca5f8994-0aa8-41ed-a0f1-b870baa50f6b",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}