Hands-on.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### If you're using your own machine, install qiskit first\n",
    "To do that, uncomment the line below and run cell"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# !pip install qiskit\n",
    "# To update qiskit:\n",
    "# !pip install -U qiskit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### And check version — should match the one below"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'qiskit-terra': '0.15.2',\n",
       " 'qiskit-aer': '0.6.1',\n",
       " 'qiskit-ignis': '0.4.0',\n",
       " 'qiskit-ibmq-provider': '0.8.0',\n",
       " 'qiskit-aqua': '0.7.5',\n",
       " 'qiskit': '0.20.1'}"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import qiskit\n",
    "qiskit.__qiskit_version__"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Let's solve MAXCUT on a simple graph\n",
    "\n",
    "Recall the MAXCUT objective is $\\max_{\\mathbb{s}}\\frac{1}{2}\\sum_{ij\\in E}(1-s_is_j)$\n",
    "\n",
    "Before going further, let's reformulate it:\n",
    "\n",
    "$\\frac{1}{2}\\sum_{ij\\in E}(1-s_is_j) = \\frac{1}{2}\\sum_{ij\\in E}(-s_is_j) + const$\n",
    "\n",
    "Therefore the above _maximization_ problem is equivalent to the problem below:\n",
    "\n",
    "$\\max_{\\mathbb{s}}\\sum_{ij\\in E}(-s_is_j)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Constructing the Hamiltonian\n",
    "\n",
    "MAXCUT is equivalent to the problem below:\n",
    "\n",
    "$\\max_{\\mathbb{s}}\\sum_{ij\\in E}(-s_is_j)$\n",
    "\n",
    "We've seen earlier that the Hamiltonian can be constructed by mapping binary variables $s_i$ onto the eigenvalues of $Z$. For the above maximization problem has the following Hamiltonian:\n",
    "\n",
    "$\\sum_{ij\\in E}(-Z_iZ_j)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Constructing QAOA circuit\n",
    "\n",
    "Cost Hamiltonian: $C = \\sum_{ij\\in E}(-Z_iZ_j)$\n",
    "\n",
    "The full QAOA circuit:\n",
    "\n",
    "$|\\psi{(\\mathbb{\\beta},\\mathbb{\\gamma})}> =  e^{-i\\beta_p B}e^{-i\\gamma_p C}\\cdots e^{-i\\beta_1 B}e^{-i\\gamma_1 C} H^{\\otimes n} |0>$\n",
    "\n",
    "Need to implement cost operator $e^{-i\\gamma_p C}$ and mixer operator $e^{-i\\beta_p B}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Let's beging by importing everything we'll need later"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "# import warnings; warnings.filterwarnings(\"ignore\")\n",
    "import numpy as np\n",
    "import networkx as nx\n",
    "import matplotlib.pyplot as plt\n",
    "from collections import defaultdict\n",
    "from operator import itemgetter\n",
    "from scipy.optimize import minimize\n",
    "from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister, execute, Aer\n",
    "from qiskit.aqua.algorithms import NumPyEigensolver\n",
    "from qiskit.quantum_info import Pauli\n",
    "from qiskit.aqua.operators import op_converter\n",
    "from qiskit.aqua.operators import WeightedPauliOperator"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Consider the following graph"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAb4AAAEuCAYAAADx63eqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAWu0lEQVR4nO3df2yV9WLH8c9pT+FgsBcjnZBBwi69UqbUXNkNLN7E6iIuhJssV2p0kpkZqUtNpmYu20I0KprNhOzejGhu7D86vVlqm+vtlkDSOxVNdGoEpEatBp0L3UDABA+N9Bc9+6MePJSe9pzn+T7P8/3xfv0L58vTv958es55nlypVCoJAIBANGR9AQAApInwAQCCQvgAAEEhfACAoBA+AEBQCB8AICiEDwAQFMIHAAgK4QMABIXwAQCCkk/6Hzg9Oq7+gyMaPlFUcWxKzYW82lY0q3PjKl25dHHS/zwAwGJZNCKX1L06jxw7o2cOHNUbn52SJI1PTV/4s0K+QSVJHeta1H1jq65bvSyJSwAAWCrLRiQSvpfe+VJP7RvW2NR5zXd6LicV8o3atbVNOzavMX0ZAAALZd0I4+Gb+YE+0bnJ6YX/8neWNDVo19b1xA8APGdDI4x+uOXIsTN6at9wXT+QJJ2bnNZT+4Y1NHLG5OUAACxiSyOMhu+ZA0c1NnU+0mvHps7r2QNHTV4OAMAitjTCWPhOj47rjc9Ozfv72vmUStLrn57S16Pjpi4JAGAJmxphLHz9B0din5GT1H8o/jkAALvY1Ahj4Rs+Ubzo46hRjE1Na/j4WUNXBACwhU2NMPYF9uLYlJFzft33G/3yjh8bOQsAYIeW2x7RZT/aFPuc4thk7DOMha+5YOaouzp/rl/0P2HkLACAHf763w7q34dOxD6nudAU+wxj4Wtb0azF+ROxpmwh36C2lZebuiQAQIYmJib06quvqq+vT/v+e1KFn9wmNUYPl6lGGHuPb/vGVbHPKEnafn38cwAA2ZiYmND+/ft1zz33aOXKldq9e7c2bNigwV89oUWLFsU621QjjIVv+dLFuvHqFuVyEQ8oTeunP7yCG1cDgGOqxe6DDz7Q22+/rQcffFCffPCezn3xvqJ+nyGXk25a12KkEUa/wH5/R6sK+caIFzKt//yXv9Xg4KDJSwIAJGCh2D300ENavXq1Tp48qc7OTj322GP6p7+4WUsWRXuHrZBvVHdHq5FrNxq+61Yv066tbVrSVN+xS5oa9MSfXaeepx/Rzp071dXVpWKxaPLSAAAx1Ro7SSqVSurt7VV7e7vWrl2rw4cP664//WnkRuza2qb2VcuM/BzWPZ3hm2++0cMPP6zBwUH19PRoy5Ytpi8PAFCjyg+oDAwMaN26ders7NT27dsvRG62kydPqru7Wx999JGef/55bdp08dcYvHs6Q9nQyBk9e+CoXv/0lHKa+eJhWflZSzeta1F3R+ucFR8cHNTOnTt16623as+ePWpubk7iMgEAs0SJnTSz8l5++WU98MADuvvuu/X444+rUCjM+XfjNiKOxMJX9vXouPoPjWj4+FkVxybVXGhS28rLtf36hZ+uy/oDgHREjV3ZQiuvmjiNiCrx8JnA+gMA8+LGTqpv5dnCifBJrD8AMMFE7MqirrysORO+MtYfANTHZOwkN1deJefCJ7H+AGAhpmNX5urKq+Rk+MpYfwDwvaRiJ7m/8io5HT6J9QcgbEnGrsyHlVfJ+fCVsf4AhCKN2El+rbxK3oRPYv0B8FdasSvzbeVV8ip8Zaw/AD5IO3aSvyuvkpfhk1h/ANyURezKfF55lbwNXxnrD4DtsoydFMbKq+R9+CTWHwD7ZB27slBWXqUgwlfG+gOQJVtiJ4W38ioFFT6J9QcgXTbFrizElVcpuPCVsf4AJMXG2Elhr7xKwYZPYv0BMMfW2JWFvvIqBR2+MtYfgChsj53EypsL4fsO6w9ALVyIXRkrb26EbxbWH4DZXIqdxMpbCOGbA+sPgGuxK2PlLYzwzYP1B4TF1dhJrLx6EL4FsP4Av7kcuzJWXn0IX41Yf4A/fIidxMqLivDVgfUHuMuX2JWx8qIjfBGw/gA3+BY7iZVnAuGLiPUH2MnH2JWx8swgfDGx/oDs+Rw7iZVnGuEzgPUHpM/32JWx8swjfAax/oBkhRI7iZWXJMJnGOsPMCuk2JWx8pJF+BLC+gOiCzF2EisvLYQvQaw/oHahxq6MlZcewpcC1h8wt9BjJ7HyskD4UsL6A2YQu++x8rJB+FLG+kOIiN3FWHnZInwZYP0hBMRubqy87BG+DLH+4BtiVx0rzx6EL2OsP7iO2C2MlWcXwmcJ1h9cQuxqw8qzE+GzCOsPNiN29WHl2YvwWYj1B1sQu/qx8uxH+CzF+kNWiF10rDw3ED7Lsf6QBmIXDyvPLYTPAaw/JIHYmcHKcw/hcwjrD3ERO3NYee4ifI5h/aFexM48Vp7bCJ+jWH+YD7FLBivPD4TPYaw/VCJ2yWLl+YPweYD1Fy5ilzxWnn8InydYf+Egdulh5fmJ8HmG9ecnYpcuVp7fCJ+HWH9+IHbZYOX5j/B5jPXnHmKXHVZeOAif51h/9iN22WPlhYXwBYL1ZxdiZwdWXpgIX0BYf9maK3a33367brvtNmKXAVZeuAhfgFh/6SF29mHlgfAFivWXHGJnL1YeJMIXPNafGcTObqw8VCJ8YP1FROzcwMrDbIQPF7D+Fkbs3MHKQzWEDxdh/V2K2LmHlYf5ED7MKfT1R+zcxMpDLQgfqgpt/RE7t7HyUCvChwX5vP6InftYeagX4UNNfFp/xM4frDxEQfhQF1fXH7HzCysPcRA+1M2V9Ufs/MTKQ1yED5HZuP6Inb9YeTCF8CEWG9YfsfMfKw8mET4Ykfb6I3ZhYOUhCYQPxiS9/ohdWFh5SArhg3Em1x+xCw8rD0kjfEhEnPVH7MLFykMaCB8SVev6I3ZhY+UhTYQPiau2/ogdJFYe0kf4kJrBwUHde++9uuaaa7R8+XLt27eP2AWMlYes5LO+APivctmNjo7q0KFDmpiY0N69e7Vjx46sLw8ZqFx5AwMDrDykqiHrC4CfJiYmtH//ft1zzz1auXKldu/erfb2dh05ckRfffWVent7tWvXLnV1dalYLGZ9uUhJqVRSb2+v2tvbtXbtWh0+fJjoIXX8qhPG1PuenQ13fUF6eC8PtiB8iMXEB1RsvOcnzOG9PNiG8KFuSXwak/XnJ1YebET4UJO0vnrA+vMDKw82I3yoKqvv2bH+3MbKg+0IHy5i05fKWX9uYeXBFYQPVsVuNtafG1h5cAnhC5TNsZsL689OrDy4iPAFxLXYzcb6swsrD64ifJ5zPXZzYf1li5UH1xE+D/kYu9lYf9lg5cEHhM8TIcRuLqy/dLDy4BPC57BQYzcb6y9ZrDz4hvA5hthVx/ozi5UHXxE+BxC72rH+zGDlwWeEz1LELh7WXzSsPISA8FmE2JnF+qsPKw+hIHwZI3bJY/3Nj5WH0BC+DBC79LH+5sbKQ4gIX0qInR1YfzNYeQgZ4UsQsbNT6OuPlYfQET7DiJ07Qlt/rDxgBuEzgNi5K5T1x8oDvkf4IiJ2fvF1/bHygEsRvjoQO7/5tv5YecDcCN8CiF14XF9/rDxgfoRvDsQOrq4/Vh6wMML3HWKHubiy/lh5QO2CDh+xQy1sX3+sPKA+wYWP2CEq29YfKw+IJojwETuYYsv6Y+UB0XkbPmKHJGW1/lh5QHxehY/YIU1prz9WHmBG4uE7PTqu/oMjGj5RVHFsSs2FvNpWNKtz4ypduXRx7POJHbKW9Ppj5QFmJRa+I8fO6JkDR/XGZ6ckSeNT0xf+rJBvUElSx7oWdd/YqutWL6vrbGIH2yS1/lh5gHmJhO+ld77UU/uGNTZ1XvOdnstJhXyjdm1t047Na+Y9k9jBBabWHysPSI7x8M1E7xOdm5xe+C9/Z0lTg3ZtXX9J/IgdXBR3/bHygGQZDd+RY2d0R887Ojd5vu7XLmlqVG/XZrX93mXEDl6od/2x8oB0GA1f14vv63effDXvrzerXoik5WP/qy/+9R+IHbxR6/pj5QHpMRa+06PjuuHp1y76EEu9GjWtV/7yWrVf/QcmLgmwRrX1x8oD0tdg6qD+gyOxz2jK5/VfJ6KHE7DVli1bNDQ0pFKppA0bNmhwcFAnT55UZ2enHnvsMQ0MDOjpp58mekAKjIVv+EQx1tqTpLGpaQ0fP2voigC7/OAHP1BPT4+ee+453XXXXVqzZo1WrVqlw4cP86tNIEXGwlccmzJ0zqSRcwAbnTx5Uj09Pbriiit0yy236JVXXtGbb76Z9WUBQTEWvuZC3tA5TUbOAWxSKpXU29ur9vZ2rV27VkNDQxoYGFBPT4927typrq4uFYvFrC8TCIKx8LWtaNbifMzjpib0P0Nv68MPP5RHtxBF4OZ7L2+u9/4AJMtY+LZvXBX7jKZFi/T748e0bds2rV+/Xo8++igRhLNmr7xq7+WV3/tj/QHpMBa+5UsX68arW5TLRXt9Lif9yfqrtHfPP+rLL7/UCy+8oG+//ZYIwklRPrHJ+gPSYSx8knR/R6sK+cZIry3kG9Xd0SpJyuVy2rRpk/bs2UME4ZRaV141rD8geVbfq3O2Uqmk9957T319ferr69OSJUt0++23q7OzU9dee61yUecmYIDpu6/Y8rR3wDfOPJ1hNiIIWyR995WsnvYO+Cqx5/ENjZzRsweO6vVPTymnmS+nl5Wfx3fTuhZ1d7SqfdWyWP8WEURW0rrHJusPMCfxJ7B/PTqu/kMjGj5+VsWxSTUXmtS28nJtv97ME9hnI4JIQ1b32GT9AfElHr4sEUEkIesnKbD+gHi8Dl8lIoi4bHuSAusPiCaY8FUigqhX1iuvGtYfUL8gw1eJCGI+tq28alh/QO2CD18lIohKtq68alh/QG0IXxVEMFyurLxqWH/A/AhfDYhgOFxbedWw/oDqCF+diKCfXF951bD+gEsRvhiIoB98WXnVsP6AixE+Q4ige3xdedWw/oAZhC8BRNB+vq+8alh/AOFLHBG0S2grrxrWH0JG+FJEBLMV6sqrhvWHUBG+jBDB9LDy5sf6Q2gInwWIYHJYebVh/SEkhM8yRNAMVl40rD+EgPBZjAhGw8qLh/UH3xE+RxDBhbHyzGL9wVeEz0FE8FKsvGSw/uAjwue40CPIyksH6w8+IXweCS2CrLx0sf7gC8LnKZ8jyMrLFusPriN8AfApgqw8O7D+4DLCFxhXI8jKsxPrDy4ifAFzJYKsPLux/uAawgdJdkaQlecW1h9cQfhwCRsiyMpzE+sPLiB8mFfaEWTl+YH1B5sRPtQs6Qiy8vzC+oOtCB8iMRlBVp7fWH+wDeFDbHEiyMoLA+sPNiF8MKrWCLLywsT6gw0IHxJTLYI333yz9u7dq48//piVFyDWH7JG+JCKUqmkd999V08++aT279+vZcuW6b777tOdd95p1ZflkR7WH7LSkPUFIAynTp3Snj179Pnnn+utt97Svn37NDExoW3btmn9+vV69NFH9eGHH4r/h4Vjy5YtGhoaUqlU0oYNGzQ4OJj1JSEQLD4kaqH38mz4sjyyx/pDmggfElPvJzaJYNh47w9pIXwwzsQnNolguFh/SBrhg1FJfC+PCIaH9YckET4Ykdb38ohgWFh/SALhQ2xZ3X2FCIaB9QfTCB8is+nuK0TQf6w/mEL4EInN99gkgv5i/cEEwoe62LTyakEE/cT6QxyEDzWzeeXVggj6hfWHqAgfFuTayqsFEfQH6w/1InyYl+srrxZE0H2sP9SD8GFOPq68WhBBt7H+UAvCh0uEsPJqQQTdxPrDQggfLgh15dWCCLqH9YdqCB8ksfLqQQTdwfrDXAhf4Fh58RBBN7D+UInwBYyVZxYRtBvrD2WEL0CsvOQRQXux/kD4AsPKSx8RtA/rL2yELxCsPDsQQbuw/sJE+ALAyrMTEbQD6y88hM9jrDx3EMHssf7CQfg8xcpzFxHMDusvDITPM6w8vxDBbLD+/Eb4PMLK8xsRTBfrz1+EzwOsvPAQwfSw/vxD+BzHygMRTB7rzy+Ez1GsPMyFCCaL9ecHwucgVh5qQQSTwfpzH+FzCCsPURFB81h/7iJ8jmDlwRQiaA7rz02Ez3KsPCSJCJrB+nML4bMYKw9pIoLxsP7cQfgsxMpD1ohgdKw/+xE+y7DyYBsiWD/Wn90InyVYeXABEawP689OhM8CrDy4iAjWhvVnH8KXIVYefEEEF8b6swfhywgrD74igtWx/uxA+FLGykNIiODcWH/ZInwpYuUhZETwYqy/7BC+FLDygIsRwe+x/tJH+BLGygPmRwRZf2kjfAlh5QH1Cz2CrL90EL4EsPKA+EKNIOsveYTPIFYekIwQI8j6Sw7hM4SVB6QjpAiy/pJB+GJi5QHZCSWCrD+zCF8MrDzAHr5HkPVnDuGLgJUH2M3nCLL+4iN8dWLlAW7xMYKsv3gIX41YeYD7fIsg6y8awlcDVh7gH18iyPqrH+GbBysPCIMPEWT91Y7wVcHKA8LkcgRZf7UhfLOw8gCUuRpB1t/8CF8FVh6AalyLIOuvOsInVh6A+rgUQdbfpYIPHysPQBwuRJD1d7Fgw8fKA2Ca7RFk/c0IMnysPABJszWCrL/AwsfKA5AFGyMY8voLJnysPAA2sCmCoa4/78PHygNgK1siGNr68zp8rDwArsg6giGtPy/Dx8oD4LIsIxjC+vMufKw8AD7JIoK+rz9vwsfKA+C7tCPo6/rzInysPAChSSuCPq4/p8PHygOAdCLo0/pzNnysPAC4VJIR9GX9ORc+Vh4A1CapCLq+/pwKHysPAKIxHUGX158T4WPlAYA5JiPo4vqzPnysPABIjokIurb+rA0fKw8A0hU3gq6sv8TDd3p0XP0HRzR8oqji2JSaC3m1rWhW58ZVunLp4jlfw8oDgGxFjWC96y9KI+JKLHxHjp3RMweO6o3PTkmSxqemL/xZId+gkqSOdS3qvrFV161eJomVBwA2ihLBhdZflEaYkkj4XnrnSz21b1hjU+c13+m5nFTIN2rX1jZt+eFlrDwAsFw9Eay2/qI0YsfmNcZ+BuPhm/mBPtG5yemF//J3mnIlffvWi9qxeQ0rDwAcUWsEK9ffT/78b/TPr31RVyOWNDVo19b1xuJnNHxHjp3RHT3v6Nzk+bpfu6hR6v+rG9S+apmpywEApGShCBaLRe38+yf17tI/Vi5f/3t3S5oa1du12UgjjIav68X39btPvpp3ula9kJx06x9epV/t+CNTlwMAyEC1CH52VYfeHTmnKNEx2Qhj4Ts9Oq4bnn7tojco67U436C3/+7mxD7JAwBIVzmCL/b9Vv/RsEm5xqbIZ5lqREOsV1foPzgS+4ycpP5D8c8BANghl8tp06ZNuvZn96qwOF6wTDXCWPiGTxRjrT1JGpua1vDxs4auCABgC5sakY99wneKY1NGzvl132/0yzt+bOQsAIAdWm57RJf9KP7X1Ipjk7HPMBa+5oKZo+7q/Ll+0f+EkbMAAHZ4sPewfvvB/8U+p7kQ/T3CMmO/6mxb0azF+XjHFfINalt5uaErAgDYwqZGGAvf9o2rYp9RkrT9+vjnAADsYlMjjIVv+dLFuvHqFkV9qn0uJ920roWvMgCAh2xqhLHwSdL9Ha0q5BsjvbaQb1R3R6vJywEAWMSWRhgN33Wrl2nX1jYtaarv2Jn7sLVxuzIA8JgtjTD2qc6y8k1Es7zzNgDATjY0IrHn8Q2NnNGzB47q9U9PKaeZLx6WlZ+1dNO6FnV3tLL0ACAwWTYi8Sewfz06rv5DIxo+flbFsUk1F5rUtvJybb8+uafrAgDckEUjEg8fAAA2MfrhFgAAbEf4AABBIXwAgKAQPgBAUAgfACAohA8AEBTCBwAICuEDAASF8AEAgvL/11epqUdzx6AAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "G = nx.Graph()\n",
    "G.add_edges_from([[0,3],[0,4],[1,3],[1,4],[2,3],[2,4]])\n",
    "nx.draw(G, pos=nx.bipartite_layout(G, [0,1,2]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### We can now build the circuit for cost operator\n",
    "\n",
    "Cost operator $ e^{-i\\gamma_p C} = e^{-i\\gamma_p \\sum_{ij\\in E}(-Z_iZ_j)}$\n",
    "\n",
    "Note that later we will see how to do this more efficiently in Qiskit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "def append_zz_term(qc, q1, q2, gamma):\n",
    "    qc.cx(q1,q2)\n",
    "    qc.rz(2*gamma, q2)\n",
    "    qc.cx(q1,q2)\n",
    "\n",
    "def get_cost_operator_circuit(G, gamma):\n",
    "    N = G.number_of_nodes()\n",
    "    qc = QuantumCircuit(N,N)\n",
    "    for i, j in G.edges():\n",
    "        append_zz_term(qc, i, j, gamma)\n",
    "    return qc\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Let's verify that our function produces the correct circuit\n",
    "\n",
    "Cost operator $e^{-i\\gamma_p C} = e^{-i\\gamma_p \\sum_{ij\\in E}(-Z_iZ_j)}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<pre style=\"word-wrap: normal;white-space: pre;background: #fff0;line-height: 1.1;font-family: &quot;Courier New&quot;,Courier,monospace\">                                                                               »\n",
       "q_0: ──■─────────────────■────■──────────────────────────────■─────────────────»\n",
       "       │                 │    │      ┌───┐    ┌───────────┐  │  ┌───┐     ┌───┐»\n",
       "q_1: ──┼─────────────────┼────┼──────┤ X ├────┤ RZ(2pi/3) ├──┼──┤ X ├─────┤ X ├»\n",
       "       │                 │    │      └─┬─┘    └───────────┘  │  └─┬─┘┌───┐└─┬─┘»\n",
       "q_2: ──┼─────────────────┼────┼────────┼─────────────────────┼────┼──┤ X ├──┼──»\n",
       "     ┌─┴─┐┌───────────┐┌─┴─┐  │        │                     │    │  └─┬─┘  │  »\n",
       "q_3: ┤ X ├┤ RZ(2pi/3) ├┤ X ├──┼────────■─────────────────────┼────■────■────┼──»\n",
       "     └───┘└───────────┘└───┘┌─┴─┐┌───────────┐             ┌─┴─┐            │  »\n",
       "q_4: ───────────────────────┤ X ├┤ RZ(2pi/3) ├─────────────┤ X ├────────────■──»\n",
       "                            └───┘└───────────┘             └───┘               »\n",
       "c: 5/══════════════════════════════════════════════════════════════════════════»\n",
       "                                                                               »\n",
       "«                                                   \n",
       "«q_0: ──────────────────────────────────────────────\n",
       "«     ┌───────────┐     ┌───┐                       \n",
       "«q_1: ┤ RZ(2pi/3) ├─────┤ X ├───────────────────────\n",
       "«     ├───────────┤┌───┐└─┬─┘┌───┐┌───────────┐┌───┐\n",
       "«q_2: ┤ RZ(2pi/3) ├┤ X ├──┼──┤ X ├┤ RZ(2pi/3) ├┤ X ├\n",
       "«     └───────────┘└─┬─┘  │  └─┬─┘└───────────┘└─┬─┘\n",
       "«q_3: ───────────────■────┼────┼─────────────────┼──\n",
       "«                         │    │                 │  \n",
       "«q_4: ────────────────────■────■─────────────────■──\n",
       "«                                                   \n",
       "«c: 5/══════════════════════════════════════════════\n",
       "«                                                   </pre>"
      ],
      "text/plain": [
       "                                                                               »\n",
       "q_0: ──■─────────────────■────■──────────────────────────────■─────────────────»\n",
       "       │                 │    │      ┌───┐    ┌───────────┐  │  ┌───┐     ┌───┐»\n",
       "q_1: ──┼─────────────────┼────┼──────┤ X ├────┤ RZ(2pi/3) ├──┼──┤ X ├─────┤ X ├»\n",
       "       │                 │    │      └─┬─┘    └───────────┘  │  └─┬─┘┌───┐└─┬─┘»\n",
       "q_2: ──┼─────────────────┼────┼────────┼─────────────────────┼────┼──┤ X ├──┼──»\n",
       "     ┌─┴─┐┌───────────┐┌─┴─┐  │        │                     │    │  └─┬─┘  │  »\n",
       "q_3: ┤ X ├┤ RZ(2pi/3) ├┤ X ├──┼────────■─────────────────────┼────■────■────┼──»\n",
       "     └───┘└───────────┘└───┘┌─┴─┐┌───────────┐             ┌─┴─┐            │  »\n",
       "q_4: ───────────────────────┤ X ├┤ RZ(2pi/3) ├─────────────┤ X ├────────────■──»\n",
       "                            └───┘└───────────┘             └───┘               »\n",
       "c: 5/══════════════════════════════════════════════════════════════════════════»\n",
       "                                                                               »\n",
       "«                                                   \n",
       "«q_0: ──────────────────────────────────────────────\n",
       "«     ┌───────────┐     ┌───┐                       \n",
       "«q_1: ┤ RZ(2pi/3) ├─────┤ X ├───────────────────────\n",
       "«     ├───────────┤┌───┐└─┬─┘┌───┐┌───────────┐┌───┐\n",
       "«q_2: ┤ RZ(2pi/3) ├┤ X ├──┼──┤ X ├┤ RZ(2pi/3) ├┤ X ├\n",
       "«     └───────────┘└─┬─┘  │  └─┬─┘└───────────┘└─┬─┘\n",
       "«q_3: ───────────────■────┼────┼─────────────────┼──\n",
       "«                         │    │                 │  \n",
       "«q_4: ────────────────────■────■─────────────────■──\n",
       "«                                                   \n",
       "«c: 5/══════════════════════════════════════════════\n",
       "«                                                   "
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "qc = get_cost_operator_circuit(G, np.pi / 3)\n",
    "qc.draw()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Now we need to build mixer operator\n",
    "\n",
    "Mixer operator: $e^{-i\\beta_p B} = e^{-i\\beta_p \\sum_{i\\in V}X_i} = e^{-i\\beta_p X_0}e^{-i\\beta_p X_1}e^{-i\\beta_p X_2}e^{-i\\beta_p X_3}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [],
   "source": [
    "def append_x_term(qc, q1, beta):\n",
    "    qc.rx(2*beta, q1)\n",
    "\n",
    "def get_mixer_operator_circuit(G, beta):\n",
    "    N = G.number_of_nodes()\n",
    "    qc = QuantumCircuit(N,N)\n",
    "    for n in G.nodes():\n",
    "        append_x_term(qc, n, beta)\n",
    "    return qc\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Let's verify the circuit\n",
    "\n",
    "Mixer operator: $e^{-i\\beta_p B} = e^{-i\\beta_p \\sum_{i\\in V}X_i} = e^{-i\\beta_p X_0}e^{-i\\beta_p X_1}e^{-i\\beta_p X_2}e^{-i\\beta_p X_3}e^{-i\\beta_p X_4}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<pre style=\"word-wrap: normal;white-space: pre;background: #fff0;line-height: 1.1;font-family: &quot;Courier New&quot;,Courier,monospace\">     ┌───────────┐\n",
       "q_0: ┤ RX(2pi/3) ├\n",
       "     ├───────────┤\n",
       "q_1: ┤ RX(2pi/3) ├\n",
       "     ├───────────┤\n",
       "q_2: ┤ RX(2pi/3) ├\n",
       "     ├───────────┤\n",
       "q_3: ┤ RX(2pi/3) ├\n",
       "     ├───────────┤\n",
       "q_4: ┤ RX(2pi/3) ├\n",
       "     └───────────┘\n",
       "c: 5/═════════════\n",
       "                  </pre>"
      ],
      "text/plain": [
       "     ┌───────────┐\n",
       "q_0: ┤ RX(2pi/3) ├\n",
       "     ├───────────┤\n",
       "q_1: ┤ RX(2pi/3) ├\n",
       "     ├───────────┤\n",
       "q_2: ┤ RX(2pi/3) ├\n",
       "     ├───────────┤\n",
       "q_3: ┤ RX(2pi/3) ├\n",
       "     ├───────────┤\n",
       "q_4: ┤ RX(2pi/3) ├\n",
       "     └───────────┘\n",
       "c: 5/═════════════\n",
       "                  "
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "qc = get_mixer_operator_circuit(G, np.pi/3)\n",
    "qc.draw()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Now we are ready to build the full circuit\n",
    "\n",
    "$|\\psi{(\\mathbb{\\beta},\\mathbb{\\gamma})}> =  e^{-i\\beta_p B}e^{-i\\gamma_p C}\\cdots e^{-i\\beta_1 B}e^{-i\\gamma_1 C} H^{\\otimes n} |0>$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_qaoa_circuit(G, beta, gamma):\n",
    "    assert(len(beta) == len(gamma))\n",
    "    p = len(beta) # infering number of QAOA steps from the parameters passed\n",
    "    N = G.number_of_nodes()\n",
    "    qc = QuantumCircuit(N,N)\n",
    "    # first, apply a layer of Hadamards\n",
    "    qc.h(range(N))\n",
    "    # second, apply p alternating operators\n",
    "    for i in range(p):\n",
    "        qc += get_cost_operator_circuit(G,gamma[i])\n",
    "        qc += get_mixer_operator_circuit(G,beta[i])\n",
    "    # finally, do not forget to measure the result!\n",
    "    qc.barrier(range(N))\n",
    "    qc.measure(range(N), range(N))\n",
    "    return qc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Let's verify that circuit looks right\n",
    "\n",
    "$|\\psi{(\\mathbb{\\beta},\\mathbb{\\gamma})}> =  e^{-i\\beta_p B}e^{-i\\gamma_p C}\\cdots e^{-i\\beta_1 B}e^{-i\\gamma_1 C} H^{\\otimes n} |0>$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<pre style=\"word-wrap: normal;white-space: pre;background: #fff0;line-height: 1.1;font-family: &quot;Courier New&quot;,Courier,monospace\">     ┌───┐                                                  ┌───────────┐     »\n",
       "q_0: ┤ H ├──■──────────────■────■────────────────────────■──┤ RX(2pi/3) ├─────»\n",
       "     ├───┤  │              │    │    ┌───┐   ┌────────┐  │  └───┬───┬───┘     »\n",
       "q_1: ┤ H ├──┼──────────────┼────┼────┤ X ├───┤ RZ(pi) ├──┼──────┤ X ├─────────»\n",
       "     ├───┤  │              │    │    └─┬─┘   └────────┘  │      └─┬─┘    ┌───┐»\n",
       "q_2: ┤ H ├──┼──────────────┼────┼──────┼─────────────────┼────────┼──────┤ X ├»\n",
       "     ├───┤┌─┴─┐┌────────┐┌─┴─┐  │      │                 │        │      └─┬─┘»\n",
       "q_3: ┤ H ├┤ X ├┤ RZ(pi) ├┤ X ├──┼──────■─────────────────┼────────■────────■──»\n",
       "     ├───┤└───┘└────────┘└───┘┌─┴─┐┌────────┐          ┌─┴─┐                  »\n",
       "q_4: ┤ H ├────────────────────┤ X ├┤ RZ(pi) ├──────────┤ X ├──────────────────»\n",
       "     └───┘                    └───┘└────────┘          └───┘                  »\n",
       "c: 5/═════════════════════════════════════════════════════════════════════════»\n",
       "                                                                              »\n",
       "«                                                                           ░ »\n",
       "«q_0: ──────────────────────────────────────────────────────────────────────░─»\n",
       "«     ┌───┐┌────────┐     ┌───┐┌───────────┐                                ░ »\n",
       "«q_1: ┤ X ├┤ RZ(pi) ├─────┤ X ├┤ RX(2pi/3) ├────────────────────────────────░─»\n",
       "«     └─┬─┘├────────┤┌───┐└─┬─┘└───┬───┬───┘  ┌────────┐ ┌───┐┌───────────┐ ░ »\n",
       "«q_2: ──┼──┤ RZ(pi) ├┤ X ├──┼──────┤ X ├──────┤ RZ(pi) ├─┤ X ├┤ RX(2pi/3) ├─░─»\n",
       "«       │  └────────┘└─┬─┘  │      └─┬─┘    ┌─┴────────┴┐└─┬─┘└───────────┘ ░ »\n",
       "«q_3: ──┼──────────────■────┼────────┼──────┤ RX(2pi/3) ├──┼────────────────░─»\n",
       "«       │                   │        │      └───────────┘  │  ┌───────────┐ ░ »\n",
       "«q_4: ──■───────────────────■────────■─────────────────────■──┤ RX(2pi/3) ├─░─»\n",
       "«                                                             └───────────┘ ░ »\n",
       "«c: 5/════════════════════════════════════════════════════════════════════════»\n",
       "«                                                                             »\n",
       "«     ┌─┐            \n",
       "«q_0: ┤M├────────────\n",
       "«     └╥┘┌─┐         \n",
       "«q_1: ─╫─┤M├─────────\n",
       "«      ║ └╥┘┌─┐      \n",
       "«q_2: ─╫──╫─┤M├──────\n",
       "«      ║  ║ └╥┘┌─┐   \n",
       "«q_3: ─╫──╫──╫─┤M├───\n",
       "«      ║  ║  ║ └╥┘┌─┐\n",
       "«q_4: ─╫──╫──╫──╫─┤M├\n",
       "«      ║  ║  ║  ║ └╥┘\n",
       "«c: 5/═╩══╩══╩══╩══╩═\n",
       "«      0  1  2  3  4 </pre>"
      ],
      "text/plain": [
       "     ┌───┐                                                  ┌───────────┐     »\n",
       "q_0: ┤ H ├──■──────────────■────■────────────────────────■──┤ RX(2pi/3) ├─────»\n",
       "     ├───┤  │              │    │    ┌───┐   ┌────────┐  │  └───┬───┬───┘     »\n",
       "q_1: ┤ H ├──┼──────────────┼────┼────┤ X ├───┤ RZ(pi) ├──┼──────┤ X ├─────────»\n",
       "     ├───┤  │              │    │    └─┬─┘   └────────┘  │      └─┬─┘    ┌───┐»\n",
       "q_2: ┤ H ├──┼──────────────┼────┼──────┼─────────────────┼────────┼──────┤ X ├»\n",
       "     ├───┤┌─┴─┐┌────────┐┌─┴─┐  │      │                 │        │      └─┬─┘»\n",
       "q_3: ┤ H ├┤ X ├┤ RZ(pi) ├┤ X ├──┼──────■─────────────────┼────────■────────■──»\n",
       "     ├───┤└───┘└────────┘└───┘┌─┴─┐┌────────┐          ┌─┴─┐                  »\n",
       "q_4: ┤ H ├────────────────────┤ X ├┤ RZ(pi) ├──────────┤ X ├──────────────────»\n",
       "     └───┘                    └───┘└────────┘          └───┘                  »\n",
       "c: 5/═════════════════════════════════════════════════════════════════════════»\n",
       "                                                                              »\n",
       "«                                                                           ░ »\n",
       "«q_0: ──────────────────────────────────────────────────────────────────────░─»\n",
       "«     ┌───┐┌────────┐     ┌───┐┌───────────┐                                ░ »\n",
       "«q_1: ┤ X ├┤ RZ(pi) ├─────┤ X ├┤ RX(2pi/3) ├────────────────────────────────░─»\n",
       "«     └─┬─┘├────────┤┌───┐└─┬─┘└───┬───┬───┘  ┌────────┐ ┌───┐┌───────────┐ ░ »\n",
       "«q_2: ──┼──┤ RZ(pi) ├┤ X ├──┼──────┤ X ├──────┤ RZ(pi) ├─┤ X ├┤ RX(2pi/3) ├─░─»\n",
       "«       │  └────────┘└─┬─┘  │      └─┬─┘    ┌─┴────────┴┐└─┬─┘└───────────┘ ░ »\n",
       "«q_3: ──┼──────────────■────┼────────┼──────┤ RX(2pi/3) ├──┼────────────────░─»\n",
       "«       │                   │        │      └───────────┘  │  ┌───────────┐ ░ »\n",
       "«q_4: ──■───────────────────■────────■─────────────────────■──┤ RX(2pi/3) ├─░─»\n",
       "«                                                             └───────────┘ ░ »\n",
       "«c: 5/════════════════════════════════════════════════════════════════════════»\n",
       "«                                                                             »\n",
       "«     ┌─┐            \n",
       "«q_0: ┤M├────────────\n",
       "«     └╥┘┌─┐         \n",
       "«q_1: ─╫─┤M├─────────\n",
       "«      ║ └╥┘┌─┐      \n",
       "«q_2: ─╫──╫─┤M├──────\n",
       "«      ║  ║ └╥┘┌─┐   \n",
       "«q_3: ─╫──╫──╫─┤M├───\n",
       "«      ║  ║  ║ └╥┘┌─┐\n",
       "«q_4: ─╫──╫──╫──╫─┤M├\n",
       "«      ║  ║  ║  ║ └╥┘\n",
       "«c: 5/═╩══╩══╩══╩══╩═\n",
       "«      0  1  2  3  4 "
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "qc = get_qaoa_circuit(G,[np.pi/3], [np.pi/2])\n",
    "qc.draw()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Running a circuit in Qiskit\n",
    "\n",
    "For the purposes of this tutorial, we will be running the circuits in simulator. To run your circuits on real IBM quantum computers, all you need to do is change the backend. However, the devices are noisy, so there might be dragons!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "#### Qiskit qubit ordering\n",
    "\n",
    "Note that Qiskit uses an ordering where zeroth qubit is the rightmost (the least significant bit in the bistring),\n",
    "so we need to invert the samples we received"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def invert_counts(counts):\n",
    "    return {k[::-1]:v for k, v in counts.items()}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Now we can run the circuit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'00000': 31, '10000': 34, '00001': 28, '10001': 34, '01001': 24, '11001': 34, '00101': 42, '10101': 27, '01101': 34, '11101': 35, '00011': 27, '10011': 30, '01011': 42, '11011': 38, '00111': 30, '10111': 25, '01111': 34, '11111': 34, '01000': 35, '11000': 34, '00100': 31, '10100': 35, '01100': 26, '11100': 37, '00010': 31, '10010': 27, '01010': 30, '11010': 36, '00110': 41, '10110': 28, '01110': 29, '11110': 21}\n"
     ]
    }
   ],
   "source": [
    "backend = Aer.get_backend('qasm_simulator')\n",
    "job = execute(qc, backend)\n",
    "result = job.result()\n",
    "print(invert_counts(result.get_counts()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### How good were these results?\n",
    "\n",
    "To answer that, we need to compute $<\\psi|C|\\psi>$. Luckily, C is a classical Hamiltonian, so we can do this easily. We are going to compute -1 times the value of the cut because by convention classical optimizers minimize (this will come up shortly)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-6"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def maxcut_obj(x,G):\n",
    "    cut = 0\n",
    "    for i, j in G.edges():\n",
    "        if x[i] != x[j]:\n",
    "            # the edge is cut\n",
    "            cut -= 1\n",
    "    return cut\n",
    "\n",
    "maxcut_obj('00011', G)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Estimate energy from the sample"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-2.9970703125"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def compute_maxcut_energy(counts, G):\n",
    "    energy = 0\n",
    "    total_counts = 0\n",
    "    for meas, meas_count in counts.items():\n",
    "        obj_for_meas = maxcut_obj(meas, G)\n",
    "        energy += obj_for_meas * meas_count\n",
    "        total_counts += meas_count\n",
    "    return energy / total_counts\n",
    "\n",
    "compute_maxcut_energy(invert_counts(result.get_counts()), G)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Putting all the pieces together\n",
    "\n",
    "As you remember, we run QAOA as a variational algorithm, so we will need a classical optimizer to optimize the parameters $\\beta$, $\\gamma$. We will use Constrained Optimization BY Linear Approximation (COBYLA) from `scipy.optimize`.\n",
    "\n",
    "We need to wrap our circuit into a black-box objective function `f` that COBYLA could optimize. Note that by convention all optimizers perform _minimization_, so we had to put a minus sign in front of our objective earlier."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_black_box_objective(G,p):\n",
    "    backend = Aer.get_backend('qasm_simulator')\n",
    "    def f(theta):\n",
    "        # let's assume first half is betas, second half is gammas\n",
    "        beta = theta[:p]\n",
    "        gamma = theta[p:]\n",
    "        qc = get_qaoa_circuit(G,beta, gamma)\n",
    "        counts = execute(qc, backend, seed_simulator=10).result().get_counts()\n",
    "        # return the energy\n",
    "        return compute_maxcut_energy(invert_counts(counts), G)\n",
    "    return f"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Now we can run the optimizer to get optimal parameters!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "     fun: -5.2978515625\n",
       "   maxcv: 0.0\n",
       " message: 'Optimization terminated successfully.'\n",
       "    nfev: 134\n",
       "  status: 1\n",
       " success: True\n",
       "       x: array([1.83749394, 2.1345936 , 1.9319531 , 2.09862998, 0.347631  ,\n",
       "       2.04618358, 6.11249139, 2.31785137, 3.56479392, 5.87260613])"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# p is the number of QAOA alternating operators\n",
    "p = 5\n",
    "obj = get_black_box_objective(G, p)\n",
    "\n",
    "init_point = np.array([0.81069872, 2.2067517 , 0.83830696, 2.15579759, 0.37060699,\n",
    "       2.42068091, 6.1575306 , 2.2453419 , 3.85060091, 6.137845  ])\n",
    "# We are going to limit the number of iterations to 2500\n",
    "res_sample = minimize(obj, init_point, method='COBYLA', options={'maxiter':2500, 'disp': True})\n",
    "res_sample"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### To get the solution strings, let's run the circuit with optimal parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'10000': 7,\n",
       " '00001': 32,\n",
       " '10001': 1,\n",
       " '01001': 3,\n",
       " '11001': 1,\n",
       " '10101': 1,\n",
       " '01101': 3,\n",
       " '11101': 24,\n",
       " '00011': 363,\n",
       " '10011': 21,\n",
       " '01011': 16,\n",
       " '11011': 2,\n",
       " '00111': 26,\n",
       " '10111': 5,\n",
       " '01111': 4,\n",
       " '11111': 1,\n",
       " '01000': 1,\n",
       " '11000': 32,\n",
       " '00100': 3,\n",
       " '10100': 23,\n",
       " '01100': 24,\n",
       " '11100': 372,\n",
       " '00010': 28,\n",
       " '10010': 2,\n",
       " '01010': 1,\n",
       " '11010': 2,\n",
       " '00110': 3,\n",
       " '01110': 3,\n",
       " '11110': 20}"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "optimal_theta = res_sample['x']\n",
    "qc = get_qaoa_circuit(G, optimal_theta[:p], optimal_theta[p:])\n",
    "counts = invert_counts(execute(qc, backend).result().get_counts())\n",
    "counts"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### How good are there results?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<BarContainer object of 5 artists>"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAQ00lEQVR4nO3dcYzfd13H8eeLDoaihs5dZ20bW5KCtggDj4ohGqHACiPrMC4pUdLoTJFsBBINtpCIhjRpQEWjjqTCtInAUhFswxSpVSQmsnKDAeu2usrGerSsB4agkpR0vP3jvsNf27v7fa9311/72fORXL7f7+f3+f5+r7vkXve97/2+30tVIUlqy9NGHUCStPgsd0lqkOUuSQ2y3CWpQZa7JDXoqlEHALj22mtr7dq1o44hSVeUe++99xtVNTbTY5dFua9du5aJiYlRx5CkK0qSr872mKdlJKlBlrskNchyl6QGWe6S1CDLXZIaZLlLUoMsd0lqkOUuSQ2y3CWpQZfFFaoLtXbn3aOO8H2P7rlx1BEkySN3SWqR5S5JDbLcJalBlrskNchyl6QGDS33JM9Lct/Ax7eTvC3JNUkOJXm4Wy4f2GdXkuNJjiW5YWk/BUnS+YaWe1Udq6rrq+p64GeA7wAfB3YCh6tqPXC42ybJBmAbsBHYAtyRZNnSxJckzWS+p2U2A/9ZVV8FtgL7uvF9wM3d+lbgrqo6U1WPAMeBTYuQVZLU03zLfRvwkW79uqo6BdAtV3Tjq4ATA/tMdmOSpEukd7kneQZwE/A3w6bOMFYzPN+OJBNJJqampvrGkCT1MJ8j99cAn6+qx7vtx5OsBOiWp7vxSWDNwH6rgZPnP1lV7a2q8aoaHxub8Z93S5Iu0nzK/Q38/ykZgIPA9m59O3BgYHxbkquTrAPWA0cWGlSS1F+vG4cl+UHgVcCbBob3APuT3Ao8BtwCUFVHk+wHHgDOArdV1ROLmlqSNKde5V5V3wF+9LyxbzL97pmZ5u8Gdi84nSTponiFqiQ1yHKXpAZZ7pLUIMtdkhpkuUtSgyx3SWqQ5S5JDbLcJalBlrskNchyl6QGWe6S1CDLXZIaZLlLUoMsd0lqkOUuSQ2y3CWpQZa7JDXIcpekBlnuktQgy12SGtSr3JM8O8lHkzyU5MEkP5fkmiSHkjzcLZcPzN+V5HiSY0luWLr4kqSZ9D1y/xPgk1X1k8ALgQeBncDhqloPHO62SbIB2AZsBLYAdyRZttjBJUmzG1ruSX4E+AXggwBV9d2q+hawFdjXTdsH3NytbwXuqqozVfUIcBzYtLixJUlz6XPk/hxgCvjLJF9I8oEkzwKuq6pTAN1yRTd/FXBiYP/JbuwcSXYkmUgyMTU1taBPQpJ0rj7lfhXwYuD9VfUi4H/pTsHMIjOM1QUDVXuraryqxsfGxnqFlST106fcJ4HJqrqn2/4o02X/eJKVAN3y9MD8NQP7rwZOLk5cSVIfQ8u9qr4OnEjyvG5oM/AAcBDY3o1tBw506weBbUmuTrIOWA8cWdTUkqQ5XdVz3luADyV5BvAV4NeY/sGwP8mtwGPALQBVdTTJfqZ/AJwFbquqJxY9uSRpVr3KvaruA8ZneGjzLPN3A7svPpYkaSG8QlWSGmS5S1KDLHdJapDlLkkNstwlqUGWuyQ1yHKXpAZZ7pLUIMtdkhpkuUtSgyx3SWqQ5S5JDbLcJalBlrskNchyl6QGWe6S1CDLXZIaZLlLUoMsd0lqUK9yT/Joki8nuS/JRDd2TZJDSR7ulssH5u9KcjzJsSQ3LFV4SdLM5nPk/vKqur6qnvxH2TuBw1W1HjjcbZNkA7AN2AhsAe5IsmwRM0uShljIaZmtwL5ufR9w88D4XVV1pqoeAY4DmxbwOpKkeepb7gV8Ksm9SXZ0Y9dV1SmAbrmiG18FnBjYd7IbO0eSHUkmkkxMTU1dXHpJ0oyu6jnvZVV1MskK4FCSh+aYmxnG6oKBqr3AXoDx8fELHpckXbxeR+5VdbJbngY+zvRplseTrATolqe76ZPAmoHdVwMnFyuwJGm4oeWe5FlJfvjJdeDVwP3AQWB7N207cKBbPwhsS3J1knXAeuDIYgeXJM2uz2mZ64CPJ3ly/oer6pNJPgfsT3Ir8BhwC0BVHU2yH3gAOAvcVlVPLEl6SdKMhpZ7VX0FeOEM498ENs+yz25g94LTSZIuileoSlKDLHdJapDlLkkNstwlqUGWuyQ1yHKXpAZZ7pLUIMtdkhpkuUtSgyx3SWqQ5S5JDbLcJalBlrskNchyl6QGWe6S1CDLXZIaZLlLUoMsd0lqkOUuSQ3qXe5JliX5QpJPdNvXJDmU5OFuuXxg7q4kx5McS3LDUgSXJM1uPkfubwUeHNjeCRyuqvXA4W6bJBuAbcBGYAtwR5JlixNXktRHr3JPshq4EfjAwPBWYF+3vg+4eWD8rqo6U1WPAMeBTYuSVpLUS98j9z8G3g58b2Dsuqo6BdAtV3Tjq4ATA/Mmu7FzJNmRZCLJxNTU1HxzS5LmMLTck7wOOF1V9/Z8zswwVhcMVO2tqvGqGh8bG+v51JKkPq7qMedlwE1JXgs8E/iRJH8NPJ5kZVWdSrISON3NnwTWDOy/Gji5mKElSXMbeuReVbuqanVVrWX6D6X/XFW/ChwEtnfTtgMHuvWDwLYkVydZB6wHjix6cknSrPocuc9mD7A/ya3AY8AtAFV1NMl+4AHgLHBbVT2x4KSSpN7mVe5V9Wng0936N4HNs8zbDexeYDZJ0kXyClVJapDlLkkNstwlqUGWuyQ1yHKXpAZZ7pLUIMtdkhpkuUtSgyx3SWqQ5S5JDbLcJalBlrskNchyl6QGWe6S1CDLXZIaZLlLUoMsd0lqkOUuSQ2y3CWpQUPLPckzkxxJ8sUkR5P8fjd+TZJDSR7ulssH9tmV5HiSY0luWMpPQJJ0oT5H7meAV1TVC4HrgS1JXgrsBA5X1XrgcLdNkg3ANmAjsAW4I8myJcguSZrF0HKvaf/TbT69+yhgK7CvG98H3NytbwXuqqozVfUIcBzYtJihJUlz63XOPcmyJPcBp4FDVXUPcF1VnQLoliu66auAEwO7T3Zj5z/njiQTSSampqYW8ClIks7Xq9yr6omquh5YDWxK8vw5pmemp5jhOfdW1XhVjY+NjfUKK0nqZ17vlqmqbwGfZvpc+uNJVgJ0y9PdtElgzcBuq4GTCw0qSeqvz7tlxpI8u1v/AeCVwEPAQWB7N207cKBbPwhsS3J1knXAeuDIIueWJM3hqh5zVgL7une8PA3YX1WfSPLvwP4ktwKPAbcAVNXRJPuBB4CzwG1V9cTSxJckzWRouVfVl4AXzTD+TWDzLPvsBnYvOJ0k6aJ4haokNchyl6QGWe6S1CDLXZIaZLlLUoMsd0lqkOUuSQ2y3CWpQZa7JDXIcpekBlnuktQgy12SGmS5S1KDLHdJapDlLkkNstwlqUGWuyQ1yHKXpAZZ7pLUoKHlnmRNkn9J8mCSo0ne2o1fk+RQkoe75fKBfXYlOZ7kWJIblvITkCRdqM+R+1ngt6rqp4CXArcl2QDsBA5X1XrgcLdN99g2YCOwBbgjybKlCC9JmtnQcq+qU1X1+W79v4EHgVXAVmBfN20fcHO3vhW4q6rOVNUjwHFg0yLnliTNYV7n3JOsBV4E3ANcV1WnYPoHALCim7YKODGw22Q3dv5z7UgykWRiamrqIqJLkmbTu9yT/BDwt8Dbqurbc02dYawuGKjaW1XjVTU+NjbWN4YkqYde5Z7k6UwX+4eq6mPd8ONJVnaPrwROd+OTwJqB3VcDJxcnriSpjz7vlgnwQeDBqvqjgYcOAtu79e3AgYHxbUmuTrIOWA8cWbzIkqRhruox52XAG4EvJ7mvG3sHsAfYn+RW4DHgFoCqOppkP/AA0++0ua2qnljs4JKk2Q0t96r6N2Y+jw6weZZ9dgO7F5BLkrQAXqEqSQ2y3CWpQZa7JDXIcpekBlnuktQgy12SGmS5S1KDLHdJapDlLkkNstwlqUGWuyQ1qM+NwyTW7rx71BHO8eieG0cdQbqseeQuSQ2y3CWpQZa7JDXIcpekBvkHVTXLPwLrqcwjd0lqkOUuSQ0aWu5J7kxyOsn9A2PXJDmU5OFuuXzgsV1Jjic5luSGpQouSZpdnyP3vwK2nDe2EzhcVeuBw902STYA24CN3T53JFm2aGklSb0MLfeq+gzwX+cNbwX2dev7gJsHxu+qqjNV9QhwHNi0OFElSX1d7Dn366rqFEC3XNGNrwJODMyb7MYkSZfQYv9BNTOM1YwTkx1JJpJMTE1NLXIMSXpqu9hyfzzJSoBuebobnwTWDMxbDZyc6Qmqam9VjVfV+NjY2EXGkCTN5GLL/SCwvVvfDhwYGN+W5Ook64D1wJGFRZQkzdfQK1STfAT4ReDaJJPAu4A9wP4ktwKPAbcAVNXRJPuBB4CzwG1V9cQSZZckzWJouVfVG2Z5aPMs83cDuxcSSpK0MF6hKkkNstwlqUGWuyQ1yHKXpAZZ7pLUIMtdkhpkuUtSgyx3SWqQ5S5JDbLcJalBlrskNchyl6QGWe6S1CDLXZIaZLlLUoMsd0lqkOUuSQ0a+p+YJF06a3fePeoI3/fonhtHHUEL4JG7JDXIcpekBi1ZuSfZkuRYkuNJdi7V60iSLrQk5Z5kGfDnwGuADcAbkmxYiteSJF1oqY7cNwHHq+orVfVd4C5g6xK9liTpPKmqxX/S5JeBLVX1G932G4GfrarbB+bsAHZ0m88Dji16kPm5FvjGiDPMl5kvjSst85WWF8x8sX6iqsZmemCp3gqZGcbO+SlSVXuBvUv0+vOWZKKqxkedYz7MfGlcaZmvtLxg5qWwVKdlJoE1A9urgZNL9FqSpPMsVbl/DlifZF2SZwDbgINL9FqSpPMsyWmZqjqb5HbgH4FlwJ1VdXQpXmsRXTaniObBzJfGlZb5SssLZl50S/IHVUnSaHmFqiQ1yHKXpAZZ7udJ8pbutglHk7xn1HmGSfJ7Sb6W5L7u47WjztRXkt9OUkmuHXWWuSR5d5IvdV/fTyX58VFnGibJe5M81OX+eJJnjzrTMElu6b7vvpfksn2L4ZVyaxXLfUCSlzN9Je0Lqmoj8AcjjtTX+6rq+u7j70cdpo8ka4BXAY+NOksP762qF1TV9cAngN8dcZ4+DgHPr6oXAP8B7Bpxnj7uB34J+Myog8zmSrq1iuV+rjcDe6rqDEBVnR5xnpa9D3g7513cdjmqqm8PbD6LKyPzp6rqbLf5WaavNbmsVdWDVTXqK9WHuWJurWK5n+u5wM8nuSfJvyZ5yagD9XR79+v3nUmWjzrMMEluAr5WVV8cdZa+kuxOcgL4Fa6MI/dBvw78w6hDNGIVcGJge7Ibu+w85f4TU5J/An5shofeyfTXYznwUuAlwP4kz6kRv190SOb3A+9m+mjy3cAfMv3NPFJDMr8DePWlTTS3ufJW1YGqeifwziS7gNuBd13SgDMYlrmb807gLPChS5ltNn0yX+aG3lrlcvGUK/eqeuVsjyV5M/CxrsyPJPke0zcHmrpU+WYyV+ZBSf6C6XPCIzdb5iQ/DawDvpgEpk8XfD7Jpqr6+iWMeI6+X2Pgw8DdXAblPixzku3A64DNoz5AedI8vs6Xqyvm1iqeljnX3wGvAEjyXOAZjP6ub3NKsnJg8/VM/1HqslVVX66qFVW1tqrWMv3N8uJRFvswSdYPbN4EPDSqLH0l2QL8DnBTVX1n1HkacsXcWuUpd+Q+xJ3AnUnuB74LbL9cjnjm8J4k1zP9q+GjwJtGmqZNe5I8D/ge8FXgN0ecp48/A64GDnW/IX22qi7r3EleD/wpMAbcneS+qrphxLHOcSXdWsXbD0hSgzwtI0kNstwlqUGWuyQ1yHKXpAZZ7pLUIMtdkhpkuUtSg/4Pu7SV/jId0hcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "energies = defaultdict(int)\n",
    "for k, v in counts.items():\n",
    "    energies[maxcut_obj(k,G)] += v\n",
    "\n",
    "x,y = zip(*energies.items())\n",
    "plt.bar(x,y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pretty good! Most of the results we get correspond to good solutions!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Best string: 00011 with cut: 6\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAb4AAAEuCAYAAADx63eqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAZSElEQVR4nO3df5BVdcHH8c/usrB3ARVBCbDNiVHDAmdsGMccmdLELJuaQXA0HRNZ/FEO/ki0TMbUZh6VJyv8w3FrtFITaFL6A2LVfjBahI8iUFo9j05jjsWqCBvyY1n2+/xxvLK77N29957vOef74/2aYUrhnj3XmXr7ufeeexqMMUYAAESisegTAAAgT4QPABAVwgcAiArhAwBEhfABAKJC+AAAUSF8AICoED4AQFQIHwAgKqMy/wn790t//av07rtSc7M0ebI0fbrU0JD5jwYAuG3HDum116TubmncOOkjH0kykaXswvePf0j33y89+GDy143vj8sDB6SpU6Wbb5Yuvlhqbc3sFAAA7jFG+uMfpeXLpbVrpZaWQ7+3b590xhnS0qXSOeccSodNDda/q7O3V7r6aumRR6S+PqmnZ+g/N25c8uwfeUT68petngIAwE3bt0vnnSf9/e/Snj1JBoYybpw0caL01FPSCSfYPQe74evtlb7wBenZZ5NnVI1SSfrhD6VFi6ydBgDAPW+8Ic2eLb39dpKLkTQ0SOPHJ0mZOdPeedgN38KF0sqV1UevrFSSnnxSmjvX2qkAANzx3nvSKack74IdPFjbYydOlLZtk6ZMsXMu9l49/dvfpJ//vPboSdLevcnLo9whCQCC9NBD0r/+VXv0pOSDL3ffbe9c7IXvBz+o7xmVbd+evNsJAAiKMdK999a3i6TkM5E//nGykWywE749e6Sf/CQ5uzTHuPdeK6cDAHDHhg3SO++kP86qVemPIdkK3+bN0qiUV0YYI/3mN1ZOBwDgjmeeSd7jS2P3bulXv7JzPnbCt2OHnQvS0/6TAQA459//tnOct96ycxw7F7A3NVk5TM/BgxrDN7oAQGDul/S11EdJ+8LiB8excpRjjkkuVk9p9MSJMm+/beGEAACu+O53pdtvN+rtTTdspk61cz52Xuo89VRpzJh0x2huli66yMrpAACK19PTo3Xr1umFF25Vb+++VMcaN0669FI752UnfE1N0vXXJxeipznGkiVWTgcAUIxy7BYuXKgpU6bozjvv1JlnTtIpp6TLzbhxyXd32mDvOr729rovQD8o6f+OOELdxx5r7XQAAPkYKnYzZ87USy+9pD/84Q+67rrrdNZZ/yNpd13HL5WkG26w94XV9sJ3zDHJV23XcbeFxvHj9fCcOZo5c6Y6OzutnRIAIBsjxe7666/Xhz/8YXV1dWn+/Plau7Zdp59+YMCdGKrR3CydeKL09a/bO3e7N3z42tekm26qPn7vfwNpQ2en7lq9Wh0dHWpvb9fixYvV3d1t9dQAAOlUGztJMsZo5cqVmjVrlqZPn66XXnpRzzwzQbNnV5+Ilpbk9q1PP53unbTDmCz89KfGTJpkzPjxxiQvgA781dRkTKlkzKmnGvPKKwMeunPnTrNo0SLT1tZm1q9fn8npAQCqs3//frN27Vpz+eWXm6OPPtqcfvrp5nvf+555/fXXKz5m+/btZt68eeZjH/uY2bhx44Df6+kx5tprkwSMHTt0IkolY1pajLnwQmP+8x/7zymb8BljTG+vMWvWGPOpTxkzapQxjY3GNDQkz3ThQmO2bRv24evXrzdtbW2mvb3d7Nq1K7PTBAAMVE/sjDGmr6/PPP7442by5Mlm6dKlZu/evRX/7K5dxqxYYczxxyd5aGpK/nPyZGPuusuY7dttP6tD7N+ItpJ9+5KrD2u4AnHXrl36xje+oc7OTnV0dGguty0CgEz09PTomWee0erVq7VmzRqddNJJmj9/vi644IIPXr4cTldXl6655hr95S9/0cMPP6zTTjut6p9tTPIF1C0t2dxxfbD8wpdCZ2en2tvbde6552r58uU64ogjij4lAPBe2thJyXt5q1at0pIlS3TZZZfpO9/5jlpq/QRLzrwIn8T6AwAbbMSuLM3KK5I34Stj/QFAbWzGTvJz5fXnXfgk1h8AjMR27Mp8XXn9eRm+MtYfABySVewk/1def16HT2L9AYhblrErC2Hl9ed9+MpYfwBikUfspLBWXn/BhE9i/QEIV16xKwtt5fUXVPjKWH8AQpB37KRwV15/QYZPYv0B8FMRsSsLeeX1F2z4ylh/AFxXZOykOFZef8GHT2L9AXBP0bEri2Xl9RdF+MpYfwCK5ErspPhWXn9RhU9i/QHIl0uxK4tx5fUXXfjKWH8AsuJi7KS4V15/0YZPYv0BsMfV2JXFvvL6izp8Zaw/APVwPXYSK28ohO99rD8A1fAhdmWsvKERvkFYfwAG8yl2EitvJIRvCKw/AL7FroyVNzLCNwzWHxAXX2MnsfJqQfhGwPoDwuZz7MpYebUhfFVi/QHhCCF2EiuvXoSvBqw/wF+hxK6MlVc/wlcH1h/gh9BiJ7HybCB8dWL9AW4KMXZlrDw7CF9KrD+geCHHTmLl2Ub4LGD9AfkLPXZlrDz7CJ9FrD8gW7HETmLlZYnwWcb6A+yKKXZlrLxsEb6MsP6A+sUYO4mVlxfClyHWH1C9WGNXxsrLD+HLAesPGFrssZNYeUUgfDlh/QEJYncIK68YhC9nrD/EiNgNxMorFuErAOsPMSB2Q2PlFY/wFYj1h9AQu8pYee4gfAVj/cF3xG5krDy3ED5HsP7gE2JXHVaemwifQ1h/cBmxqw0rz12Ez0GsP7iC2NWOlec+wuco1h+KQuzqx8rzA+FzHOsPeSB26bDy/EL4PMD6QxaInR2sPP8QPo+w/pAWsbOHlecvwucZ1h9qRezsY+X5jfB5ivWH4RC7bLDywkD4PMb6Q3/ELlusvHAQvgCw/uJF7LLHygsP4QsE6y8exC4/rLwwEb7AsP7CROzyxcoLG+ELEOsvDMSuGKy88BG+gLH+/EPsisPKiwfhCxzrz33ErnisvLgQvkiw/txC7NzAyosT4YsI669YQ8VuwYIFmjdvHrErACsvXoQvQqy//BA797DyQPgixfrLDrFzFysPEuGLHuvPDmLnNlYe+iN8YP3Vidj5gZWHwQgfPsD6Gxmx8wcrD5UQPgzA+jscsfMPKw/DIXwYUuzrj9j5iZWHahA+VBTb+iN2fmPloVqEDyMKef0RO/+x8lArwoeqhLT+iF04WHmoB+FDTXxdf8QuLKw8pEH4UDNf1h+xCxMrD2kRPtTNxfVH7MLFyoMthA+puLD+iF34WHmwifDBirzXH7GLAysPWSB8sCbr9Ufs4sLKQ1YIH6yzuf6IXXxYecga4UMm0qw/YhcvVh7yQPiQqWrXH7GLGysPeSJ8yFyl9UfsILHykD/Ch9x0dnZq0aJF+vjHP65JkyZp7dq1xC5irDwUZVTRJ4Dw9V92u3fv1osvvqienh6tWLFCl1xySdGnhwL0X3lr1qxh5SFXjUWfAMLU09OjdevWaeHChZoyZYruvPNOzZo1S1u2bNH27du1cuVK3XrrrVq8eLG6u7uLPl3kxBijlStXatasWZo+fbo2b95M9JA7XuqENbW+Z+fCt74gP7yXB1cQPqRi4wMqLn7nJ+zhvTy4hvChZll8GpP1FyZWHlxE+FCVvC49YP2FgZUHlxE+VFTUdXasP7+x8uA6wocBXLqonPXnF1YefEH44FTsBmP9+YGVB58Qvki5HLuhsP7cxMqDjwhfRHyL3WCsP7ew8uArwhc432M3FNZfsVh58B3hC1CIsRuM9VcMVh5CQPgCEUPshsL6ywcrDyEhfB6LNXaDsf6yxcpDaAifZ4hdZaw/u1h5CBXh8wCxqx7rzw5WHkJG+BxF7NJh/dWHlYcYED6HEDu7WH+1YeUhFoSvYMQue6y/4bHyEBvCVwBilz/W39BYeYgR4csJsXMD6y/BykPMCF+GiJ2bYl9/rDzEjvBZRuz8Edv6Y+UBCcJnAbHzVyzrj5UHHEL46kTswhLq+mPlAYcjfDUgdmELbf2x8oChEb4RELv4+L7+WHnA8AjfEIgdfF1/rDxgZITvfcQOQ/Fl/bHygOpFHT5ih2q4vv5YeUBtogsfsUO9XFt/rDygPlGEj9jBFlfWHysPqF+w4SN2yFJR64+VB6SXffjefVfatEnauVNqbpYmT5ZOO00aNcr6jyJ2yFPe64+VB9iRXfheeEFavlx68klpzBipr09qaJCMkUaPlpYskRYvTkKYArFD0bJef6w8wC774du7V5o/X/rtb6V9+5LgDaX8P9zvf1+68sqafgSxg2uyWn+sPMA+u+Hbu1c680zp5ZeT/16N1lbp1lulb31r2D9G7OADW+uPlQdkx274vvQl6amnqo9eWWur9NBD0oIFA/42sYOP0q4/Vh6QLXvh27JFOv302qNXNnWq9MYb6jlwgNghCLWuP1YekA974bvsMunRR6WDB+t6eG+ppPvmzNF/Pf88sUMwql1/rDwgP3bC192dfDpz3766D9En6Z8nnqjGp58mdghOpfXHygPy12jlKFu2JJcspDyRj7z5JtFDkObOnautW7fKGKOZM2eqs7NTXV1dmj9/vm6//XatWbNGd999N9EDcmAnfO++m1yfl1a97w8CHjjyyCPV0dGhBx98UF/5yld0/PHH67jjjtPmzZt5aRPIkZ3wjR6dXJyeVlNT+mMADuvq6lJHR4cmTJigc845R0888YQ2bNhQ9GkBUbETvg99qO4PtfRnJkywcDKAe4wxWrlypWbNmqXp06dr69atWrNmjTo6OtTe3q7Fixeru7u76NMEomAnfKecIh15ZKpD7JfUceCAli1bpm3btinQ785GhIZ7L2+o9/4AZMtO+BoapJtuSi5Er9PoMWM0+6GHtGfPHp1//vmaMWMGEYTXBq+8Su/lld/7Y/0B+bB3Hd/OnclF6PV8QKWpSTrrLOn9f9s1xmjTpk1avXq1Vq9erVKppAULFmj+/Pn6xCc+oQYb7ycCGar3ujxX7vcHhMzuV5atWiVdfrm0Z09tj5s4MbkkYtq0w36LCMIntq7Lc+1u70BI7N+d4YEHpBtuqG75NTVJRx0lbdggnXzyiH+cCMJltr99hfUHZCOb+/GtXy9de6305ptJAAffmqilJbnu7+yzpQcfHHLpjYQIwhVZf/sK6w+wK7sb0RqT3Hl9+XLpd7+Tdu9O7ro+YYL01a9KV12VvCdo5UcRQRQjr+/YZP0B9mQXvoIQQeShqO/YZP0B6QUXvv6IILJQ9J0UWH9AOkGHrz8iiLRcu5MC6w+oTzTh648IolZFr7xKWH9A7aIMX39EEMNxbeVVwvoDqhd9+PojgujP1ZVXCesPqA7hq4AIxsuXlVcJ6w8YHuGrAhGMh28rrxLWH1AZ4asREQyT7yuvEtYfcDjClwIRDEMoK68S1h8wEOGzhAj6J9SVVwnrD0gQvgwQQfeFvvIqYf0BhC9zRNAtsa28Slh/iBnhyxERLFasK68S1h9iRfgKQgTzw8obHusPsSF8DiCC2WHlVYf1h5gQPscQQTtYefVh/SEGhM9hRLA+rLx0WH8IHeHzBBEcGSvPLtYfQkX4PEQED8fKywbrDyEifJ6LPYKsvHyw/hASwheQ2CLIyssX6w+hIHyBCjmCrLxisf7gO8IXgZAiyMpzA+sPPiN8kfE1gqw8N7H+4CPCFzFfIsjKcxvrD74hfJDkZgRZeX5h/cEXhA+HcSGCrDw/sf7gA8KHYeUdQVZeGFh/cBnhQ9WyjiArLyysP7iK8KEuNiPIygsb6w+uIXxILU0EWXlxYP3BJYQPVlUbQVZenFh/cAHhQ2YqRfCss87SihUr9PLLL7PyIsT6Q9EIH3JhjNGf/vQn3XXXXVq3bp2OOuooXXnllbrooouculge+WH9oSiNRZ8A4vDWW29p+fLlevXVV/Xcc89p7dq16unp0fnnn68ZM2Zo2bJl2rZtm/j3sHjMnTtXW7dulTFGM2fOVGdnZ9GnhEiw+JCpkd7Lc+FieRSP9Yc8ET5kptZPbBLBuPHeH/JC+GCdjU9sEsF4sf6QNcIHq7K4Lo8Ixof1hywRPliR13V5RDAurD9kgfAhtaK+fYUIxoH1B9sIH+rm0revEMHwsf5gC+FDXVz+jk0iGC7WH2wgfKiJSyuvGkQwTKw/pEH4UDWXV141iGBYWH+oF+HDiHxbedUgguFg/aFWhA/D8n3lVYMI+o/1h1oQPgwpxJVXDSLoN9YfqkH4cJgYVl41iKCfWH8YCeHDB2JdedUggv5h/aESwgdJrLxaEEF/sP4wFMIXOVZeOkTQD6w/9Ef4IsbKs4sIuo31hzLCFyFWXvaIoLtYfyB8kWHl5Y8Iuof1FzfCFwlWnhuIoFtYf3EifBFg5bmJCLqB9RcfwhcwVp4/iGDxWH/xIHyBYuX5iwgWh/UXB8IXGFZeWIhgMVh/YSN8AWHlhY0I5ov1Fy7CFwBWXnyIYH5Yf+EhfJ5j5YEIZo/1FxbC5ylWHoZCBLPF+gsD4fMQKw/VIILZYP35j/B5hJWHehFB+1h//iJ8nmDlwRYiaA/rz0+Ez3GsPGSJCNrB+vML4XMYKw95IoLpsP78QfgcxMpD0Yhg/Vh/7iN8jmHlwTVEsHasP7cRPkew8uADIlgb1p+bCJ8DWHnwERGsDuvPPYSvQKw8hIIIjoz15w7CVxBWHkJFBCtj/bmB8OWMlYeYEMGhsf6KRfhyxMpDzIjgQKy/4hC+HLDygIGI4CGsv/wRvoyx8oDhEUHWX94IX0ZYeUDtYo8g6y8fhC8DrDwgvVgjyPrLHuGziJUHZCPGCLL+skP4LGHlAfmIKYKsv2wQvpRYeUBxYokg688uwpcCKw9wR+gRZP3ZQ/jqwMoD3BZyBFl/6RG+GrHyAL+EGEHWXzqEr0qsPMB/oUWQ9VcfwlcFVh4QnlAiyPqrHeEbBisPiEMIEWT9VY/wVcDKA+LkcwRZf9UhfIOw8gCU+RpB1t/wCF8/rDwAlfgWQdZfZYRPrDwAtfEpgqy/w0UfPlYegDR8iCDrb6Bow8fKA2Cb6xFk/SWiDB8rD0DWXI0g6y+y8LHyABTBxQjGvP6iCR8rD4ALXIpgrOsv+PCx8gC4ypUIxrb+gg4fKw+AL4qOYEzrL8jwsfIA+KzICMaw/oILHysPQEiKiGDo6y+Y8LHyAIQu7wiGuv6CCB8rD0Bs8opgiOvP6/Cx8gAgnwiGtP68DR8rDwAOl2UEQ1l/3oWPlQcA1ckqgr6vP6/Cx8oDgPrYjqDP68+L8LHyAMAemxH0cf05Hz5WHgBkx0YEfVt/2YfvtdekZ5+Vdu6UmpulY4+Vzj1XGjdu2Iex8gAgX2kjWOv6M0bauFF65RWpu1saO1b66Eelz3xGamy0+cwGyiZ8Bw9Kv/61dM890qZN0qhR0oEDyTNpbpZ6e6VLLpGuu06aMeOwh7PyAKBY9UawmvW3e7f0s58liXj77SSAvb1SU1Pyq1SSbrhBuuIKadKkbJ6cXTt2GDN7tjHjxhmTPJ+hf40aZUypZMyyZcb09RljjOnr6zOPP/64mTx5slm6dKnZu3ev9dMDANSmr6/PbNy40dx4442mra3NnHTSSea2224zW7duNX3v///3YOvXrzdtbW2mvb3d7Nq164O/v3mzMRMnGjN27PCJKJWSjDz1lP3nY3fx7dolffKT0j//KfX0VPeY1lbpiivU9e1vs/IAwHGmhiU4eP1NnDhXn/50sviqVSpJv/iF9PnP23sO9sJnjDRnjvT889L+/TU9tHf0aN0yZoyarr6a9/IAwBPVRrCzs1MLF35Tb721QT09Y2v+Oa2tSVpOPtnOedsL38aN0mc/K733Xl0PP3DkkWp+553kBV4AgFdGiuAtt+zXf/93ow4eHF3zsRsbpQsvlB57zM652gvfBRdIv/xlsvzqMX689Oij0he/aOV0AADFGBzBlpbxev31Tdq3r7XuY7a0SG++KU2YkP787ITvnXekadNqfonzMGeckVz6AAAIgjFG99zzv7rttjYdOFD/21itrdIdd0g33pj+nOxcKfHnPyc5TmvLlvTHAAA4o6GhQTt3npgqepK0Z4/0+9/bOSc74du5s/6XOPvbsyf9MQAATunqsnOcHTvsHGeUlaOUSpKF+z3t6+tTyfIdhAEARfu+pCWpj9Ja/1uEA9gJ37RpyWX3KbVMmybzxhsWTggA4Ir77pO++c10HwNpbEy+zswGOy91nnyyNGVKumOUStJVV1k5HQCAOy68MP0xSiWpvT39cSRb4WtokJYuHfGLp4dljLR4sZXTAQC4Y+pU6eyz070j1taWfDGYDfa+//rii+t/7OjR0nnnJXduAAAE5+abk9VWj7FjpVtusXcu9sI3dmxyAXutz6ypKQnej35k7VQAAG6ZM0e65praP6BSKkmf+5x06aX2zsXuHY/OOSe510S18Rs9WjruOOm556Sjj7Z6KgAAt9xzj3T55dXHr7U1ycpjj1m5cOAD9m/1N29ecpXhmWcmF7WPHuJ72caOTeJ42WXS5s3Ji7cAgKA1NEj33y898IB0/PHJx0KGCtr48ckLgXfcIT3xxNAZSXUeVm9LNNirr0orVkjr1ye3LGpulo45Rlq0KLkRbZoPwwAAvGVM8g2V990nbduW3KqoVJJOOEFaskQ699zs7lmQbfgAAHCM/Zc6AQBwGOEDAESF8AEAokL4AABRIXwAgKgQPgBAVAgfACAqhA8AEBXCBwCICuEDAETl/wE7OoL2ci00IgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# get the best solution:\n",
    "best_cut, best_solution = min([(maxcut_obj(x,G),x) for x in counts.keys()], key=itemgetter(0))\n",
    "print(f\"Best string: {best_solution} with cut: {-best_cut}\")\n",
    "\n",
    "# color the graph nodes by part\n",
    "colors = ['r' if best_solution[node] == '0' else 'b' for node in G]\n",
    "nx.draw(G,node_color = colors, pos=nx.bipartite_layout(G, [0,1,2]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Using statevector simulator\n",
    "\n",
    "What if you want to get the full QAOA state? Luckily, Qiskit provides `statevector_simulator` backend. \n",
    "\n",
    "### What is a state vector?\n",
    "\n",
    "State vector is a vector of amplitudes, e.g.: $\\frac{1}{\\sqrt{2}}\\bigg(|00\\rangle + |11\\rangle\\bigg) \\rightarrow (\\frac{1}{\\sqrt{2}}, 0, 0, \\frac{1}{\\sqrt{2}})$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.70710678+0.j, 0.        +0.j, 0.        +0.j, 0.70710678+0.j])"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "qc = QuantumCircuit(2)\n",
    "qc.h(0)\n",
    "qc.cx(0,1)\n",
    "sv = execute(qc, backend=Aer.get_backend('statevector_simulator')).result().get_statevector()\n",
    "sv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Less luckily, qubit ordering considerations still apply. Below are some helper functions we're going to use."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Endianness conversion tools from https://github.com/Qiskit/qiskit-terra/issues/1148#issuecomment-438574708\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "def state_num2str(basis_state_as_num, nqubits):\n",
    "    return '{0:b}'.format(basis_state_as_num).zfill(nqubits)\n",
    "\n",
    "def state_str2num(basis_state_as_str):\n",
    "    return int(basis_state_as_str, 2)\n",
    "\n",
    "def state_reverse(basis_state_as_num, nqubits):\n",
    "    basis_state_as_str = state_num2str(basis_state_as_num, nqubits)\n",
    "    new_str = basis_state_as_str[::-1]\n",
    "    return state_str2num(new_str)\n",
    "\n",
    "def get_adjusted_state(state):\n",
    "    nqubits = np.log2(state.shape[0])\n",
    "    if nqubits % 1:\n",
    "        raise ValueError(\"Input vector is not a valid statevector for qubits.\")\n",
    "    nqubits = int(nqubits)\n",
    "\n",
    "    adjusted_state = np.zeros(2**nqubits, dtype=complex)\n",
    "    for basis_state in range(2**nqubits):\n",
    "         adjusted_state[state_reverse(basis_state, nqubits)] = state[basis_state]\n",
    "    return adjusted_state"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Running using full-state (`statevector`) simulator\n",
    "\n",
    "Note that a few things are different when using the full statevector simulator. First, we must not measure. Second, as the result, we get amplitudes instead of samples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_qaoa_circuit_sv(G, beta, gamma):\n",
    "    assert(len(beta) == len(gamma))\n",
    "    p = len(beta) # infering number of QAOA steps from the parameters passed\n",
    "    N = G.number_of_nodes()\n",
    "    qc = QuantumCircuit(N,N)\n",
    "    # first, apply a layer of Hadamards\n",
    "    qc.h(range(N))\n",
    "    # second, apply p alternating operators\n",
    "    for i in range(p):\n",
    "        qc += get_cost_operator_circuit(G,gamma[i])\n",
    "        qc += get_mixer_operator_circuit(G,beta[i])\n",
    "    # no measurement in the end!\n",
    "    return qc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Second, we need to get which amplitudes correspond to which computational basis states (= binary strings)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "def state_to_ampl_counts(vec, eps=1e-15):\n",
    "    \"\"\"Converts a statevector to a dictionary\n",
    "    of bitstrings and corresponding amplitudes\n",
    "    \"\"\"\n",
    "    qubit_dims = np.log2(vec.shape[0])\n",
    "    if qubit_dims % 1:\n",
    "        raise ValueError(\"Input vector is not a valid statevector for qubits.\")\n",
    "    qubit_dims = int(qubit_dims)\n",
    "    counts = {}\n",
    "    str_format = '0{}b'.format(qubit_dims)\n",
    "    for kk in range(vec.shape[0]):\n",
    "        val = vec[kk]\n",
    "        if val.real**2+val.imag**2 > eps:\n",
    "            counts[format(kk, str_format)] = val\n",
    "    return counts\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Let's try plugging in the optimal parameters we found using the QASM simulator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'00000': (0.013008244489519906+0.013030397554903457j),\n",
       " '00001': (0.14334369756512017-0.0797375314273633j),\n",
       " '00010': (0.14334369756511992-0.07973753142736324j),\n",
       " '00011': (0.046149503038565465-0.6020701181608085j),\n",
       " '00100': (-0.05384089700023571+0.0005250354845161692j),\n",
       " '00101': (-0.03502587128054775-0.014045081674946976j),\n",
       " '00110': (-0.03502587128054772-0.014045081674946893j),\n",
       " '00111': (-0.07540147572810121+0.1250784478518241j),\n",
       " '01000': (-0.05384089700023556+0.0005250354845162039j),\n",
       " '01001': (-0.03502587128054775-0.014045081674946997j),\n",
       " '01010': (-0.035025871280547606-0.014045081674946927j),\n",
       " '01011': (-0.07540147572810113+0.12507844785182418j),\n",
       " '01100': (-0.07540147572810117+0.1250784478518242j),\n",
       " '01101': (-0.03502587128054777-0.014045081674946863j),\n",
       " '01110': (-0.03502587128054768-0.014045081674946933j),\n",
       " '01111': (-0.05384089700023567+0.0005250354845161339j),\n",
       " '10000': (-0.05384089700023568+0.0005250354845163011j),\n",
       " '10001': (-0.03502587128054772-0.014045081674946967j),\n",
       " '10010': (-0.03502587128054764-0.014045081674946959j),\n",
       " '10011': (-0.07540147572810112+0.12507844785182412j),\n",
       " '10100': (-0.07540147572810121+0.12507844785182415j),\n",
       " '10101': (-0.035025871280547634-0.014045081674946922j),\n",
       " '10110': (-0.03502587128054772-0.014045081674946901j),\n",
       " '10111': (-0.05384089700023566+0.0005250354845162397j),\n",
       " '11000': (-0.07540147572810112+0.12507844785182418j),\n",
       " '11001': (-0.03502587128054766-0.014045081674946919j),\n",
       " '11010': (-0.035025871280547675-0.014045081674946868j),\n",
       " '11011': (-0.053840897000235706+0.0005250354845161703j),\n",
       " '11100': (0.046149503038565354-0.6020701181608085j),\n",
       " '11101': (0.14334369756512017-0.07973753142736309j),\n",
       " '11110': (0.14334369756512005-0.07973753142736321j),\n",
       " '11111': (0.013008244489519937+0.013030397554903476j)}"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "backend = Aer.get_backend('statevector_simulator')\n",
    "qc = get_qaoa_circuit_sv(G, optimal_theta[:p], optimal_theta[p:])\n",
    "sv = get_adjusted_state(execute(qc, backend).result().get_statevector())\n",
    "state_to_ampl_counts(sv)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Let's verify that we get the same energy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-5.296260485608114"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def compute_maxcut_energy_sv(sv, G):\n",
    "    \"\"\"Compute objective from statevector\n",
    "    For large number of qubits, this is slow. \n",
    "    \"\"\"\n",
    "    counts = state_to_ampl_counts(sv)\n",
    "    return sum(maxcut_obj(np.array([int(x) for x in k]), G) * (np.abs(v)**2) for k, v in counts.items())\n",
    "\n",
    "compute_maxcut_energy_sv(sv, G)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Yay! Energy matches"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Now, we can run the full optimization loop with the statevector simulator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_black_box_objective_sv(G,p):\n",
    "    backend = Aer.get_backend('statevector_simulator')\n",
    "    def f(theta):\n",
    "        # let's assume first half is betas, second half is gammas\n",
    "        beta = theta[:p]\n",
    "        gamma = theta[p:]\n",
    "        qc = get_qaoa_circuit_sv(G,beta, gamma)\n",
    "        sv = execute(qc, backend).result().get_statevector()\n",
    "        # return the energy\n",
    "        return compute_maxcut_energy_sv(get_adjusted_state(sv), G)\n",
    "    return f"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Copying and pasting the code for the `qasm_simulator`..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "     fun: -5.946219918156263\n",
       "   maxcv: 0.0\n",
       " message: 'Maximum number of function evaluations has been exceeded.'\n",
       "    nfev: 2500\n",
       "  status: 2\n",
       " success: False\n",
       "       x: array([2.18869004, 2.0161459 , 1.98647176, 1.8248621 , 0.54904817,\n",
       "       2.26821077, 6.08799377, 2.38220732, 3.75526057, 6.03938487])"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# p is the number of QAOA alternating operators\n",
    "p = 5\n",
    "obj = get_black_box_objective_sv(G, p)\n",
    "\n",
    "init_point = np.array([0.81069872, 2.2067517 , 0.83830696, 2.15579759, 0.37060699,\n",
    "       2.42068091, 6.1575306 , 2.2453419 , 3.85060091, 6.137845  ])\n",
    "# We are going to limit the number of iterations to 2500\n",
    "res_sv = minimize(obj, init_point, method='COBYLA', options={'maxiter':2500, 'disp': True})\n",
    "res_sv"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "We get **much better solution** if we have access to the full state vector! In fact, if we run it long enough, we will get arbitratily close to the true optimum (to verify it yourself, try setting `'maxiter':10000` above)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Best energy with samples: -5.30, with full state vector: -5.95\n"
     ]
    }
   ],
   "source": [
    "print(f\"Best energy with samples: {res_sample['fun']:.2f}, with full state vector: {res_sv['fun']:.2f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "This is because with samples we get stuck in a lower quality point:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Angles found with samples: beta=[1.83749394 2.1345936  1.9319531  2.09862998 0.347631  ], gamma=[2.04618358 6.11249139 2.31785137 3.56479392 5.87260613]\n",
      "Angles found with full state: beta=[2.18869004 2.0161459  1.98647176 1.8248621  0.54904817], gamma=[2.26821077 6.08799377 2.38220732 3.75526057 6.03938487]\n"
     ]
    }
   ],
   "source": [
    "print(f\"Angles found with samples: beta={res_sample['x'][:p]}, gamma={res_sample['x'][p:]}\")\n",
    "print(f\"Angles found with full state: beta={res_sv['x'][:p]}, gamma={res_sv['x'][p:]}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that these are two different points!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "An important point to make is that this is **not a local minimum** that COBYLA gets stuck in when using the sampling. We can easily verify that by running COBYLA with full state from that point"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "     fun: -5.379488009620951\n",
       "   maxcv: 0.0\n",
       " message: 'Maximum number of function evaluations has been exceeded.'\n",
       "    nfev: 50\n",
       "  status: 2\n",
       " success: False\n",
       "       x: array([1.84912811, 2.1430627 , 1.94267686, 2.08362536, 0.36741251,\n",
       "       2.04556667, 6.12254338, 2.30958241, 3.57688979, 5.86524965])"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# p is the number of QAOA alternating operators\n",
    "p = 5\n",
    "obj = get_black_box_objective_sv(G, p)\n",
    "\n",
    "init_point = res_sample['x']\n",
    "# 'rhobeg' controls initial changes to the variables. \n",
    "# We set it low to make sure we are searching a small neighborhood of the optimal point\n",
    "res_sv2 = minimize(obj, init_point, method='COBYLA', options={'maxiter':50, 'disp': True, 'rhobeg': 0.001})\n",
    "res_sv2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The lesson here is that our simple deterministic optimizers (like COBYLA) assume that *every function value they get is true*. This causes them to thing that a point is a local minimum, when in fact it is not. We have to bake the assumption of stochaisticity into the optimizers we use!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### That was pretty involved! Is there an easier way?\n",
    "\n",
    "So far we have effectively reimplemented the entire stack of tools used in Qiskit. Wasn't too bad, huh? Note that you can fairly easily do the above in any framework you prefer (Google cirq, Rigetti PyQuil, ...), all you need to learn is how to apply gates. What follows is more specific to IBM Qiskit (though not entirely).\n",
    "\n",
    "Luckily, you don't have to reimplement QAOA every time you want to use it (though I do recommend doing it at least once to understand what's going on!)\n",
    "\n",
    "Qiskit provides a number of tools for optimization, some as simple to use as plugging your graph into a declarative interface. Today we're going to look at some of them. I have copied and pasted some of the code from Qiskit here for educational purposes."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Recap\n",
    "\n",
    "To solve MAXCUT we had to do the following:\n",
    "\n",
    "1. Choose a problem class (maxcut) and a problem instance (graph)\n",
    "2. Come up with a general Hamiltonian for the problem class\n",
    "3. Build specific Hamiltonian for the problem instance\n",
    "4. From the Hamiltonian, construct QAOA circuit\n",
    "5. Wrap QAOA circuit in a black box function\n",
    "6. Optimize QAOA parameters\n",
    "7. Get solution by running QAOA circuit with optimal parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Qiskit to help!\n",
    "\n",
    "Let's see how Qiskit can help us with steps 3 and 4:\n",
    "\n",
    "\n",
    "1. Choose a problem class (maxcut) and a problem instance (graph)\n",
    "2. Come up with a general Hamiltonian for the problem class\n",
    "3. **Build specific Hamiltonian for the problem instance**\n",
    "4. **From the Hamiltonian, construct QAOA circuit**\n",
    "5. Wrap QAOA circuit in a black box function\n",
    "6. Optimize QAOA parameters\n",
    "7. Get solution by running QAOA circuit with optimal parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "### Building the circuit in Qiskit\n",
    "\n",
    "First thing we're going to try is to use Qiskit tools to build circuit instead of building it by hand. To do that, we will have to construct a Hamiltonian object and pass it to QAOA class initializer."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Building the Hamiltonian in Qiskit\n",
    "\n",
    "Qiskit uses weird notation that we will not get into"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "# copied from qiskit.optimization.ising.max_cut\n",
    "\n",
    "def get_cost_operator(weight_matrix):\n",
    "    \"\"\"Generate Hamiltonian for the max-cut problem of a graph.\n",
    "    Args:\n",
    "        weight_matrix (numpy.ndarray) : adjacency matrix.\n",
    "    Returns:\n",
    "        WeightedPauliOperator: operator for the Hamiltonian\n",
    "        float: a constant shift for the obj function.\n",
    "    \"\"\"\n",
    "    num_nodes = weight_matrix.shape[0]\n",
    "    pauli_list = []\n",
    "    shift = 0\n",
    "    for i in range(num_nodes):\n",
    "        for j in range(i):\n",
    "            if weight_matrix[i, j] != 0:\n",
    "                x_p = np.zeros(num_nodes, dtype=np.bool)\n",
    "                z_p = np.zeros(num_nodes, dtype=np.bool)\n",
    "                z_p[i] = True\n",
    "                z_p[j] = True\n",
    "                pauli_list.append([0.5 * weight_matrix[i, j], Pauli(z_p, x_p)])\n",
    "                shift -= 0.5 * weight_matrix[i, j]\n",
    "    return WeightedPauliOperator(paulis=pauli_list), shift"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### First, let's construct a Hamiltonian \n",
    "\n",
    "Recall the MAXCUT objective is $\\max_{\\mathbb{s}}\\frac{1}{2}\\sum_{ij\\in E}(1-s_is_j) = \\max_{\\mathbb{s}}\\frac{1}{2}\\sum_{ij\\in E}(-s_is_j) + const$. Qiskit conveniently takes care of the constant in form of offset.\n",
    "\n",
    "Corresponding Hamiltonian is $\\frac{1}{2}\\sum_{ij\\in E}(-Z_iZ_j)$, but Qiskit by convention always searches for the _ground_ state, so the sign is flipped."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "IZIIZ\t(0.5+0j)\n",
      "IZIZI\t(0.5+0j)\n",
      "IZZII\t(0.5+0j)\n",
      "ZIIIZ\t(0.5+0j)\n",
      "ZIIZI\t(0.5+0j)\n",
      "ZIZII\t(0.5+0j)\n",
      "\n"
     ]
    }
   ],
   "source": [
    "w = nx.adjacency_matrix(G, nodelist=range(G.number_of_nodes())) # have to pass nodelist so the columns are ordered properly\n",
    "C, offset = get_cost_operator(w)\n",
    "print(C.print_details())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looks correct!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Verifying that Hamiltonian is correct\n",
    "\n",
    "Recall that since \"Hamiltonian\" is just a fancy word for matrix*, we can find its ground state (lowest eigenvalue eigenvector) classically with any eigensolver. We can use Qiskit `NumPyEigensolver`, which in is just a wrapper around `np.linalg.eig` and `scisparse.linalg.eigs`\n",
    "\n",
    "*not really"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Let's run NumPyEigensolver on our Hamiltonian (Hermitian matrix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'eigenvalues': array([-3.+0.j]), 'eigenstates': ListOp([VectorStateFn(Statevector([0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j,\n",
       "             0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j,\n",
       "             0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j,\n",
       "             0.+0.j, 0.+0.j, 0.+0.j, 1.+0.j, 0.+0.j, 0.+0.j, 0.+0.j,\n",
       "             0.+0.j, 0.+0.j, 0.+0.j, 0.+0.j],\n",
       "            dims=(2, 2, 2, 2, 2)), coeff=1.0, is_measurement=False)], coeff=1.0, abelian=False)}"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ee = NumPyEigensolver(C, k=1) # k=1 means lowest eigenvalue and eigenvector\n",
    "\n",
    "result = ee.run()\n",
    "result"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Now, we need to convert the eigenvectors of $2^n\\times 2^n$ matrix into solutions to our MAXCUT problem (vectors of size $n$)\n",
    "\n",
    "Luckily, Qiskit has a function for that!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Solution: [0 0 0 1 1] with cut: [-6.+0.j]\n"
     ]
    }
   ],
   "source": [
    "from qiskit.optimization.applications.ising.common import sample_most_likely\n",
    "print(f\"Solution: {sample_most_likely(result['eigenstates'][0])} with cut: {result['eigenvalues']+offset}\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And it's a correct solution. So we've built our Hamiltonian correctly. Note that you can use the same process (building a Hamiltonian, running ExactEigensolver on it) in prototyping your own algorithms."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Q: Can I build my Hamiltonian as a matrix?\n",
    "\n",
    "Yes you can! Let's try doing it it for the same five-node graph:\n",
    "\n",
    "$C=\\frac{1}{2}\\sum_{ij\\in E}Z_iZ_i = \\frac{1}{2}(Z_0Z_3 + Z_0Z_4 + Z_1Z_3 + Z_1Z_4 + Z_2Z_3 + Z_2Z_4) = \\frac{1}{2}(Z\\otimes I \\otimes I \\otimes Z \\otimes I + Z\\otimes I \\otimes I \\otimes I \\otimes Z + I\\otimes Z \\otimes I \\otimes Z \\otimes I + I\\otimes Z \\otimes I \\otimes I \\otimes Z + I\\otimes I \\otimes Z \\otimes Z \\otimes I + I\\otimes I \\otimes Z \\otimes I \\otimes Z)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "I = np.eye(2,2)\n",
    "Z = np.array([[1.,0.], [0.,-1.]])\n",
    "\n",
    "Z0Z3 = np.kron(np.kron(np.kron(np.kron(Z,I), I), Z), I)\n",
    "Z0Z4 = np.kron(np.kron(np.kron(np.kron(Z,I), I), I), Z)\n",
    "Z1Z3 = np.kron(np.kron(np.kron(np.kron(I,Z), I), Z), I)\n",
    "Z1Z4 = np.kron(np.kron(np.kron(np.kron(I,Z), I), I), Z)\n",
    "Z2Z3 = np.kron(np.kron(np.kron(np.kron(I,I), Z), Z), I)\n",
    "Z2Z4 = np.kron(np.kron(np.kron(np.kron(I,I), Z), I), Z)\n",
    "\n",
    "C_numpy = 0.5 * (Z0Z3 + Z0Z4 + Z1Z3 + Z1Z4 + Z2Z3 + Z2Z4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "PS If you know a better way of doing this -- please let me know!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Let's check that the Hamiltonian we've built is correct\n",
    "\n",
    "No shenanigans here -- just numpy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Solution: [1. 1. 0. 0. 0.] with cut: -6.0\n"
     ]
    }
   ],
   "source": [
    "# Get eigenvalues and eigenvectors\n",
    "eigval, eigvec = np.linalg.eig(C_numpy)\n",
    "# Sort eigenvalues\n",
    "idx = eigval.argsort()\n",
    "# Get lowest\n",
    "lowest_eigval = eigval[idx][0]\n",
    "lowest_eigvec = eigvec[idx][0]\n",
    "# Only need qiskit to convert eigenvector into the solution to our problem\n",
    "print(f\"Solution: {sample_most_likely(lowest_eigvec)} with cut: {lowest_eigval+offset}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And this is correct!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Constructing QAOA circuit for our Hamiltonian\n",
    "\n",
    "Qiskit provides a class that constructs QAOA circuit for us from a Hamiltonian. \n",
    "\n",
    "The abstraction Qiskit uses for this is \"Variational Form\". Roughly speaking, variational form is just a parameterized circuit that prepares some parameterized state $\\psi(\\theta)$. For QAOA, as you might recall, this parameterized state is:\n",
    "\n",
    "$$  \\begin{aligned} \n",
    "    |\\psi{(\\mathbb{\\theta})}> & = |\\psi{(\\mathbb{\\beta},\\mathbb{\\gamma})}>\\\\\n",
    "    &=  e^{-i\\beta_p B}e^{-i\\gamma_p C}\\cdots e^{-i\\beta_1 B}e^{-i\\gamma_1 C} H^{\\otimes n} |0>.\n",
    "  \\end{aligned} $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiskit.aqua.algorithms.minimum_eigen_solvers.qaoa.var_form import QAOAVarForm\n",
    "\n",
    "p = 5\n",
    "var_form = QAOAVarForm(C.to_opflow(), p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Constructing QAOA circuit for our Hamiltonian (cont'd)\n",
    "\n",
    "To construct circuit, we need to pass the parameters $\\mathbb{\\beta}, \\mathbb{\\gamma}$ to the variational form."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'00000': 118, '00001': 18, '10000': 109, '10001': 6, '10010': 8, '10011': 5, '10100': 6, '10101': 12, '10110': 9, '10111': 84, '11000': 2, '11001': 31, '11010': 33, '11011': 28, '11100': 16, '11101': 26, '11110': 31, '11111': 130, '00010': 21, '00011': 25, '00100': 18, '00101': 34, '00110': 31, '01000': 92, '01001': 3, '01010': 2, '01011': 8, '01100': 2, '01101': 6, '01110': 6, '01111': 104}\n"
     ]
    }
   ],
   "source": [
    "N = G.number_of_nodes()\n",
    "\n",
    "# Lower and upper bounds: beta \\in {0, pi}, gamma \\in {0, 2*pi}\n",
    "lb = np.zeros(2*p)\n",
    "ub = np.hstack([np.full(p, np.pi), np.full(p, 2*np.pi)])\n",
    "random_point = np.random.uniform(lb, ub, 2*p)\n",
    "\n",
    "# Build the circuit\n",
    "qc = var_form.construct_circuit(random_point)\n",
    "\n",
    "# Don't forget to measure! This requires a little hack\n",
    "if not qc.cregs:\n",
    "    c = ClassicalRegister(N, name='c')\n",
    "    qc.add_register(c)\n",
    "qc.measure(qc.qregs[0], qc.cregs[0])\n",
    "# And run\n",
    "backend = Aer.get_backend('qasm_simulator')\n",
    "counts = execute(qc, backend).result().get_counts()\n",
    "print(counts)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Let's compare Qiskit circuit with the one we constructed earlier\n",
    "\n",
    "Note that Qiskit uses Hamiltonian $C = \\frac{1}{2}\\sum_{ij\\in E}(-Z_iZ_j)$, so we will have to add a factor of $\\frac{1}{2}$ to our parameters gamma to match our circuit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "collapsed": true
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<pre style=\"word-wrap: normal;white-space: pre;background: #fff0;line-height: 1.1;font-family: &quot;Courier New&quot;,Courier,monospace\">     ┌───┐                                                  ┌───────────┐     »\n",
       "q_0: ┤ H ├──■──────────────■────■────────────────────────■──┤ RX(2pi/3) ├─────»\n",
       "     ├───┤  │              │    │    ┌───┐   ┌────────┐  │  └───┬───┬───┘     »\n",
       "q_1: ┤ H ├──┼──────────────┼────┼────┤ X ├───┤ RZ(pi) ├──┼──────┤ X ├─────────»\n",
       "     ├───┤  │              │    │    └─┬─┘   └────────┘  │      └─┬─┘    ┌───┐»\n",
       "q_2: ┤ H ├──┼──────────────┼────┼──────┼─────────────────┼────────┼──────┤ X ├»\n",
       "     ├───┤┌─┴─┐┌────────┐┌─┴─┐  │      │                 │        │      └─┬─┘»\n",
       "q_3: ┤ H ├┤ X ├┤ RZ(pi) ├┤ X ├──┼──────■─────────────────┼────────■────────■──»\n",
       "     ├───┤└───┘└────────┘└───┘┌─┴─┐┌────────┐          ┌─┴─┐                  »\n",
       "q_4: ┤ H ├────────────────────┤ X ├┤ RZ(pi) ├──────────┤ X ├──────────────────»\n",
       "     └───┘                    └───┘└────────┘          └───┘                  »\n",
       "c: 5/═════════════════════════════════════════════════════════════════════════»\n",
       "                                                                              »\n",
       "«                                                                          \n",
       "«q_0: ─────────────────────────────────────────────────────────────────────\n",
       "«     ┌───┐┌────────┐     ┌───┐┌───────────┐                               \n",
       "«q_1: ┤ X ├┤ RZ(pi) ├─────┤ X ├┤ RX(2pi/3) ├───────────────────────────────\n",
       "«     └─┬─┘├────────┤┌───┐└─┬─┘└───┬───┬───┘  ┌────────┐ ┌───┐┌───────────┐\n",
       "«q_2: ──┼──┤ RZ(pi) ├┤ X ├──┼──────┤ X ├──────┤ RZ(pi) ├─┤ X ├┤ RX(2pi/3) ├\n",
       "«       │  └────────┘└─┬─┘  │      └─┬─┘    ┌─┴────────┴┐└─┬─┘└───────────┘\n",
       "«q_3: ──┼──────────────■────┼────────┼──────┤ RX(2pi/3) ├──┼───────────────\n",
       "«       │                   │        │      └───────────┘  │  ┌───────────┐\n",
       "«q_4: ──■───────────────────■────────■─────────────────────■──┤ RX(2pi/3) ├\n",
       "«                                                             └───────────┘\n",
       "«c: 5/═════════════════════════════════════════════════════════════════════\n",
       "«                                                                          </pre>"
      ],
      "text/plain": [
       "     ┌───┐                                                  ┌───────────┐     »\n",
       "q_0: ┤ H ├──■──────────────■────■────────────────────────■──┤ RX(2pi/3) ├─────»\n",
       "     ├───┤  │              │    │    ┌───┐   ┌────────┐  │  └───┬───┬───┘     »\n",
       "q_1: ┤ H ├──┼──────────────┼────┼────┤ X ├───┤ RZ(pi) ├──┼──────┤ X ├─────────»\n",
       "     ├───┤  │              │    │    └─┬─┘   └────────┘  │      └─┬─┘    ┌───┐»\n",
       "q_2: ┤ H ├──┼──────────────┼────┼──────┼─────────────────┼────────┼──────┤ X ├»\n",
       "     ├───┤┌─┴─┐┌────────┐┌─┴─┐  │      │                 │        │      └─┬─┘»\n",
       "q_3: ┤ H ├┤ X ├┤ RZ(pi) ├┤ X ├──┼──────■─────────────────┼────────■────────■──»\n",
       "     ├───┤└───┘└────────┘└───┘┌─┴─┐┌────────┐          ┌─┴─┐                  »\n",
       "q_4: ┤ H ├────────────────────┤ X ├┤ RZ(pi) ├──────────┤ X ├──────────────────»\n",
       "     └───┘                    └───┘└────────┘          └───┘                  »\n",
       "c: 5/═════════════════════════════════════════════════════════════════════════»\n",
       "                                                                              »\n",
       "«                                                                          \n",
       "«q_0: ─────────────────────────────────────────────────────────────────────\n",
       "«     ┌───┐┌────────┐     ┌───┐┌───────────┐                               \n",
       "«q_1: ┤ X ├┤ RZ(pi) ├─────┤ X ├┤ RX(2pi/3) ├───────────────────────────────\n",
       "«     └─┬─┘├────────┤┌───┐└─┬─┘└───┬───┬───┘  ┌────────┐ ┌───┐┌───────────┐\n",
       "«q_2: ──┼──┤ RZ(pi) ├┤ X ├──┼──────┤ X ├──────┤ RZ(pi) ├─┤ X ├┤ RX(2pi/3) ├\n",
       "«       │  └────────┘└─┬─┘  │      └─┬─┘    ┌─┴────────┴┐└─┬─┘└───────────┘\n",
       "«q_3: ──┼──────────────■────┼────────┼──────┤ RX(2pi/3) ├──┼───────────────\n",
       "«       │                   │        │      └───────────┘  │  ┌───────────┐\n",
       "«q_4: ──■───────────────────■────────■─────────────────────■──┤ RX(2pi/3) ├\n",
       "«                                                             └───────────┘\n",
       "«c: 5/═════════════════════════════════════════════════════════════════════\n",
       "«                                                                          "
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "qc_our = get_qaoa_circuit_sv(G,[np.pi/3], [np.pi/2])\n",
    "qc_our.draw()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "scrolled": false,
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<pre style=\"word-wrap: normal;white-space: pre;background: #fff0;line-height: 1.1;font-family: &quot;Courier New&quot;,Courier,monospace\">     ┌───┐                                                                    »\n",
       "q_0: ┤ H ├────────────────────────────────────────────────────────────────────»\n",
       "     ├───┤                         ┌───┐┌────────┐     ┌───┐             ┌───┐»\n",
       "q_1: ┤ H ├─────────────────────────┤ X ├┤ RZ(pi) ├─────┤ X ├─────────────┤ X ├»\n",
       "     ├───┤┌───┐┌────────┐┌───┐┌───┐└─┬─┘├────────┤┌───┐└─┬─┘┌───────────┐└─┬─┘»\n",
       "q_2: ┤ H ├┤ X ├┤ RZ(pi) ├┤ X ├┤ X ├──┼──┤ RZ(pi) ├┤ X ├──┼──┤ RX(2pi/3) ├──┼──»\n",
       "     ├───┤└─┬─┘└────────┘└─┬─┘└─┬─┘  │  └────────┘└─┬─┘  │  └───────────┘  │  »\n",
       "q_3: ┤ H ├──┼──────────────┼────■────┼──────────────■────┼─────────────────■──»\n",
       "     ├───┤  │              │         │                   │                    »\n",
       "q_4: ┤ H ├──■──────────────■─────────■───────────────────■────────────────────»\n",
       "     └───┘                                                                    »\n",
       "«     ┌───┐┌────────┐     ┌───┐             ┌───┐┌────────┐┌───┐┌───────────┐\n",
       "«q_0: ┤ X ├┤ RZ(pi) ├─────┤ X ├─────────────┤ X ├┤ RZ(pi) ├┤ X ├┤ RX(2pi/3) ├\n",
       "«     └─┬─┘├────────┤┌───┐└─┬─┘┌───────────┐└─┬─┘└────────┘└─┬─┘└───────────┘\n",
       "«q_1: ──┼──┤ RZ(pi) ├┤ X ├──┼──┤ RX(2pi/3) ├──┼──────────────┼───────────────\n",
       "«       │  └────────┘└─┬─┘  │  └───────────┘  │              │               \n",
       "«q_2: ──┼──────────────┼────┼─────────────────┼──────────────┼───────────────\n",
       "«       │              │    │                 │              │  ┌───────────┐\n",
       "«q_3: ──┼──────────────■────┼─────────────────■──────────────■──┤ RX(2pi/3) ├\n",
       "«       │                   │  ┌───────────┐                    └───────────┘\n",
       "«q_4: ──■───────────────────■──┤ RX(2pi/3) ├─────────────────────────────────\n",
       "«                              └───────────┘                                 </pre>"
      ],
      "text/plain": [
       "     ┌───┐                                                                    »\n",
       "q_0: ┤ H ├────────────────────────────────────────────────────────────────────»\n",
       "     ├───┤                         ┌───┐┌────────┐     ┌───┐             ┌───┐»\n",
       "q_1: ┤ H ├─────────────────────────┤ X ├┤ RZ(pi) ├─────┤ X ├─────────────┤ X ├»\n",
       "     ├───┤┌───┐┌────────┐┌───┐┌───┐└─┬─┘├────────┤┌───┐└─┬─┘┌───────────┐└─┬─┘»\n",
       "q_2: ┤ H ├┤ X ├┤ RZ(pi) ├┤ X ├┤ X ├──┼──┤ RZ(pi) ├┤ X ├──┼──┤ RX(2pi/3) ├──┼──»\n",
       "     ├───┤└─┬─┘└────────┘└─┬─┘└─┬─┘  │  └────────┘└─┬─┘  │  └───────────┘  │  »\n",
       "q_3: ┤ H ├──┼──────────────┼────■────┼──────────────■────┼─────────────────■──»\n",
       "     ├───┤  │              │         │                   │                    »\n",
       "q_4: ┤ H ├──■──────────────■─────────■───────────────────■────────────────────»\n",
       "     └───┘                                                                    »\n",
       "«     ┌───┐┌────────┐     ┌───┐             ┌───┐┌────────┐┌───┐┌───────────┐\n",
       "«q_0: ┤ X ├┤ RZ(pi) ├─────┤ X ├─────────────┤ X ├┤ RZ(pi) ├┤ X ├┤ RX(2pi/3) ├\n",
       "«     └─┬─┘├────────┤┌───┐└─┬─┘┌───────────┐└─┬─┘└────────┘└─┬─┘└───────────┘\n",
       "«q_1: ──┼──┤ RZ(pi) ├┤ X ├──┼──┤ RX(2pi/3) ├──┼──────────────┼───────────────\n",
       "«       │  └────────┘└─┬─┘  │  └───────────┘  │              │               \n",
       "«q_2: ──┼──────────────┼────┼─────────────────┼──────────────┼───────────────\n",
       "«       │              │    │                 │              │  ┌───────────┐\n",
       "«q_3: ──┼──────────────■────┼─────────────────■──────────────■──┤ RX(2pi/3) ├\n",
       "«       │                   │  ┌───────────┐                    └───────────┘\n",
       "«q_4: ──■───────────────────■──┤ RX(2pi/3) ├─────────────────────────────────\n",
       "«                              └───────────┘                                 "
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from qiskit.compiler import transpile\n",
    "var_form_qiskit = QAOAVarForm(C.to_opflow(), 1)\n",
    "qc_qiskit = var_form_qiskit.construct_circuit([2*np.pi, 2*np.pi/3])\n",
    "circuit_qiskit = transpile(qc_qiskit, optimization_level=0,basis_gates=['rz', 'rx', 'h', 'cx'])\n",
    "circuit_qiskit.draw()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Yay! The two match"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Confirming two circuits produce the same unitary (up to global phase)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True\n"
     ]
    }
   ],
   "source": [
    "def equal_up_to_global_phase(sv1, sv2):\n",
    "    nnz = np.nonzero((sv1 != 0) & (sv2 != 0))  \n",
    "    phase = sv1[nnz] / sv2[nnz]\n",
    "    return np.all(np.isclose(phase, phase[0]))\n",
    "\n",
    "# to get rid of numerical errors\n",
    "def remove_vals_smaller_than_eps(a, eps=10*np.finfo(np.float).eps):\n",
    "    a.real[abs(a.real) < eps] = 0.0\n",
    "    a.imag[abs(a.imag) < eps] = 0.0\n",
    "    return a\n",
    "\n",
    "backend = Aer.get_backend(\"unitary_simulator\")\n",
    "unitary_our = remove_vals_smaller_than_eps(execute(qc_our, backend=backend).result().get_unitary())\n",
    "unitary_qiskit = remove_vals_smaller_than_eps(execute(qc_qiskit, backend=backend).result().get_unitary())\n",
    "print(equal_up_to_global_phase(unitary_our, unitary_qiskit))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Q: what if I want to use the Hamiltonian that I've built in numpy?\n",
    "\n",
    "You can do that too! All you need to do is to convert it to Qiskit object"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "type of C_numpy is <class 'numpy.ndarray'>\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<qiskit.aqua.algorithms.minimum_eigen_solvers.qaoa.var_form.QAOAVarForm at 0x7fab70a8c5e0>"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print(f\"type of C_numpy is {type(C_numpy)}\")\n",
    "\n",
    "from qiskit.aqua.operators import MatrixOperator\n",
    "C_converted = MatrixOperator(C_numpy)\n",
    "\n",
    "# Now we can use it!\n",
    "var_form = QAOAVarForm(C_converted.to_opflow(), p)\n",
    "var_form"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Q: Okay, but we still have to wrap it into objective and pass it to optimizer, which is a lot of work! Can Qiskit do that for me too?\n",
    "\n",
    "Yes! The reason I'm leaving it to the end of the tutorial is that it looks nice and clean, but it might not be flexible enough for a researcher like yourself. Moreover, the cleanliness of the interface hides a lot of underlying complexity (as you've seen earlier) and introduces some limitations. One of such limitations is that you have to work with optimizers supported by Qiskit (though it's fairly easy to add a new one). The main point is that qiskit is open source, so you can always modify it to suit your needs! "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Evolved Hamiltonian is not composed of only MatrixOps, converting to Matrix representation, which can be expensive.\n",
      "Evolved Hamiltonian is not composed of only MatrixOps, converting to Matrix representation, which can be expensive.\n",
      "Evolved Hamiltonian is not composed of only MatrixOps, converting to Matrix representation, which can be expensive.\n",
      "Evolved Hamiltonian is not composed of only MatrixOps, converting to Matrix representation, which can be expensive.\n",
      "Evolved Hamiltonian is not composed of only MatrixOps, converting to Matrix representation, which can be expensive.\n",
      "Evolved Hamiltonian is not composed of only MatrixOps, converting to Matrix representation, which can be expensive.\n",
      "Evolved Hamiltonian is not composed of only MatrixOps, converting to Matrix representation, which can be expensive.\n",
      "Evolved Hamiltonian is not composed of only MatrixOps, converting to Matrix representation, which can be expensive.\n",
      "Evolved Hamiltonian is not composed of only MatrixOps, converting to Matrix representation, which can be expensive.\n",
      "Evolved Hamiltonian is not composed of only MatrixOps, converting to Matrix representation, which can be expensive.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "final maxcut energy: (-5.3056640625+0j)\n",
      "solution: [0 0 1 1 1]\n"
     ]
    }
   ],
   "source": [
    "from qiskit.aqua.algorithms.minimum_eigen_solvers import QAOA\n",
    "from qiskit.aqua.components.optimizers import COBYLA\n",
    "from qiskit.optimization.applications.ising import max_cut\n",
    "\n",
    "cobyla = COBYLA(maxiter=1000)\n",
    "qaoa = QAOA(C_converted, cobyla, p=5)\n",
    "backend = Aer.get_backend('qasm_simulator')\n",
    "result = qaoa.run(backend)\n",
    "\n",
    "x = sample_most_likely(result['eigenstate'])\n",
    "print('final maxcut energy:', result['eigenvalue'] + offset)\n",
    "print('solution:', x)"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}