tutorials/02_QAOA.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# QAOA - Quantum Approximate Optimization Algorithm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this notebook I would like to show you one of the most prominent algorithms in the quantum computing in the last few years."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this tutorial you will learn the following:\n",
    "- How to use pyQuil - python platform for quantum computing.\n",
    "- How to solve MaxCut problem using QAOA implementation in pyQuil.\n",
    "- What does it mean in practice that quantum computing is probabilistic.\n",
    "- Basic information about how QAOA works."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is probably the most theoretical notebook in this series, especially the second part. If you get discouraged by this and don't get everything, don't worry. Take your time, try reading different resources or just go to the next part and come back later. These ideas are that difficult, but are pretty complicated at times and they may just need some time to sink in :)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Prerequisites"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before we can go into QAOA itself, you should prepare a couple of things:\n",
    "1. Download and install [Forest SDK](https://qcs.rigetti.com/sdk-downloads)\n",
    "2. Install [pyQuil](https://github.com/rigetticomputing/pyquil) and [Grove](https://github.com/rigetticomputing/grove)\n",
    "3. Learn the basics of quantum computing by yourself - there are plenty of great resources about the basics in the web and I decided I don't want to reinvent the wheel in this tutorial. I recommend going to [pyQuil documentation](http://pyquil.readthedocs.io/en/stable/).\n",
    "4. You can also check out [this tutorial](https://github.com/markf94/rigetti_training_material) by guys from [ProteinQure](https://proteinqure.com). It was designed for the workshops, not self-learning, but there are many hands-on exercises and it will be a good test of your pyQuil skills :)\n",
    "\n",
    "Even though I tried to make this tutorial as accessible to the not-quantum-computing people as possible, you probably need some level of familiarity with the following concepts:\n",
    "- qubits\n",
    "- quantum gates\n",
    "- superposition\n",
    "- entanglement\n",
    "\n",
    "You don't need to go very in-depth - I think spending about 30 minutes with the pyQuil docs should be enough for the beginning. And if something is still confusing later on you can always come back to the documentation when you need it :)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Terminology"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need this line to have a proper [bra-ket](https://en.wikipedia.org/wiki/Bra–ket_notation) formatting in markdown.\n",
    "$$\n",
    "\\newcommand{\\ket}[1]{\\left| #1 \\right\\rangle}\n",
    "\\newcommand{\\bra}[1]{\\left\\langle #1 \\right|}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I try to use as easy language as possible and by easy I mean not very specific neither for physicists, programmers nor quantum computing people. But sometimes it's just not possible, so here I will try to explain you some important concepts in layman terms. These explanation won't give you deep understanding, rather a general overview."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`Operator` - well... if you are not very mathematically accurate, you can treat \"operator\" as a fancy name for a matrix. If you are - I'm not the best person to explain it to you ;) But in essence, operators represents some kind of action on a quantum state and there is a matrix associated with every operator, so these terms are sometimes used interchangeably."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`Hamiltonian` is an operator describing the quantum system. It's eigenvalues represent the energy values that the system can have."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`Pauli matrix` is one of the three 2x2 matrices commonly used in quantum computing. They have very interesting properties ([wiki](https://en.wikipedia.org/wiki/Pauli_matrices)) and they correspond to some very basic operations on qubits. You will see a lot of them while doing quantum computing.\n",
    "\n",
    "- $\\sigma_x$ -> X gate. It flips $\\ket 0$ and $\\ket 1$, so it's equivalent to logical NOT gate.\n",
    "- $\\sigma_y$ -> Y gate. It maps $\\ket 0$ into $i \\ket 1$ and $\\ket 1$ into $-i \\ket 0$\n",
    "- $\\sigma_z$ -> Z gate. It doesn't affect $\\ket 0$ state, but changes $\\ket 1$ -> - $\\ket 1$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## QAOA - intro"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "QAOA is an algorithm for solving a broad range of optimization problems using NISQ (Noisy Intermediate-Scale Quantum) devices. \n",
    "\n",
    "Let's start easy and instead of going straight to the TSP we will solve MaxCut problem. This part is based on the [Rigetti's tutorial on QAOA](https://grove-docs.readthedocs.io/en/latest/qaoa.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### MaxCut - explanation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the MaxCut problem we start with a graph and we want to divide it into two subgraphs in such way, that the edges between them would have the highest possible sum."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's take look on an example. We will color the nodes with blue and red to distinguish into which group they fall - let's say blue represents 0 and red represents 1.\n",
    "\n",
    "We start easy with a graph with two nodes and one connection:\n",
    "![](images/02_2_nodes.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can color it in the following ways:\n",
    "![](images/02_2_nodes_cuts.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, a little more advanced example:\n",
    "![](images/02_3_nodes_cuts.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To understand it better I recommend drawing a 4-node graph on a paper and try to solve it by hand.\n",
    "\n",
    "If you still don't get it, check out [this section of the Rigetti's tutorial](https://grove-docs.readthedocs.io/en/latest/qaoa.html#our-first-np-hard-problem)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### QAOA solves MaxCut - quick take"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's take a graph with four nodes.\n",
    "We encode it as a list of tuples, where each tuple represents an edge between given nodes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "first_graph = [(0, 1), (0, 2), (0, 3)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![](images/02_4_nodes.png)\n",
    "The optimal solution for this graph is either 1000 or 0111:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's try to solve MaxCut problem for this graph using maxcut_solver from the grove library.\n",
    "\n",
    "Keep in mind that runtime of this algorithm are non-deterministic - it should finish running under a minute, but sometimes it takes longer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from grove.pyqaoa.maxcut_qaoa import maxcut_qaoa\n",
    "import pyquil.api as api\n",
    "qvm_connection = api.QVMConnection()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture \n",
    "#%%capture supresses printing. \n",
    "#get_angles() prints out a lot of stuff that we don't care about right now.\n",
    "maxcut_solver = maxcut_qaoa(graph=first_graph)\n",
    "betas, gammas = maxcut_solver.get_angles()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have the values of betas and gammas we can run the actual program.\n",
    "\n",
    "Don't worry about understanding what betas and gammas are for now. We will get back to it later."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "angles = np.hstack((betas, gammas))\n",
    "param_prog = maxcut_solver.get_parameterized_program()\n",
    "prog = param_prog(angles)\n",
    "qubits = [0, 1, 2, 3]\n",
    "measurements = qvm_connection.run_and_measure(prog, qubits, trials=1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And see what are the results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[((1, 0, 0, 0), 231),\n",
       " ((0, 1, 1, 1), 220),\n",
       " ((1, 1, 0, 0), 84),\n",
       " ((1, 0, 0, 1), 84),\n",
       " ((0, 0, 1, 1), 81),\n",
       " ((1, 0, 1, 0), 75),\n",
       " ((0, 1, 1, 0), 68),\n",
       " ((0, 1, 0, 1), 66),\n",
       " ((0, 0, 0, 0), 19),\n",
       " ((1, 1, 1, 1), 17),\n",
       " ((0, 0, 1, 0), 14),\n",
       " ((0, 1, 0, 0), 11),\n",
       " ((0, 0, 0, 1), 10),\n",
       " ((1, 0, 1, 1), 7),\n",
       " ((1, 1, 0, 1), 7),\n",
       " ((1, 1, 1, 0), 6)]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from collections import Counter\n",
    "measurements = [tuple(measurement) for measurement in measurements]\n",
    "measurements_counter = Counter(measurements)\n",
    "measurements_counter.most_common()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So as we can see, the correct results are at the top of our lists - so everything seems good.\n",
    "\n",
    "But why?\n",
    "\n",
    "We will go through the code step by step in the next section"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### QAOA solves MaxCut - details"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Finding the right angles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WARNING: Fast method for expectation will be used. Noise\n",
      "                     models will be ineffective\n",
      "\tParameters: [1.54679531 6.139004  ] \n",
      "\tE => -1.4940854241926789\n",
      "\tParameters: [1.54679531 6.139004  ] \n",
      "\tE => -1.497553504498376\n",
      "\tParameters: [1.54679531 6.139004  ] \n",
      "\tE => -1.5152091246046904\n",
      "\tParameters: [1.48395675 5.93756793] \n",
      "\tE => -1.6630606037921771\n",
      "\tParameters: [1.48395675 5.93756793] \n",
      "\tE => -1.6335549243815815\n",
      "\tParameters: [1.40178325 5.39081289] \n",
      "\tE => -2.0092685257645435\n",
      "\tParameters: [1.30269167 5.30927924] \n",
      "\tE => -2.216995479920419\n",
      "\tParameters: [1.08879888 4.17500233] \n",
      "\tE => -2.261790080569641\n",
      "\tParameters: [1.08879888 4.17500233] \n",
      "\tE => -2.233258800349405\n",
      "\tParameters: [1.08879888 4.17500233] \n",
      "\tE => -2.2509385947021006\n",
      "\tParameters: [1.08879888 4.17500233] \n",
      "\tE => -2.2570753231984773\n",
      "\tParameters: [1.08879888 4.17500233] \n",
      "\tE => -2.255386716253838\n",
      "\tParameters: [1.11283539 4.14986029] \n",
      "\tE => -2.242086582256904\n",
      "\tParameters: [1.11283539 4.14986029] \n",
      "\tE => -2.281399329541506\n",
      "\tParameters: [1.17084493 4.33609311] \n",
      "\tE => -2.2712914971049964\n",
      "\tParameters: [1.13687189 4.12471825] \n",
      "\tE => -2.297163937095549\n",
      "\tParameters: [1.13687189 4.12471825] \n",
      "\tE => -2.293988488017738\n",
      "\tParameters: [1.1609084  4.09957621] \n",
      "\tE => -2.3040751738756113\n",
      "\tParameters: [1.1609084  4.09957621] \n",
      "\tE => -2.3093140630074114\n",
      "\tParameters: [1.1609084  4.09957621] \n",
      "\tE => -2.310056495263365\n",
      "\tParameters: [1.1849449  4.07443417] \n",
      "\tE => -2.3097331072788796\n",
      "\tParameters: [1.1849449  4.07443417] \n",
      "\tE => -2.3153114845422595\n",
      "\tParameters: [1.17279654 4.10257319] \n",
      "\tE => -2.316293546851724\n",
      "\tParameters: [1.17279654 4.10257319] \n",
      "\tE => -2.3162662044047053\n",
      "\tParameters: [1.18108349 4.09101161] \n",
      "\tE => -2.316416885858272\n",
      "\tParameters: [1.17929381 4.10469866] \n",
      "\tE => -2.316450415597711\n",
      "\tParameters: [1.17649259 4.10021416] \n",
      "\tE => -2.3164730994694076\n",
      "Values of betas: [1.17649259]\n",
      "Values of gammas: [4.10021416]\n"
     ]
    }
   ],
   "source": [
    "# We initialize the maxcut_qaoa object with our graph\n",
    "maxcut_solver = maxcut_qaoa(graph=first_graph)\n",
    "# The QAOA algorithm tries to find the optimal values of betas and gammas.\n",
    "# This line is where all the optimization takes place.\n",
    "betas, gammas = maxcut_solver.get_angles()\n",
    "print(\"Values of betas:\", betas)\n",
    "print(\"Values of gammas:\", gammas)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This time we have not supressed the output and with every iteration of the algorithm you can see two lines being printed.\n",
    "<br>\n",
    "The first one is parameters - these are values of betas and gammas.\n",
    "<br>\n",
    "The second one is the cost value - it's what we try to minize. As you can probably see, the minimization process isn't ideally smooth. Value is generally getting smaller and smaller, but not necessarily every step yields an improvement.\n",
    "<br>\n",
    "We will deal with what exactly this energy means later on."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Creating a circuit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1.17649259 4.10021416]\n"
     ]
    }
   ],
   "source": [
    "# We create an array of angles with correct format\n",
    "angles = np.hstack((betas, gammas))\n",
    "print(angles)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "H 0\n",
      "H 1\n",
      "H 2\n",
      "H 3\n",
      "CNOT 0 1\n",
      "RZ(4.100214159979846) 1\n",
      "CNOT 0 1\n",
      "X 0\n",
      "PHASE(2.050107079989923) 0\n",
      "X 0\n",
      "PHASE(2.050107079989923) 0\n",
      "CNOT 0 2\n",
      "RZ(4.100214159979846) 2\n",
      "CNOT 0 2\n",
      "X 0\n",
      "PHASE(2.050107079989923) 0\n",
      "X 0\n",
      "PHASE(2.050107079989923) 0\n",
      "CNOT 0 3\n",
      "RZ(4.100214159979846) 3\n",
      "CNOT 0 3\n",
      "X 0\n",
      "PHASE(2.050107079989923) 0\n",
      "X 0\n",
      "PHASE(2.050107079989923) 0\n",
      "H 0\n",
      "RZ(-2.3529851869388656) 0\n",
      "H 0\n",
      "H 1\n",
      "RZ(-2.3529851869388656) 1\n",
      "H 1\n",
      "H 2\n",
      "RZ(-2.3529851869388656) 2\n",
      "H 2\n",
      "H 3\n",
      "RZ(-2.3529851869388656) 3\n",
      "H 3\n",
      "\n",
      "Number of gates: 37\n"
     ]
    }
   ],
   "source": [
    "# We take a template for quil program from the maxcut_solver.\n",
    "param_prog = maxcut_solver.get_parameterized_program()\n",
    "# We initialize this program with the angles we have found\n",
    "prog = param_prog(angles)\n",
    "# Now we can print the program. \n",
    "# Some of the values you see here are the angles we calculated earlier.\n",
    "print(prog)\n",
    "print(\"Number of gates:\", len(prog))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Running the quantum program"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# These are just the ids of qubits we want to use.\n",
    "# It's not very important if you don't use the real QPU.\n",
    "qubits = [0, 1, 2, 3]\n",
    "# Here we connect to the Forest API and run our program there.\n",
    "# We do that 1000 times and after each one we measure the output.\n",
    "measurements = qvm_connection.run_and_measure(prog, qubits, trials=1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Analyzing the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[((1, 0, 0, 0), 232),\n",
       " ((0, 1, 1, 1), 223),\n",
       " ((0, 1, 0, 1), 86),\n",
       " ((0, 0, 1, 1), 83),\n",
       " ((0, 1, 1, 0), 78),\n",
       " ((1, 0, 0, 1), 74),\n",
       " ((1, 1, 0, 0), 63),\n",
       " ((1, 0, 1, 0), 62),\n",
       " ((1, 1, 1, 1), 34),\n",
       " ((0, 0, 0, 0), 27),\n",
       " ((1, 1, 0, 1), 9),\n",
       " ((0, 1, 0, 0), 8),\n",
       " ((0, 0, 0, 1), 7),\n",
       " ((1, 1, 1, 0), 7),\n",
       " ((0, 0, 1, 0), 4),\n",
       " ((1, 0, 1, 1), 3)]"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Since list of 1000 elements is hard to analyze, we use Counter\n",
    "from collections import Counter\n",
    "# This is just a hack - we can't use Counter on a list of lists but we can on a list of tuples.\n",
    "measurements = [tuple(measurement) for measurement in measurements]\n",
    "measurements_counter = Counter(measurements)\n",
    "# This line gives us the results in the diminishing order\n",
    "measurements_counter.most_common()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Probabilistic nature of quantum computing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One of the most basic rules of the quantum mechanics is that the act of observation changes the state of the system that we observe. \n",
    "<br>\n",
    "Even though our algorithm produces some state which is a combination of all the states you see above, we can't measure it directly. What we can do (and we did) is to measure it repeatedly and get some distribution of the results.\n",
    "<br>\n",
    "Fortunately, we use a simulator and we can look how the state we produced looks exactly:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(0.0813147891-0.1300523955j)|0000> + (0.4653262596+0.1064676374j)|0001> + (0.0804387018-0.0324794329j)|0010> + (-0.2535746518+0.1050939054j)|0011> + (0.0804387018-0.0324794329j)|0100> + (-0.2535746518+0.1050939054j)|0101> + (-0.2535746518+0.1050939054j)|0110> + (0.0804387018-0.0324794329j)|0111> + (0.0804387018-0.0324794329j)|1000> + (-0.2535746518+0.1050939054j)|1001> + (-0.2535746518+0.1050939054j)|1010> + (0.0804387018-0.0324794329j)|1011> + (-0.2535746518+0.1050939054j)|1100> + (0.0804387018-0.0324794329j)|1101> + (0.4653262596+0.1064676374j)|1110> + (0.0813147891-0.1300523955j)|1111>\n"
     ]
    }
   ],
   "source": [
    "wf = qvm_connection.wavefunction(prog)\n",
    "print(wf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is hard to read and interpret, so let's print it in a different form."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Probability amplitudes for all the possible states:\n",
      "0000 (0.08131478909571628-0.13005239547041625j)\n",
      "0001 (0.4653262595505163+0.10646763736089705j)\n",
      "0010 (0.08043870178752226-0.03247943285344676j)\n",
      "0011 (-0.253574651788108+0.10509390542847408j)\n",
      "0100 (0.08043870178752226-0.032479432853446726j)\n",
      "0101 (-0.253574651788108+0.10509390542847413j)\n",
      "0110 (-0.253574651788108+0.10509390542847415j)\n",
      "0111 (0.08043870178752228-0.03247943285344669j)\n",
      "1000 (0.08043870178752228-0.03247943285344669j)\n",
      "1001 (-0.253574651788108+0.10509390542847415j)\n",
      "1010 (-0.253574651788108+0.10509390542847413j)\n",
      "1011 (0.08043870178752226-0.032479432853446726j)\n",
      "1100 (-0.253574651788108+0.10509390542847408j)\n",
      "1101 (0.08043870178752226-0.03247943285344676j)\n",
      "1110 (0.4653262595505163+0.10646763736089705j)\n",
      "1111 (0.08131478909571628-0.13005239547041625j)\n"
     ]
    }
   ],
   "source": [
    "print(\"Probability amplitudes for all the possible states:\")\n",
    "for state_index in range(maxcut_solver.nstates):\n",
    "    print(maxcut_solver.states[state_index], wf[state_index])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These numbers are still not very easy to interpret, but we can calculate the actual probabilities based on this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Probabilities of measuring given states:\n",
      "0000 0.023525720493274362\n",
      "0001 0.22786388563248594\n",
      "0010 0.007525298303743494\n",
      "0011 0.07534483298766928\n",
      "0100 0.007525298303743491\n",
      "0101 0.07534483298766928\n",
      "0110 0.0753448329876693\n",
      "0111 0.007525298303743493\n",
      "1000 0.007525298303743493\n",
      "1001 0.0753448329876693\n",
      "1010 0.07534483298766928\n",
      "1011 0.007525298303743491\n",
      "1100 0.07534483298766928\n",
      "1101 0.007525298303743494\n",
      "1110 0.22786388563248594\n",
      "1111 0.023525720493274362\n"
     ]
    }
   ],
   "source": [
    "print(\"Probabilities of measuring given states:\")\n",
    "states_with_probs = []\n",
    "for state_index in range(maxcut_solver.nstates):\n",
    "    states_with_probs.append([maxcut_solver.states[state_index], np.real(np.conj(wf[state_index])*wf[state_index])])\n",
    "    print(states_with_probs[-1][0], states_with_probs[-1][1])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These results correspond nicely with the values we got from the measurement. Not ideally, but you can increase the `trial` parameter to get closer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def plot_state_histogram(states_with_probs):\n",
    "    states = np.array(states_with_probs)[:,0]\n",
    "    probs = np.array(states_with_probs)[:,1].astype(float)\n",
    "    n = len(states_with_probs)\n",
    "    plt.barh(range(n), probs, tick_label=states)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAD8CAYAAACRkhiPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAEyZJREFUeJzt3X/sZXdd5/HnyyGtEtNhyjTSTJEZ4ixk2G6iXgtqdNHizJTN2rrCpk0WRnciiUCIISSWoKkBm7hi1FQRfzYtJkv5kSizODIZC27BtDJ3oLZMydhxqDIN2MiU8kuKhff+cU+7X7+93/me7/d77znn+73PR3LzvfdzPufz/Zwzd+Y9555zXjdVhSRpMX1b3xOQJPXHIiBJC8wiIEkLzCIgSQvMIiBJC8wiIEkLzCIgSQvMIiBJC8wiIEkL7Bl9T2A1O3furN27d/c9DUnaVE6ePPkvVXXZav0GXwR2797NeDzuexqStKkk+cc2/fw4SJIWmEVAkhaYRUCSFphFQJIWmEVAkhaYRUCSFphFQJIWmEVAkhbY4G8Wu//hx9h941/0OoeHfu2/9Pr7JWlePBKQpAW2ahFIcmuSR5J8aknbK5OcSvKtJKMl7c9O8pEkX0nyu8vGuTnJZ5N8ZbabIElarzZHArcBB5e1fQr4b8Bdy9q/Dvwy8KYp4/wf4Ko1zk+SNEernhOoqruS7F7W9mmAJMv7fhX4WJLvmTLOPdPWkST1Z5DnBJK8Jsk4yfibX3us7+lI0pY1yCJQVX9YVaOqGm175va+pyNJW9Ygi4AkqRsWAUlaYG0uEX03cDfwgiTnkhxO8lNJzgE/CPxFkmNL+j8E/CbwM03/fU37rzfrPLNp/5U5bI8kaQ1SVX3P4YJGo1H59ZKStDZJTlbVaLV+fhwkSQvMIiBJC8wAOUkaoK6CK7vMDvr+JPcnOZPklnjrsCT1rsvsoHcCPwfsbR7Lx5QkdWzVIlBVdwHnl7V9uqpOT+n71ar6GJNi8JQklwOXVNU9Nbkc6V3AdRuauSRpw7o6MbwLOLfk9bmmTZLUo0FeHWSAnCR1o6si8DBwxZLXVzRtUxkgJ0nd6KQIVNXngC8leUlzVdCrgQ908bslSStb9T6BJjvopcDOJvvnJiYnin8HuIxJdtC9VXWg6f8QcAlwUZLrgP1V9QDwWiZXGn0H8JfNQ5LUozbfLHbDCov+bIX+u1doHwP/sfXMJElzN/g7hq/ctZ1xR3fOSdKiGeTVQZKkbgz+SKBtdlBXORuStJV4JCBJC2ymAXLNsjc3IXGnkxy40DiSpH7NNECu+SrJ64EXNev8XpJtFxhHktSjmQbIAdcCd1TV41X1GeAMcNVK40iS+jXrcwK7gM8ueb2uoDizgySpG4M8MWx2kCR1Y9ZF4GHguUteXzAoTpLUr1kXgSPA9UkuTrKHyTeIfXzGv0OSNCNtLhF9N3A38IIk55IcTvJTTZjcDzIJkDsGUFWngPcCDwAfAl5XVd9caZz5bJIkqa1Mvu1xuEajUY3H476nIUmbSpKTVTVard8gTwxLkrphEZCkBbZlAuTUH8P7pM2r1ZHADPODDjZtZ5LcOLvNkCStR9uPg25jg/lBTYbQO4BrgH3ADU1fSVJPWn0cVFV3Jdm9rO3TAJPvjf93nsoPAj6T5Kn8IOBMVZ1t1ruj6fvAeicvSdqYeZwYXik/aCa5QpKk2Rnk1UEGyElSN+ZRBFbKD2qdK2SAnCR1Yx5FYKX8oBPA3iR7klzE5OTxkTn8fklSS61ODDe5Py8FdjaZQTcx+YKY3wEuY5IfdG9VHaiqU0mezA96gn+fH/R64BiwDbi1yRqSJPWk7dVBN6yw6M9W6H8zcPOU9qPA0dazkyTN1eDvGL5y13bG3pEqSXMxyKuDJEndGPyRQNvsIPNrJGntPBKQpAXW5pvFpoXHXZrkeJIHm587mvYXJrk7yeNJ3rTaOJKkfrU5EriNp4fH3QjcWVV7gTub1zC5bPQNwG+0HEeS1KNVi0BV3cXkH/elrgVub57fDlzX9H2kqk4A/9ZyHElSj9Z7TuC7qupzzfPPA981o/kAZgdJUlc2fGK4Jt9UP9Nvqzc7SJK6sd4i8M9JLgdofj4yuylJkrqy3iJwBDjUPD8EfGA205EkdanNJaLvBu4GXpDkXJLDwK8BP5HkQeBlzWuSPKcJmHsj8EtN/0suMI4kqUeZfKQ/XKPRqMbjcd/TkKRNJcnJqhqt1s87hiVpgVkEJGmBbZkAOfXH8D5p82p1JDDD/KCDSU4nOZPkxuW/R5LUrbYfB93GBvODkmwD3gFcA+wDbkiyb33TliTNQqsiMKP8oKuAM1V1tqq+AdzRjCFJ6slGTgyvNT9oF/DZJa/PNW2SpJ7M5OqgWecHGSAnSd3YSBFYa37Qw8Bzl7y+oml7GgPkJKkbGykCa80POgHsTbInyUXA9c0YkqSetLpPoMn9eSmws8kGuolJXtB7mwygfwT+e9P3OcAYuAT4VpJfAPZV1ZeSvB44BmwDbq2qUzPeHknSGrQqAlV1wwqLrp7S9/NMPuqZNs5R4Gjr2UmS5mrwdwxfuWs7Y+9IlaS5MDtIkhbY4I8EzA4aPrODpM3LIwFJWmDzCJBLkluakLj7knzfknU+lOSLST44+02RJK3VPALkrgH2No/XAO9css7bgVetd7KSpNmaeYBc0/6umrgHeNaTdxZX1Z3Alzc8a0nSTMwjQG7DQXFmB0lSNwYZIGd2kCR1Yx4Bcq2D4iRJ/ZpHgNwR4NXNVUIvAR5b8rGRJGlAZh4gxyQb6OXAGeBrwM8uGeejwAuB72zGOVxVx2azKZKktcrk4/zhGo1GNR6P+56GJG0qSU5W1Wi1ft4xLEkLzCIgSQtsywTIGWImSWvXdXbQoab/g0kOTftdkqTudJYdlORSJlcVvRi4CrjpycIhSepHl9lBB4DjVXW+qh4FjvP0wiJJ6lCX2UEbzhSSJM3WILODDJCTpG50mR3UOlPIADlJ6kaX2UHHgP1JdjQnhPc3bZKknnSWHVRV55O8DTjR9HtrVS0/2SxJ6lCrIlBVN6yw6OopfQt43Qrj3Arc2np2kqS5Gvwdw1fu2s7Yu4ElaS7MDpKkBTb4IwGzgyRpfjwSkKQF1jZA7mCS000o3I1N2+ub15Vk55K+L0xyd5LHk7xp2ThPC6KTJPVn1SKQZBvwDibBcPuAG5LsA/4GeBmTy0OXOg+8AfiNKcPdhnlBkjQYbY4ErgLOVNXZqvoGcAdwbVV9sqoeWt65qh6pqhPAv01ZNi2ITpLUkzZFoPPgN7ODJKkbgzwxbHaQJHWjTRFoHfwmSdpc2hSBE8DeJHuSXARczyQkTpK0ya1aBKrqCeD1TBI/Pw28t6pOJXlDEyZ3BXBfkj8GSPKcpv2NwC8lOZfkkmbZu4G7gRc07Yfns1mSpDYyyXsbrtFoVOPxuO9pSNKmkuRkVY1W6zfIE8OSpG5YBCRpgW2ZADn1x/A+afPqOjvoaeNIkvrTWXbQBcaRJPWky+ygqeNsbPqSpI3oMjuo8wwiSdKFDfLqIAPkJKkbXWYHtR7HADlJ6kaX2UFmEEnSwHSWHbTSOPPZLElSG61uFquqo8DRZW23ALdM6ft5JoWh1TiSpP4M/o7hK3dtZ+wdqZI0F4O8OkiS1I3BHwmYHTR8ZgdJm5dHApK0wOYRIJcktzTL7kvyfUuWfSjJF5N8cPabIklaq3kEyF0D7G0erwHeuWTZ24FXbXzakqRZmHmAHJNQuHfVxD3As5JcDlBVdwJfntHcJUkbNI8AuQ0HxZkdJEndGOSJYbODJKkb8wiQm1XgnCRpzuYRIHcEeHVzldBLgMeq6nMzmKskacZmHiDHJBvoLHAG+CPgtU+OleSjwPuAq5tguQMz3RpJ0pqkqvqewwWNRqMaj8d9T0OSNpUkJ6tqtFq/QZ4YliR1wyIgSQtsywTIGWImSWvXdXbQoSQPNo9Ds98cSdJarHoksCQ76CeY3P17IskRJtlBHwT+etkqS7ODXswkO+jFSS4FbgJGQAEnkxypqkdnsymSpLXqMjvoAHC8qs43//AfBw7OZjMkSevRZXbQhjOFJEmzNcirgwyQk6RudJkd1HocA+QkqRtdZgcdA/Yn2ZFkB7C/aZMk9aSz7KCqOg+8jUlROQG8tWmTJPVk8NlBF1++ty4/9Nur9vNmMUn6/9pmBw3+juErd21n7D/wkjQXg7w6SJLUjcEfCbTNDlJ//ChO2rw8EpCkBbaRALk9Sf62aXtPc/koSX40ySeSPJHkFcvG+VCSLyb54Ow3RZK0VqsWgSUBctcA+4AbkuwD/hfwW1X1PcCjwOFmlX8Cfgb431OGezvwqo1PW5I0C+sOkAN+HHh/0+d24DqAqnqoqu4DvrV8oKq6E/jyLCYuSdq4jQTIfbG5kWxp20yYHSRJ3RjkiWGzgySpGxsJkHtWkmcsa5MkbSIbCZD7CPDk1T+HgA/MZ4qSpHlZd4Ac8IvAG5OcAZ4N/AlAkh9oguVeCfxBklNPjpXko8D7gKuTnEtyYNYbJElqb/ABcqPRqMbjcd/TkKRNpW2A3CBPDEuSumERkKQFtmUC5Awxk6S16zo76FCSB5vHodlvjiRpLTrLDkpyKXAT8GImURQ3Nd81LEnqSZfZQQeA41V1vqoeBY4DB2ewDZKkdeoyO2ilcSRJPRnk1UEGyElSN7rMDlppnKcxQE6SutFldtAxYH+SHc0J4f1NmySpJ51lB1XVeeBtTIrKCeCtTZskqSetbharqqPA0WVtZ5lcObS87wkmH/VMG+dW4Na1T1OSNA+Dv2P4yl3bGXs3sCTNxSCvDpIkdWPwRwJts4MkaSvpKg/NIwFJWmDzCJC7uHl9plm+u2l/dpKPJPlKkt+d1wZJktqbR4DcYeDRpv23mn4AXwd+GXjTTLdAkrRuMw+Qa5bd3jx/P5PvE05VfbWqPsakGEiSBmAeAXJP9W+WP8bkZrLWzA6SpG4M8sSw2UGS1I15BMg91b9Zvh34wkxmK0maqXkEyB1pXtMs/3BV1eymLEmalVVvFquqJ5I8GSC3Dbi1qk4l+UXgjiS/CnySJkCu+fmnTbDceSZFA4AkDwGXABcluQ7YX1UPzHKDJEntZej/SR+NRjUej/uehiRtKklOVtVotX6DPDEsSeqGRUCSFpgBci10FeQkSV3rLDuoWfbmpv10kgPz2CBJUnudZQc161wPvAg4CPxeM7YkqSedZQc17XdU1eNV9RngDFO+nlKS1J0us4NWGkeS1JNBXh1kgJwkdaPL7KCVxnkaA+QkqRtdZgcdAa5vrh7aA+wFPj6bzZAkrUdn2UHNOu8FHgCeAF5XVd+c+RZJklprdbNYVR0Fji5rO8uUq3uq6uvAK1cY52bg5rVPU5I0D4O/Y/jKXdsZe8euJM3FIK8OkiR1wyIgSQvMIiBJC8wiIEkLzCIgSQvMIiBJC8wiIEkLzCIgSQvMIiBJCyyTbLfhSvJl4HTf8xigncC/9D2JAXK/TOd+mW4r75fnVdVlq3UafGwEcLqqRn1PYmiSjN0vT+d+mc79Mp37xY+DJGmhWQQkaYFthiLwh31PYKDcL9O5X6Zzv0y38Ptl8CeGJUnzsxmOBCRJc9JbEUhyMMnpJGeS3Dhl+cVJ3tMs/9sku5cse3PTfjrJgS7nPW/r3S9Jdif51yT3No/f73ru89Riv/xokk8keSLJK5YtO5TkweZxaPm6m9kG98s3l7xfjnQ36/lrsV/emOSBJPcluTPJ85Ys27Lvl6mqqvMHk+8q/gfg+cBFwN8B+5b1eS3w+83z64H3NM/3Nf0vBvY042zrYzsGtl92A5/qext63C+7gf8EvAt4xZL2S4Gzzc8dzfMdfW9T3/ulWfaVvrehx/3yY8Azm+c/v+Tv0ZZ9v6z06OtI4CrgTFWdrapvAHcA1y7rcy1we/P8/cDVSdK031FVj1fVZ4AzTPmu401qI/tlK1t1v1TVQ1V1H/CtZeseAI5X1fmqehQ4DhzsYtId2Mh+2cra7JePVNXXmpf3AFc0z7fy+2WqvorALuCzS16fa9qm9qmqJ4DHgGe3XHez2sh+AdiT5JNJ/m+SH5n3ZDu0kT/zRX+/XMi3JxknuSfJdbOdWq/Wul8OA3+5znU3vc1wx7Da+Rzw3VX1hSTfD/x5khdV1Zf6npgG63lV9XCS5wMfTnJ/Vf1D35PqUpL/AYyA/9z3XPrS15HAw8Bzl7y+ommb2ifJM4DtwBdarrtZrXu/NB+PfQGgqk4y+Uz0P8x9xt3YyJ/5or9fVlRVDzc/zwJ/DXzvLCfXo1b7JcnLgLcAP1lVj69l3a2kryJwAtibZE+Si5ic4Fx+dcIR4Mkz868APlyTMzdHgOubq2T2AHuBj3c073lb935JclmSbQDN/+z2MjmptRW02S8rOQbsT7IjyQ5gf9O2Fax7vzT74+Lm+U7gh4EH5jbTbq26X5J8L/AHTArAI0sWbeX3y3R9nZEGXg78PZP/sb6laXsrkz8UgG8H3sfkxO/HgecvWfctzXqngWv6Prs+hP0C/DRwCrgX+ATwX/velo73yw8w+fz2q0yOGE8tWfd/NvvrDPCzfW/LEPYL8EPA/UyunLkfONz3tnS8X/4K+Ofm78u9wJFFeL9Me3jHsCQtMO8YlqQFZhGQpAVmEZCkBWYRkKQFZhGQpAVmEZCkBWYRkKQFZhGQpAX2/wAwoia98Qil9wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_state_histogram(states_with_probs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It's a good time to experiment a little bit by yourself. Here are some ideas:\n",
    "\n",
    "- Run this code on different graphs. \n",
    "- Print out what's the distribution of the results. How \"convinced\" the algorithm is that given solution is correct?\n",
    "- Try running it several times on one graph. Do you get similar distributions every time?\n",
    "- How big graphs you can solve with reasonable calculation times?\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## QAOA - how does it work?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Black box approach"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The easiest way to look at QAOA is to think only about its inputs and outputs.\n",
    "<br>\n",
    "To solve the MaxCut problem the only thing we needed was the graph, which we passed to the instance of `maxcut_qaoa` class. But if you look into [its body](https://github.com/rigetticomputing/grove/blob/master/grove/pyqaoa/maxcut_qaoa.py), you will see that inside we use `qaoa` class with a little bit more advanced version of the following line:\n",
    "\n",
    "```\n",
    "QAOA_inst = QAOA(connection, qubits, steps=steps, cost_ham=cost_operators,\n",
    "                     ref_ham=driver_operators, \n",
    "                     init_betas=initial_beta,\n",
    "                     init_gammas=initial_gamma,\n",
    "                     minimizer_kwargs=minimizer_kwargs)\n",
    "```\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`connection` is just a regular `api.QVMConnection()` - so we tell QAOA how it should talk to the Forest cloud."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`qubits` are the same qubits that we used earlier."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`steps` is an integer, which specifies how long our quantum program will be. By default it's equal to 1 and that's the value we used in our MaxCut example. We will return to this later."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`cost_ham` (`cost_operators`) is how we encode our problem. This is where all the information about our graph sits and it's encoded as a combination of Pauli operators."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`ref_ham` (`driver_operators`) is how we encode what changes are possible in the realm of our problem. In `maxcut_qaoa` we say that we can flip every qubit from 0 to 1. Again - we use only Pauli operators to create it."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`init_gammas` and `init_betas` are just the initial values for the $\\gamma$ and $\\beta$ angles."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`minimizer_kwargs` is a dictionary with the minimizer arguments. As you can see in an example below, we specify the minimization method and its options. \n",
    "```\n",
    "minimizer_kwargs = {'method': 'Nelder-Mead',\n",
    "                            'options': {'ftol': 1.0e-2, 'xtol': 1.0e-2,\n",
    "                                        'disp': False}}\n",
    "\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You may wonder why it's a `minimizer` since we want to find a maximum cut. This is due to a nice property of optimization problems, that maximization is just a minimization with of the same problem but with a minus in the cost function. Since most optimization algorithms are stated in terms of minimization, we stick to this convention."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After we have all that we need to run the actual algorithm and find the right angles. Then we create a quantum program which gives us solution to the problem we specified with the `cost_operators`:\n",
    "```\n",
    "betas, gammas = QAOA_inst.get_angles()\n",
    "QAOA_inst.get_parametrized_program()\n",
    "```\n",
    "\n",
    "Or actually, we can go straight for the answer with:\n",
    "```\n",
    "QAOA_inst.get_string(betas, gammas)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is probably not the simplest black box interface, but honestly, it's not extremely complicated either. The only non-obvious part is what should the cost and driver operators be. Let's see in the next section!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### How to encode our problem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since in MaxCut we want to create to subsets of nodes, given node is either in set $A$ or $B$.\n",
    "<br>\n",
    "The total cost of our solution is: \n",
    "<br>\n",
    "<br>\n",
    "$C_{total} = \\sum C_{ij} = \\sum \\frac{1}{2} w_{ij} (1 - z_i z_j) $,\n",
    "<br>\n",
    "<br>\n",
    "where the sum goes through all the pairs of nodes and values $z_i$ are equal to 1 if given node is in the set $A$ and -1 if it's in the $B$. $w_{ij}$ is the weight of given edge, which in our current example is equal to 1.\n",
    "<br>\n",
    "\n",
    "It turns out that we can construct a matrix (Hamiltonian), which correspons to this cost function using only Pauli operators:\n",
    "<br>\n",
    "<br>\n",
    "$H_{cost} = \\sum \\frac{1}{2} w_{ij} (\\mathbb{1} - \\sigma^z_i \\sigma^z_j$)\n",
    "<br>\n",
    "\n",
    "If you wonder what $\\mathbb{1}$ is, it's an identity matrix, which can also be considered as special a Pauli matrix.\n",
    "\n",
    "How does it look like in the code?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pyquil.paulis import PauliTerm, PauliSum\n",
    "import networkx as nx\n",
    "\n",
    "maxcut_graph = nx.Graph()\n",
    "for edge in first_graph:\n",
    "    maxcut_graph.add_edge(*edge)\n",
    "graph = maxcut_graph.copy()\n",
    "        \n",
    "cost_operators = []\n",
    "for i, j in graph.edges():\n",
    "    cost_operators.append(PauliTerm(\"Z\", i, 0.5)*PauliTerm(\"Z\", j) + PauliTerm(\"I\", 0, -0.5))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The code above is copied from the `maxcut_qaoa` class. Let's see how some of these operators look:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(0.5+0j)*Z0*Z1 + (-0.5+0j)*I\n",
      "(0.5+0j)*Z0*Z2 + (-0.5+0j)*I\n"
     ]
    }
   ],
   "source": [
    "print(cost_operators[0])\n",
    "print(cost_operators[1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is how we encode our `cost_operators`, so now let's go to the `driver_operators`.\n",
    "They specify how our state can change. In this basic case, let's say that we allow every qubit to just change from 0 to 1 or otherwise. This means that we need a bunch of X gates (which correspond to $\\sigma_x$ Pauli matrices):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "driver_operators = []\n",
    "for i in graph.nodes():\n",
    "    driver_operators.append(PauliSum([PauliTerm(\"X\", i, -1.0)]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You don't need to bother with the `PauliSum` function, it's here just so that the input type is correct, but it doesn't play any other role."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(-1+0j)*X0\n"
     ]
    }
   ],
   "source": [
    "print(driver_operators[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And that's it!\n",
    "<br>\n",
    "Now let's try to actually run the QAOA by ourselves:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "from grove.pyqaoa.qaoa import QAOA\n",
    "\n",
    "initial_beta = 0\n",
    "initial_gamma = 0\n",
    "minimizer_kwargs = {'method': 'Nelder-Mead',\n",
    "                            'options': {'ftol': 1.0e-2, 'xtol': 1.0e-2,\n",
    "                                        'disp': False}}\n",
    "\n",
    "QAOA_inst = QAOA(qvm_connection, qubits, steps=1, \n",
    "                     cost_ham=cost_operators,\n",
    "                     ref_ham=driver_operators, \n",
    "                     init_betas=initial_beta,\n",
    "                     init_gammas=initial_gamma,\n",
    "                     minimizer_kwargs=minimizer_kwargs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WARNING: Fast method for expectation will be used. Noise\n",
      "                     models will be ineffective\n",
      "Values of betas: [0.39137304]\n",
      "Values of gammas: [0.95239311]\n",
      "And the most common measurement is... \n",
      "(0, 1, 1, 1)\n"
     ]
    }
   ],
   "source": [
    "betas, gammas = QAOA_inst.get_angles()\n",
    "print(\"Values of betas:\", betas)\n",
    "print(\"Values of gammas:\", gammas)\n",
    "print(\"And the most common measurement is... \")\n",
    "most_common_result, _ = QAOA_inst.get_string(betas, gammas)\n",
    "print(most_common_result)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Awesome! Now you have all the tools to use QAOA, except the most important one..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Knowing how it actually works ;)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### What are betas and gammas?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Imagine you're in the mountains and you want to find the lowest valley. In this analogy if we used the classical gradient descent method, you will start from some point and climb downwards until we can't find anything better. Most of the classical algorithms work this way - you are a moving around the mountain range to find this valley.\n",
    "<br>\n",
    "Here it's a little bit different. You start with some simple landscape and standing in the lowest point. Then, very, very slowly, the mountains start erecting from the ground. If this means that the minimum shifts, you do a little step. In the end, the mountains are the same as they were in the previous example and you are standing in the lowest valley.\n",
    "<br>\n",
    "What's happening here is the adiabatic transformation from some initial hamiltonian to the hamiltonian representing our problem. It is proven that if it's done slowly enough, the ground state (state of the minimal energy) of the initial hamiltonian will transform smoothly into the ground state of the final hamiltonian. It is also how the quantum annealers work (e.g. D-Wave).\n",
    "<br>\n",
    "One of the problems is that using Rigetti's quantum computer (or any other gate based machine) we cannot do this transformation in infinite number of infinitely small steps. We need to set how many steps we want and this is what `steps` parameter does. The more steps we have, the higher chance we will get the good result, but we also need more gates and running the algorithm gets longer and longer. Also, it means we have more betas and gammas to optimize. In this metaphore, betas and gammas tell you how exactly should you add new mountains in each step."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For those of you, who prefer equations over methaphores - below is the short version of how it works. We start with some definitions and then go to the explanations:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "${\\displaystyle U(C,\\gamma) = e^{-i \\gamma C} = \\prod_{\\alpha} e^{-i \\gamma C_{\\alpha}} }$, where $C_{\\alpha}$ are the cost operators we defined earlier.\n",
    "<br>\n",
    "${\\displaystyle U(B,\\beta) = e^{-i \\beta B} = \\prod_{j} e^{-i \\beta B_{j}} }$, where $B_{j}$ are the driver operators we defined earlier.\n",
    "<br>\n",
    "$\\ket{s}$ is the initial state.\n",
    "<br>\n",
    "And finally: $\\ket{\\beta, \\gamma} = U(B, \\beta) U(C, \\gamma) \\ket{s}$\n",
    "<br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is for the case where `steps=2`, but it generalizes nicely to higher number of steps, e.g.:\n",
    "<br>\n",
    "$\\ket{\\beta, \\gamma} = U(B, \\beta_1) U(B, \\beta_0) U(C, \\gamma_1) U(C, \\gamma_0)\\ket{s}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you are not familiar with this notation here is a step-by-step description of what does it mean in practice:\n",
    "- Let's take some initial state $\\ket{s}$\n",
    "- Evolve this state using the operator $U(C, \\gamma_0)$ \n",
    "- Evolve this state further using the rest of operators, from right to left.\n",
    "\n",
    "\"Evolving\" is just a fancy name for \"changing\" that we use when talking about the states in quantum computing."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And now, if the values of angles $\\beta$ and $\\gamma$ are the right ones, we get a state which is a ground state of the Hamiltonian defining our problem, which should correspond to our solution.  \n",
    "If you wonder why we call $\\gamma$ and $\\beta$ angles, that's because $U(C, \\gamma)$ and $U(B, \\beta)$ correspond to rotations by angles $\\gamma$ and $\\beta$ in some high dimensional space."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### How to find the right angles? The hybrid part"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "QAOA is an iterative algorithm, which in every iteration uses a quantum and a classical part."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Quantum**\n",
    "1. Prepare a quantum state based on the current values of $\\gamma$ and $\\beta$.\n",
    "2. Measure the state.\n",
    "3. Repeat many times to get some statistics.\n",
    "\n",
    "**Classical**\n",
    "\n",
    "Here we use some classical optimization method, e.g. Nelder-Mead algorithm, and in each iteration of this algorithm we perform the quantum operations to evaluate how good our current angles are.\n",
    "<br>\n",
    "The algorithm ends when we can't further improve the angles."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Summary and next steps"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What we learned?\n",
    "\n",
    "- how to solve MaxCut problem using QAOA\n",
    "- how to use QAOA class from grove\n",
    "- how to encode an optimization problem for QAOA\n",
    "- what the QAOA is and how it works\n",
    "\n",
    "What's next?\n",
    "\n",
    "- Solving TSP using QAOA"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Additional resources"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- [Pauli matrices](https://en.wikipedia.org/wiki/Pauli_matrices)\n",
    "- [Original QAOA paper](https://arxiv.org/abs/1411.4028)\n",
    "- [QAOA Grove docs](https://grove-docs.readthedocs.io/en/latest/qaoa.html)\n",
    "- [Source code of MaxCut QAOA in Grove](https://github.com/rigetticomputing/grove/blob/master/grove/pyqaoa/maxcut_qaoa.py)\n",
    "- [Presentation from one of the creators of QAOA](https://youtu.be/J8y0VhnISi8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}