{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"tags": [
"remove_cell"
]
},
"source": [
"# Solving combinatorial optimization problems using QAOA"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"> ℹ️ Originally from [the Qiskit textbook](https://qiskit.org/textbook/ch-applications/qaoa.html), this notebook demonstrates how to take existing Qiskit work and run it on a Rigetti backend.\n",
"\n",
"In this tutorial, we introduce combinatorial optimization problems, explain approximate optimization algorithms, explain how the Quantum Approximate Optimization Algorithm (QAOA) works and present the implementation of an example that can be run on a simulator or on a real quantum system."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import networkx as nx\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Combinatorial Optimization Problem\n",
"\n",
"Combinatorial optimization problems involve finding an optimal object out of a finite set of objects. We would focus on problems that involve \n",
"finding \"optimal\" bitstrings composed of 0's and 1's among a finite set of bitstrings. One such problem corresponding to a graph is the Max-Cut problem. \n",
"\n",
"### Max-Cut problem\n",
"A Max-Cut problem involves partitioning nodes of a graph into two sets, such that the number of edges between the sets is maximum. The example below\n",
"has a graph with four nodes and some of the ways in which it can be partitioned into two sets, \"red\" and \"blue\" is shown.\n",
"\n",
"![](https://qiskit.org/textbook/ch-applications/images/qaoa_maxcut.svg)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For 4 nodes, as each node can be assigned to either the \"red\" or \"blue\" sets, there are $2^4=16$ possible assigments.\n",
"\n",
"Out of which we have to find one that gives maximum number of edges between the sets \"red\" and \"blue\". The number of such edges between two sets in the figure, as we go from left to right, are 0, 2, 2, and 4. We can see, after enumerating all possible $2^4=16$ assignments, that the rightmost figure is the assignment that gives the maximum number of edges between the two sets. Hence if we encode \"red\" as 0 and \"blue\" as 1, the bitstrings \"0101\" and \"1010\" that represent the assignment of nodes to either set are the solutions. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As you may have realized, as the number of nodes in the graph increases, the number of possible assignments \n",
"that you have to examine to find the solution increases exponentially."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## QAOA\n",
"\n",
"QAOA (Quantum Approximate Optimization Algorithm) introduced by Farhi et al.[1] is a quantum algorithm that attempts to \n",
"solve such combinatorial problems. \n",
"\n",
"It is a variational algorithm that uses a unitary $U(\\boldsymbol{\\beta}, \\boldsymbol{\\gamma})$ characteized by the parameters \n",
"$(\\boldsymbol{\\beta}, \\boldsymbol{\\gamma})$ to prepare a quantum state $\\lvert \\psi(\\boldsymbol{\\beta}, \\boldsymbol{\\gamma}) \\rangle$. \n",
"The goal of the algorithm\n",
"is to find optimal parameters $(\\boldsymbol{\\beta}_{opt}, \\boldsymbol{\\gamma}_{opt})$ such that the \n",
"quantum state $\\lvert \\psi(\\boldsymbol{\\beta}_{opt}, \\boldsymbol{\\gamma}_{opt}) \\rangle$ encodes the solution to the problem. \n",
"\n",
"The unitary $U(\\boldsymbol{\\beta}, \\boldsymbol{\\gamma})$ has a specific form and is composed of two unitaries \n",
"$U(\\boldsymbol{\\beta}) = e^{-i \\boldsymbol{\\beta} H_B}$ and $U(\\boldsymbol{\\gamma}) = e^{-i \\boldsymbol{\\gamma} H_P}$\n",
"where $H_B$ is the mixing Hamiltonian and $H_P$ is the problem Hamiltonian. Such a choice of unitary drives its inspiration \n",
"from a related scheme called quantum annealing.\n",
"\n",
"The state is prepared by applying these unitaries as alternating blocks of the two unitaries applied $p$ times such that \n",
"\n",
"$$\\lvert \\psi(\\boldsymbol{\\beta}, \\boldsymbol{\\gamma}) \\rangle = \\underbrace{U(\\boldsymbol{\\beta}) U(\\boldsymbol{\\gamma}) \n",
" \\cdots U(\\boldsymbol{\\beta}) U(\\boldsymbol{\\gamma})}_{p \\; \\text{times}} \n",
"\\lvert \\psi_0 \\rangle$$\n",
"\n",
"where $\\lvert \\psi_0 \\rangle$ is a suitable initial state. \n",
"\n",
"We will demonstrate these steps using the Max-Cut problem discussed above. For that we would first define\n",
"the underlying graph of the problem shown above."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAb4AAAEuCAYAAADx63eqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAA20UlEQVR4nO3deVyVdd4//tdZgAOHRZFFXAgQNxCTxQPZTHPXlN3ZXppLrmX7ZrdZ9qu5p7mbKVwzzcqlcElNYZqZyuo31nSPdwscQVTcUhRERTggysEDh7Nc1/cP8xqZXNg/1znn9Xw85h8DfFHz4MXnc30+70sjy7IMIiIiH6EVHYCIiKg7sfiIiMinsPiIiMinsPiIiMinsPiIiMinsPiIiMinsPiIiMin6EUHIM/glmRUW+1wuiX46bSIDjVAp9WIjkVE1GYsProsq92JbfursbXkFMpqbIBGhgYayJAhy0BCZDBuT4nB6KRohBj8RMclImoVDSe30L9zuiVsNFdg7Q/lcEsy9DoNDHodtBet8CRJht3lhsstQ6fVYNqoOEwyxcJPx91zIlI3Fh+1YGmwY07eHpTVnoPRXw99K4rM5ZZgc7gQHxGMBWOHIyrE0A1JiYjah8VHCkuDHY+tL0KdzYEQgx4aTeuf4cmyjAa7C+FGf6yYks7yIyLV4r4UATi/vTknbw/qbA6EBvq1qfQAQKPRIDTQD3U2B17M2wOXW+qipEREHcPiIwDARnMFymptCDF07LxTiEGPo7XnsNFc0UnJiIg6F7c6CVa7E/cs/x4Beu1ln+k5m2w4/NUanCnfB7/AYMTdcD+ikjIv+bEut4Rml4S/PXU9T3sSkepwxUfYtr/659Obl/+/w5FtH0Gj0yHrqcUYfMcjKN22Hrbaykt+rF6nhVuSsW1/dVdFJiJqNxYfYWvJKeh1l3+m53Y0o/bwTsT96l7o/A0I6zcQvRJHwLLvh8t+jl6nwdaSU10Rl4ioQ1h8Ps4tySirscGg1132Y5rOVEOj0SIwPFr5M2NkPzReZsUHAAa9DkdrbHBL3EknInVh8fm4aqsd0MgtLqf/O7fDDl1AIBoaGnDmzBk4nQ7oAoLgdtgv+zla7fkJL9XWy38MEZEIHFnm45xuCRpc+eqCzt8AR+M5nDp1CkajETU1NWg4Xg63y42GhgYEBQVBp/vlilGr0cDJaw1EpDIsPh/np9NCxpW3IwPCIuFwNKOX0Q+9+vYFIGPfgX9A7tUHZ86cQWXlSQQEBMBoDIbRaERgYCA0Gg0kWeYIMyJSHRafj4sONQCyBpJ0+e3OM9YGhMWnoKbo7+gROR02y3HUl+/FtQ/+fzBG9IEsS2hsbILNdg5VVVVwOp0IDAqCX2Aw7GeqIYdd0+YL8UREXYX3+AgPrd2BY6dtCPL/5e9Bzc3NOHbsGGL7ROPo39fjzLH98DMYEfebsZe/x+dy4UyDDX72swgp3gCNRoPMzExkZWXBZDKhR48eXfwdERFdHouP8OedJ7D0m8MIDfz3y+Yyjh2rQEhICMLDw9v0Na1NTjz324G4N7Uvjh07hoKCAhQUFKCoqAj9+/dHZmYmMjMzMWLECPj7+3feN0NEdBUsPkKD3Ym7LzG5xWq1ora2FvHx8W3aqrzS5BaXy4WSkhKlCI8cOYJrr71WKcLExERuixJRl2LxEQBg7Y/l+OC7owg1nB9QLUluHDlyFP369UVgYFCrv44sy7DanZj5qwRMvS7uqh9vtVpRWFioFGFjYyOysrKUIoyIiOjAd0VE9EssPgJw/lrDzHWFKK+1ITTQD9XV1ZAkN2Ji+rTp61ibnIiPMGL11IxWvcvv3508eRIFBQXIz89HYWEhIiMjleeDqampCAwMbPPXJCK6GIuPFBfex1d91obaqhNIHJAAna51B3+74n18kiRh//79ymrw4MGDGDp0qLIiHDJkCLRaXpcgorZh8VEL1fVNuO0PG6EJjUZUz9A2vYE9ISIY87vwDeyNjY3YuXMn8vPzUVBQgLq6OowcOVIpwpiYmC75e4nIu7D4qIVPP/0UeZ/8Bb999L+x9sdjP7+1QQODXtfinp8kybC73HC5Zeh0Gky/Lg6TTLHt2t5sL4vFomyLms1mhISEKM8GMzIyEBwc3G1ZiMhzsPhIUV9fj3HjxmHZsmUYPHgwGuxObNtfja0lp3C0xgYZMrQ/T2TRQIOESCNuT4nBLUnRwt+7J0kSSktLkZ+fj/z8fOzduxeJiYnKanDYsGGXHKtGRL6HxUeKN954A35+fpgzZ84v/plbOj9w2umW4KfTIjrUAN0VBluL1tzcjF27dikrwsrKSqSnpysHZfr3789rE0Q+isVHAICSkhLMmTMHeXl5XrlFWFdXB7PZrBShXq9XtkVNJhPCwsJERySibsLiI0iShClTpmDKlCn4z//8T9FxupwsyygrK1NOi+7cuRNxcXFKEQ4fPpzTZIi8GIuPsGXLFvzjH//Ae++955Pbf06nEyUlJcpp0bKyMowYMUJ5PpiQkOCT/16IvBWLz8edPn0a48ePx6pVqxAfHy86jipcmCZz4aCMw+GAyWRSirBXr16iIxJRB7D4fNzvfvc7REVF4ZlnnhEdRbVOnDihrAYLCwsRHR3dYpqMwdA19xaJqGuw+HxYUVERfv/73yM3N5ejwFrJ7XYr02Ty8/Nx6NAhJCcnK88HBw8ezGkyRCrH4vNRTqcTkyZNwlNPPYX/+I//EB3HYzU2NipDtvPz81FfXw+TyaQUYe/evUVHJKJ/w+LzUWvXrsXOnTuxZMkSHtzoRNXV1S2myYSFhSnbounp6TAajaIjEvk8Fp8PqqqqwoMPPoh169ahb9++ouN4LUmScOjQIeXaRElJCQYPHqysBpOTkzlNhkgAFp8PmjNnDgYPHoyZM2eKjuJT7HY7iouLlSKsqqpCRkaGUoT9+vXj6puoG7D4fMz333+PhQsXYvPmzbykLdjp06eVEiwoKIC/v7+yLTpy5EiEhoaKjkjklVh8PqS5uRkPPPAAXn75ZWRlZYmOQxe5ME3mwt3BXbt2IS4uDllZWcjKykJKSgr8/MQOAifyFiw+H7JixQqUlZUhOztbdBS6CofD0WKaTHl5OVJTU5UVYXx8PLdFidqJxecjKioqMGPGDGzatAlRUVGi41Ab1dfXY8eOHUoRulwu5dlgZmYmwsPDRUck8hgsPh8gyzKeffZZZGZmYvLkyaLjUAfJsozjx48rzwYLCwvRp08fpQRTU1MREBAgOiaRarH4fMA333yDlStXYsOGDdDr9aLjUCdzu93Yt2+fsho8fPgwhg0bpswWHThwIKfJEF2ExeflGhsbMXbsWPzpT39Camqq6DjUDc6dO4eioiLlIn1DQ4MyTSYrK4tb3eTzWHxe7u2338aZM2fw2muviY5Cgpw6dUrZFjWbzejZs2eLaTJBQUGiIxJ1KxafFzty5Agef/xxbN68mYcfCMD5aTIHDx5UinDfvn0YMmSI8nwwKSmJ02TI67H4vJQsy3j00Udx6623YuzYsaLjkEo1NTWhuLhYeT5YU1OjTJPJysriSDvySiw+L7V161Zs3rwZa9as4cEGarXa2lqYzWalCA0Gg1KCGRkZnCZDXoHF54WsVivGjRuHt956C0lJSaLjkIeSZRlHjhxRDsns3r0bCQkJyrYop8mQp2LxeaF58+ZBlmXMnTtXdBTyIg6HA7t371aK8Pjx40hLS1OKMC4ujtNkyCOw+LzM/v378fzzzyM3N5fbUtSlzp49C7PZrBShLMvK3UGTyYSePXuKjkh0SSw+LyJJEqZPn44HHngAd9xxh+g45ENkWUZFRYVSgkVFRejbt68yZHvEiBF8GwipBovPi/z5z3/GV199hZUrV3LLiYRyuVzYu3evcm2itLQUKSkpykGZxMREHroiYVh8XqKurg7jx4/H+++/jwEDBoiOQ9RCQ0MDCgsLlSK02WzKNJnMzExOk6FuxeLzEq+99hp69OiBWbNmiY5CdFWVlZXKlYkdO3YgIiJCKcG0tDROk6EuxeLzAsXFxXjllVeQl5fHHxjkcS5Mk7nwEt4DBw5g6NChykGZoUOHcluUOhWLz8O5XC48+OCDePTRR/Hb3/5WdByiDmtsbERxcbFyUKa2thYjR45Ung/26dNHdETycCw+D/fRRx8hPz8fy5Yt44EW8koWi0W5NlFQUICgoCBlNZiRkYGQkBDREcnDsPg8mMViwcSJE5GTk4PY2FjRcYi6nCRJKC0tVUpwz549GDBggLIaHDZsGN85SVfF4vNgc+fORXx8PB577DHRUYiEcDgc2LVrl3JQ5sSJE0hPT1eKMDY2ljsh9AssPg+Vn5+PN998E1u2bEFAQIDoOESqUFdXhx07dijPBzUajVKCJpMJPXr0EB2RVIDF54EcDgfGjx+PF154Addff73oOESqJMsyjh07pqwGd+7cif79+yvXJjhNxnex+DzQ6tWr8dNPP2HBggWioxB5DKfTib179ypFePToUVx77bVKESYmJnJb1Eew+DzMyZMnMXXqVGzYsAG9e/cWHYfIY1mtVhQWFipF2NTU1GLIdmRkpOiI1EVYfB5ElmXMmjULaWlpmDZtmug4RF7l5MmTyrPBHTt2IDIyUhmynZqaisDAQNERqZOw+DzI//7v/2L58uXYuHEjXwBK1IUkScL+/fuVaxMHDhxAUlKSsiIcMmQIp8l4MBafh2hqasK4cePwhz/8Aenp6aLjEPmUxsZGFBUVKUVYV1eHkSNHKkUYExMjOiK1AYvPQ7zzzjuorq7G66+/LjoKkc+zWCzKs0Gz2YyQkBDlkExGRgaCg4NFR6QrYPF5gLKyMjzyyCPYvHkzevXqJToOEV1EkiQcPnxYeT64d+9eJCYmKqvB5ORkTpNRGRafysmyjCeeeAI33ngjxo8fLzoOEV1Fc3Nzi2kylZWVLabJ9O/fn9cmBGPxqdxXX32F9evXY926ddDpdKLjEFEb1dXVwWw2K0Wo1+uVbVGTyYSwsDDREX0Oi0/Fzp07h7Fjx2LBggVISUkRHYeIOkiWZZSVlSnbosXFxbjmmmuUbdHhw4dzmkw3YPGp2MKFC9Hc3IxXXnlFdBQi6gJOpxMlJSXKS3jLy8sxYsQIpQgTEhK4LdoFWHwq9dNPP+GZZ55Bbm4ut0KIfITVam0xZNvhcMBkMilFyMNtnYPFp0KSJOGhhx7Cvffei7vvvlt0HCISQJZlnDhxQrk7WFhYiOjoaOWQTGpqKgwGg+iYHonFp0J//etf8emnn2L16tWcDkFEAAC32419+/YpRXjo0CEkJycrB2UGDx7MnxetxOJTmbNnz+KBBx7AO++8g0GDBomOQ0QqZbPZUFRUpJwWra+vh8lkUoqQQ+wvj8WnMq+//jqCgoIwe/Zs0VGIyINUVVUpq0Gz2YzQ0FBlyHZ6ejqMRqPoiKrB4lORPXv24KWXXkJeXh7/T0pE7SZJEg4dOtRimsygQYNaTJNR071gtySj2mqH0y3BT6dFdKgBOm3XnWZl8amE2+3GlClTMH36dIwePVp0HCLyIna7HcXFxUoRVldXIz09XSnCfv36dfu1CavdiW37q7G15BTKamyARoYGGsiQIctAQmQwbk+JweikaIQYOvdtNCw+ldi0aRO2b9+Od999l/d2iKhLnT59WtkWzc/PR0BAgHJadOTIkQgNDe2yv9vplrDRXIG1P5TDLcnQ6zQw6HXQXrTCkyQZdpcbLrcMnVaDaaPiMMkUCz9d5xzeYfGpQG1tLSZMmIDVq1cjLi5OdBwi8iGyLOPo0aNKCe7atQtxcXEtpsl01vs/LQ12zMnbg7LaczD666FvRZG53BJsDhfiI4KxYOxwRIV0/AoHi08FXnnlFfTt2xdPPvmk6ChE5OMcDgf27NmjrAjLy8uRmpqqrAjj4+PbtStlabDjsfVFqLM5EGLQt+lryLKMBrsL4UZ/rJiS3uHyY/EJZjab8frrryM3N5eXUYlIderr62E2m5UidLlcypWJzMxMhIeHX/VrON0SZq4rRHmtDaGB7V89WpuciI8wYvXUjFatFi+HxSeQ0+nEhAkT8Nxzz+GGG24QHYeI6IpkWcbx48eVu4NFRUWIiYlRtkVTU1MREBDwi89b+2M5PviuDKGXWen99PlqnD22H26nA/7GMPQz3Yre1/7yZ6Isy7DanZj5qwRMvS6u3d8Hi0+gnJwclJSUYPHixaKjEBG12YVpMheGbB8+fBgpKSlKEQ4cOBDnHG7cs/x7BOi1l12l2WorEdgjElq9HxpPn8Kejxcg+f5nEdI77hcf63JLaHZJ+NtT17f7tCeLT5BTp05h8uTJWL9+Pfr06SM6DhFRh507d67FNJmGhgb0yhiD0oCBiAgztuqQTGNdFUo2zUfCbycicsjIS36MtcmJ5347EPel9WtXThafILNnz0ZycjIeeugh0VGIiLrEqVOnMHNdISrONKHZ1gC9Xg+j0YhgoxFBRmOL2aKlf/8I1Xu/h+Rywhgdi2snvgid/6XPPTQ6XIjrZcQH0y5djFejb9dnUYf83//9H8rKyvDmm2+KjkJE1GWionujUReC2D49oNHEwG63w2az4XRdHU6cPAmDwYDgYCOMxmAMuOVBDLh5EqyVR1B//CdodJevJ4Neh6M1NrgluV0TXjjKu5vZ7XYsWLAAc+fO5ZuWicirVVvtgEaGVquBRqNBYGAgIiIicM0112DQoEGIiIiA2y2hsvIkysvLodFqEdZvIJqtdTi1638v+3W12vMTXqqt9nbl4oqvm+Xk5GDYsGEwmUyioxARdSmnW4IGl16RabVaBAcHIzg4GLIsQ5Yl5Z/JsgT72Zorfm2tRgOnW7rix1z2c9v1WdQuFRUVyMvLw6xZs0RHISLqcn46LWRc+RhJQ50Fp/b+gF49wiBLEs6U7UXNATN6XDP0ip8nyXK7R5hxxddNZFlGdnY2Hn74YURFRYmOQ0TU5aJDDYCsgSTJLWZxXiDLMiyWGrhO7kPhym8BWUZAaC8k3DQBvRJHXPbrStL5gdbRoe0b+sHi6yZff/01zpw5g/Hjx4uOQkTULXRaDeIjjTh22oYg/1/WTX19PbQBQciY+mqbRpjZXW4kRBrb/eoibnV2g8bGRrz11luYO3euqt6BRUTU1W5PiYHL/cvtTrfbDYvFgpiYmDbP/nS5ZdyeEtPuTCy+brBixQpkZWXh2muvFR2FiKhbjU6Khk6rgevfDqJYLBaEhoa2eUaxyy1Bp9XglqTodmdi8XWxw4cP44svvsAzzzwjOgoRUbcLMfhh2qg42BwuXJiX0tTUhIaGBkRGRrbpa8myDJvDhemj4jr0cloWXxeSJAnZ2dl44okn0LNnT9FxiIiEmGSKRXxEMBrs58uvqqoKUVFRbX7002B3ISEiGJNMsR3Kw+LrQlu3boXL5cI999wjOgoRkTB+Oi0WjB2OcKM/Tp0+C2g06NEjrNWfL8syrE1OhBv9MX/s8A69kghg8XUZq9WKZcuWYe7cuS3m0RER+aKoEAPm35mIpupjCOkZeckDL5fickuw2s+/h68zXkIL8DpDl3nnnXdw8803Y+jQK1/CJCLyFZtyVmBGYjCiRw3Gmh/K0ehwQ6/TwKDXtbjnJ0ky7C43XG4ZOp0GM3+VgEmm2A6v9C5g8XWBffv2Yfv27cjNzRUdhYhIFfbs2YMff/wReXl5MBqNuDe1L7btr8bWklM4WmODDBlajQaSfP5yekKkEbenxOCWpOgOHWS5FL6WqJNJkoSpU6di0qRJGDNmjOg4RETCud1uTJkyBdOmTcOtt976y38unR847XRL8NNpER1qaPfl9Nbgiq+T5eXlISgoCLfddpvoKEREqpCXl4fQ0FCMHj36kv9cp9WgT4/AbsvD4utEdXV1WLVqFVasWNHmSQRERN7o9OnTWLVqFVauXKman4s8btiJlixZgjvvvBMJCQmioxARqcLbb7+Nu+66S1U/F7ni6yQ7d+5EUVERD7QQEf1MrT8XueLrBC6XC9nZ2Zg9ezaCgoJExyEiEu7Cz8X/+q//Ut3PRRZfJ9i4cSN69+6NG2+8UXQUIiJV2Lx5M6KionDTTTeJjvIL3OrsoOrqaqxduxZr165VzYNbIiKRLBYLPvzwQ+Tk5Kjy5yJXfB20aNEijB8/Hv369RMdhYhIFZYsWYL7778fsbEdGybdVVh8HfDDDz/g0KFDmD59uugoRESqYDabsXfvXjz00EOio1wWi6+dHA4H5s+fjxdffBH+/v6i4xARCed0OjF//nzMnj27zS+Y7U4svnZas2YNBg0ahFGjRomOQkSkChs2bEC/fv1www03iI5yRTzc0g4nTpzA5s2bsXHjRtFRiIhUoaqqCuvXr8e6detUeaDlYlzxtZEsy5g/fz6mTZuG6Oho0XGIiFRh8eLFmDBhAvr27Ss6ylWx+Nro22+/RVVVFSZNmiQ6ChGRKlw46Ddt2jTRUVqFxdcGjY2NWLRoEebOnQu9nrvEREQOhwMLFizwqIN+LL42WL16NdLT05GWliY6ChGRKqxbtw6JiYkeddCPy5ZWOnr0KD777DNs3rxZdBQiIlWorKzEpk2b8NFHH4mO0iZc8bWCLMvIzs7GI488gvDwcNFxiIhUYeHChZg8eTJiYmJER2kTFl8rfPnll2hsbMTYsWNFRyEiUoXt27fj2LFjmDx5sugobcatzqtoaGjA0qVLsWjRImi1/D2BiMhut2PhwoV49dVX4efnJzpOm/En+VW89957uOGGG5CcnCw6ChGRKuTk5CA5ORkmk0l0lHbhiu8KDhw4gK+//hp5eXmioxARqUJFRQXy8vKwadMm0VHajSu+y5AkCdnZ2XjmmWcQGhoqOg4RkXCyLGPBggWYPn06oqKiRMdpNxbfZfz1r3+FXq/H7bffLjoKEZEqfPvtt6iursbEiRNFR+kQbnVewpkzZ/Dee+/h3Xff5YEWIiIATU1NWLRoEV5//XWPn1zFn+qXsGzZMowZMwYDBw4UHYWISBVWr16NtLQ0r5hc5dm13QV2796N/Px8HmghIvpZWVkZ/va3v3nN5Cqu+C7idrvx5ptv4vnnn0dQUJDoOEREwl14FdvMmTPRq1cv0XE6BYvvIps3b0avXr1w8803i45CRKQK27ZtQ319PcaNGyc6SqfhVufPLBYLPvzwQ3z44Yeqf3swEVF3sNlsWLJkCbKzs6HT6UTH6TRc8f1syZIluP/++xEbGys6ChGRKqxcuRJZWVkYPny46CidisUHwGw2Y+/evZgxY4boKEREqlBaWoovvvgCTz/9tOgonc7ni8/hcCA7Oxtz5syBwWAQHYeISDhZljFv3jw89thjXvkqNp8vvo8++gjx8fH49a9/LToKEZEqfPnll2hubsZ9990nOkqX8OnDLZWVldiwYYPHvT2YiKirXHgV2+LFi712cpV3flet5KlvDyYi6ioXXsWWlJQkOkqX8dni2759OyoqKjzy7cFERF3hp59+wtdff42nnnpKdJQu5ZPFZ7fbsWDBArz00kse+fZgIqLOduFVbE899RTCwsJEx+lSPll8H374IYYPH46RI0eKjkJEpAqfffYZAODOO+8UnKTr+dzhlvLycnzyySf4+OOPRUchIlKF+vp6LF++HMuWLfPaAy0X8/7v8CIX7qY8/PDDiIiIEB2HiEgVli9fjptvvhmDBw8WHaVb+NSK78Kw1QceeEB0FCIiVdi/fz/++c9/+tSr2Hym+Gw2G9566y3MmzfPq4atEhG1lyRJePPNN/Hss88iJCREdJxu4zNbnStWrMCoUaO8btgqEVF7/eUvf0FAQADGjBkjOkq38okV36FDh/DVV19hy5YtoqMQEanCmTNn8P777+O9997zuVexef2K78LdlCeffBI9evQQHYeISBWWLVuG2267DYmJiaKjdDuvX/F9/vnnkGUZd911l+goRESqsGfPHvz4448+daDlYl5dfPX19XjnnXd85m4KEdHVuN1uZGdnY9asWTAajaLjCOHVbfDOO+9g9OjRPnM3hYjoanJzcxEaGorRo0eLjiKM1674SkpK8N133yE3N1d0FCIiVTh9+jRWr16NVatW+dyBlot55YrvwlL+ueeeQ3BwsOg4RESq8Pbbb+Ouu+5CfHy86ChCeeWKLy8vDyEhIbj11ltFRyEiUoWdO3eiqKiIu2DwwuI7ffo0Vq1a5fNLeSKiC1wuF7KzszF79mwEBQWJjiOc1211LlmyBPfcc4/PL+WJiC74+OOPERUVhRtvvFF0FFXwqhVfYWEhiouLuZQnIvqZxWJBTk4OcnJyuAv2M69Z8TmdTsybNw8vvPACAgMDRcchIlKFt956C2PHjkVsbKzoKKrhNcW3ceNG9O3bF7/5zW9ERyEiUgWz2Yx9+/ZhxowZoqOoildsdVZVVWHdunVYt24dl/JERAAcDoeyC2YwGETHURWvWPEtWrQIEydORN++fUVHISJShQ0bNiA2NhY33HCD6Ciq4/HF99133+HIkSOYOnWq6ChERKpQVVWFjz76CHPmzBEdRZU8uviam5sxf/58vPjii/D39xcdh4hIFRYtWoQJEyagT58+oqOokkcX35o1a5CUlISsrCzRUYiIVOGHH35AaWkppk2bJjqKanns4ZaKigps2bIFmzZtEh2FiEgVHA4Hd8FawSNXfLIsY8GCBZgxYwaioqJExyEiUoW1a9di4MCBGDVqlOgoquaRK75//OMfsFgsmDBhgugoRESqcPLkSXz88cfYsGGD6Ciq53ErvsbGRixatAgvv/wy9HqP7G0iok63cOFCTJkyBb179xYdRfU8rvhWrVoFk8mEESNGiI5CRKQK27dvR0VFBR588EHRUTyCRy2Zjhw5gs8//xybN28WHYWISBXsdjsWLlyIV199FX5+fqLjeASPWfHJsozs7Gw8/vjjCA8PFx2HiEgVcnJyMGzYMJhMJtFRPIbHrPi++OILNDc349577xUdhYhIFSoqKpCXl8drXW3kEcVntVqxdOlSvPXWW9BqPWaRSkTUZWRZxvz583mtqx08okXeffdd3HjjjUhKShIdhYhIFb799lte62on1a/49u/fj2+//ZZvVSci+tmFa12vv/46r3W1g6pXfJIkITs7G88++yxCQ0NFxyEiUoXVq1cjPT0daWlpoqN4JFX/qvDJJ58gICAAY8aMER2FiEgVjh49ik8//ZTXujpAtcVXV1eHFStW4P333+db1YmI8K8DLTNnzkSvXr1Ex/FYqt3qXLp0Ke644w4MGDBAdBQiIlX4+9//DqvVinHjxomO4tFUueIrLi6G2WxGXl6e6ChERKpgs9mwZMkSzJs3DzqdTnQcj6a6FZ/L5UJ2djZmz56NoKAg0XGIiFRhxYoVuO666zB8+HDRUTye6lZ8H3/8MaKionDTTTeJjkJEpAqlpaX48ssvsWXLFtFRvIKqis9isSAnJwc5OTk80EJEhJZzinv27Ck6jldQ1Vbn4sWL8cADDyA2NlZ0FCIiVfjiiy/gcDg4p7gTqab48vPzceDAAUyfPl10FCIiVWhoaMDSpUsxd+5czinuRKr4N+lwODBv3jy8+OKLCAgIEB2HiEgV3nvvPfzmN7/hnOJO1i3P+NySjGqrHU63BD+dFtGhBui0/3qGt27dOiQmJuL666/vjjhERKp38OBBfP3117zW1QW6rPisdie27a/G1pJTKKuxARoZGmggQ4YsAwmRwbg9JQYpPSVs2rQJGzZs6KooREQe5cKc4qeeeopzirtApxef0y1ho7kCa38oh1uSoddpEOSvg/aiFZ4kyTh22oal3xzG6doa3HDHI+gVyfdJEREBwKeffgqtVos777xTdBSvpJFlWe6sL2ZpsGNO3h6U1Z6D0V8Pve7KjxAbGhpQbalBREw/JEQGY8HY4YgKMXRWHCIij1NfX49x48Zh2bJlGDx4sOg4XqnTDrdYGux4bH0RymttCDX4XbX0JElCVVUV+sREIyzQD+W1Njy2vgiWBntnRSIi8jjLly/HzTffzNLrQp1SfE63hDl5e1BncyA00K9Vl89ra2sRFBSEoCAjNBoNQgP9UGdz4MW8PXC5pc6IRUTkUfbt24ft27fjiSeeEB3Fq3VK8W00V6Cs1oYQQ+seGTY3N+Ps2bOIjo5u8echBj2O1p7DRnNFZ8QiIvIYF794OyQkRHQcr9bhwy1WuxNrfyiH0V93yZVe5c5/oHrv97DVnETkUBMGj5mBqqoqREREQK9v+ddrNBoY/fVY80M57k3tixCDX0fjERF5hE8++QQGgwG33Xab6Cher8Mrvm37q38+vXnpL+VvDEP/6+5A75Tzd/Tq661wu92XnTmn12nhlmRs21/d0WhERB7hwou3X3rpJc4p7gYdLr6tJaeg113+P1TE4HREDEyF3hAMWZJgsVQjJibmiv9x9ToNtpac6mg0IiKPsGzZMowZMwaJiYmio/iEDm11uiUZZTU2BPm37qWItsZG9AkOQWBg4BU/zqDX4WiNDW5JbjHhhYjI2+zevRsFBQXIzc0VHcVndGjFV221Axq5xeX0y3G6nGhubkZU1NUvqmu15ye8VFt5tYGIvJfb7UZ2djZmzZoFo9EoOo7P6FDxOd0SNGjdiszlckGWJFRXV6O+vh4ul+vKwTQaOHmtgYi82JYtW9CjRw/ccsstoqP4lA5tdfrptJDRusEvgQYDeoaHIzDQgAarFVVVVfDz80NwsBFGYzCCggKh0fyrhyVZht9VLsETEXmq2tpafPDBB1i1ahUPtHSzDhVfdKgBkDWQpMtvd8qSG7IkQZbPP68LCwlBj7AwQKNFU1MTbDYbampq0NxsR2BgEIxGI4KMQdBAf/7rExF5obfffht333034uPjRUfxOR2e1fnQ2h04dtqGIP9Ld+ix7z9Fxfeftviz2OvvwjXX39Xiz9xuNxobG3Hu3Dk0NDZDbz+D+yNrkJWVhczMTPTq1asjMYmIVKOoqAi///3vkZube9XDftT5Olx8f955Aku/OYzQwM67bG5tcmJKWi/0PHsY+fn5KCwsRO/evZGVlYWsrCyMGDGCL6wlIo/kcrkwadIkPP7447jppptEx/FJHS6+BrsTdy//HgF67VUHU7eGyy2h2SXhb09dr0xucbvd2LdvH/Lz81FQUIDDhw8jJSVFWQ0OHDiQe+RE5BHWr18Ps9mMpUuX8ueWIJ3yWqK1P5bjg++OItTQugHVlyPLMqx2J2b+KgFTr4u77MedO3cOhYWFKCgoQH5+Pmw2m1KCmZmZiIiIaHcGIqKuYrFYMHHiROTk5CA2NlZ0HJ/VKcXndEuYua7w/CuJOrDlaW1yIj7CiNVTM9q0eqysrER+fr6yLRoVFYXMzExkZWUhNTUVBgMPyRCReC+//DKuueYaPP7446Kj+LROexHthffx1dkcCDHo27Tyk2UZDXYXwo3+WDElvUMvo5UkCfv371e2RX/66ScMGzZMeT6YmJgIrZbXJIioe5nNZvzxj3/Eli1b+Mu4YELfwA6cf6Znc7iQEBGM+V3wBnabzYaioiJlRdjQ0KCsBjMzMxEZGdmpfx8R0b9zOByYMGECZs2ahRtuuEF0HJ/XqcUHnN/23GSuwJofyn9+a4MGBr2uxT0/SZJhd7nhcsvQ6TSYfl0cJpliO+VwzNVUVlaioKAABQUFMJvNiIiIUFaDqampPFpMRJ0uJycHJSUlWLx4segohC4ovgsa7E5s21+NrSWncLTGBhkytBoNJFmGBhokRBpxe0oMbkmKFvbePUmScODAAeWQzMGDB5GcnKysCAcNGsRtUSLqkFOnTmHy5MlYv349+vTpIzoOoQuL72Ju6fzAaadbgp9Oi+hQgyrfutDY2KhsixYUFODs2bMwmUzKirA1A7aJiC72wgsvYMiQIZg5c6boKPSzbik+T1VVVaVsixYUFCA8PFwpwbS0NG6LEtEVff/991i4cCE2b94Mf39/0XHoZyy+VpIkCQcPHlRWgwcOHMDQoUOVQzJDhgzhtigRKRwOB8aNG4eXXnoJo0aNEh2HLsLia6fGxkYUFxcrp0Xr6uqUbdHMzEz07t1bdEQiEmjlypUoLS3F/PnzRUehf8Pi6yQWi0UpQbPZjLCwMGVbND09HUFBQaIjElE3OXHiBKZNm4YNGzbwl2AVYvF1AUmScOjQIWVbdN++fRgyZIhShNwWJfJesixj1qxZSE1NxfTp00XHoUtg8XWDpqamFtuitbW1MJlMyrWJmJgY0RGJqJP885//xLJly7Bp0yb4+Ym5qkVXxuITwGKxwGw2KyvCkJCQFtuiRqNRdEQiage73Y5x48bhd7/7HUwmk+g4dBksPsEkSUJpaamyGty7dy8GDRqkFGFSUhK3RYk8xLvvvosTJ07gjTfeEB2FroDFpzJ2u13ZFi0oKIDFYkFGRoZShJz8QKROFRUVmDFjBjZt2sRhFyrH4lO5mpoamM1mZaya0WhUrkxkZGQgODhYdEQinyfLMp5++mmMGjUKDz74oOg4dBUsPg9yYVv0QgmWlJRg4MCByiGZ5ORk6HQ60TGJfM4333yDFStWYOPGjdDr9aLj0FWw+DxYc3Mzdu3apTwfrKqqarEt2rdvX9ERibxeY2Mjxo4diz/+8Y9IS0sTHYdagcXnRU6fPq2sBgsKCmAwGJQSzMjIQEhIiOiIRF5n6dKlqK2txf/8z/+IjkKtxOLzUrIs48iRI0oR7t69G4mJicrzwWHDhnFblKiDjh49ikcffRRbtmxBeHi46DjUSiw+H+FwOFpsi1ZWViIjI0N5PtivXz9oNOp7VRSRWsmyjMcffxw33XQTxo8fLzoOtQGLz0fV1dUpl+jz8/Ph7+/fYls0NDRUdEQiVfvqq6+wbt06rF+/nrsnHobFR5BlGWVlZUoJ7tq1CwkJCcq2aEpKCk+qEV3EZrNh7NixmDdvHoYPHy46DrURi49+weFwYPfu3crzwePHjyM9PV1ZEfbv35/bouTTFi9ejHPnzuG///u/RUehdmDx0VWdOXMGO3bsUFaEOp1OWQ2aTCZui5JPOXz4MJ588kls2bIFPXv2FB2H2oHFR21yYVv0wmqwuLgY8fHxyiGZlJQUTqQnryVJEh555BHcfvvtuO+++0THoXZi8VGHOBwOlJSUKHcHjx07hrS0NGVFeM0113BblLzG559/ji1btmDNmjUcHu/BWHzUqc6ePdtitigA5dmgyWRCWFiY4IRE7WO1WjFu3Di89dZbSEpKEh2HOoDFR11GlmUcO3ZMKcGdO3fimmuuUVaDw4cP57YoeYx58+ZBkiS8/PLLoqNQB7H4qNs4nU5lWzQ/Px/l5eVITU1VVoRxcXHcFiVVOnDgAJ577jnk5eXxMJcXYPGRMPX19S1Oi0qSpJTgyJEjeWKOVEGSJDz00EO47777cNddd4mOQ52AxUeqIMsyjh8/rpRgUVER+vfvr2yLXnvttfD39xcdk3zQX/7yF3z22WdYvXo1D7R4CRYfqZLL5UJJSYnyfPDo0aMYMWKEsiKMj4/ntih1ufr6eowbNw7vvPMOBg0aJDoOdRIWH3kEq9WKwsJCZUXodDqVu4Mmk4mT8alL/OlPf4K/vz/mzJkjOgp1IhYfeRxZlnHixIkW26J9+vRRVoMjRozgtih12N69e/HCCy8gNzeX77L0Miw+8ngulwv79u1TLtGXlpbi2muvVZ4PDhgwgNui1CaSJGHq1KmYNGkSxowZIzoOdTIWH3mdhoYG7NixQ3k+2NzcjMzMTGVrlNuidDW5ubnYtm0bVqxYwV+avBCLj7zeiRMnlBIsLCxETEyMUoIjRoxAQECA6IikInV1dRg/fjzef/99DBgwQHQc6gIsPvIpbrdb2RbNz89HaWkpUlJSlOeDiYmJ/A3fx/3hD39AWFgYZs2aJToKdREWH/m0c+fOKadFCwoKYLPZWswWjYiIEB2RutHu3bvx8ssvIy8vD0FBQaLjUBdh8RFd5OTJky22RaOiopRt0dTUVBgMBtERqYu43W5MnjwZM2bMwOjRo0XHoS7E4iO6DLfbjQMHDijboocOHcKwYcNabItykof32LhxI7777jssX76c291ejsVH1Eo2mw1FRUVKEZ47dw4mk0m5NhEZGSk6IrVTbW0txo8fjw8++ABxcXGi41AXY/ERtVNlZaWyLbpjxw5EREQoq8G0tDRui3qQV199Fb1798bTTz8tOgp1AxYfUSeQJAkHDhxQivDgwYNITk5Wng8OGjSI26IqVVhYiNdeew25ubkIDAwUHYe6AYuPqAs0Nja22Ba1Wq0wmUxKEUZFRYmOSDg/9WfixIl44okncNNNN4mOQ92ExUfUDaqqqpTVoNlsRq9evZQSTEtL40pDkHXr1qGwsBBvv/02D7T4EBYfUTeTJAkHDx5U7g4eOHAAQ4cOVQ7JDBkyhNui3cBisWDixIlYs2YN+vfvLzoOdSMWH5FgjY2N2Llzp7IirKurazFbNDo6WnRErzR37lzExcXh8ccfFx2FuhmLj0hlLBaL8mzQbDajR48eSgmmp6dzokgnyM/PxxtvvIHc3FzOavVBLD4iFZMkCYcOHVKKcP/+/RgyZIhybYLbom3ncDgwYcIEPP/88/j1r38tOg4JwOIj8iBNTU0oLi5WirC2trbFJfqYmBjREVXvww8/xN69e7F48WLRUUgQFh+RB7NYLCgoKFD+FxISoqwG09PTYTQaRUdUlcrKSkyZMgXr169Hnz59RMchQVh8RF5CkiQcPnxYOSSzd+9eDB48WHk+mJSU5PPborNnz0ZSUhIefvhh0VFIIBYfkZey2+0ttkVramqQkZGhrAh9bcXz3XffYdGiRdi8eTP8/f1FxyGBWHxEPqKmpgZms1lZERqNRqUEMzIyvHpbtLm5GQ888ADmzp2L6667TnQcEozFR+SDJElCaWmpUoIlJSUYOHCgckgmOTkZOp1OdMxOs3LlSpSWlmL+/Pmio5AKsPiICM3NzSguLlaKsKqqCiNHjlSeD/bt21d0xHY7ceIEpk2bho0bN3IYAAFg8RHRJZw+fVopwYKCAgQGBiolmJGRgZCQENERFW5JRrXVDqdbgp9Oi+hQA3Ta83M3ZVnGrFmzkJaWhmnTpglOSmrB4iOiK5JlGUeOHFFKcPfu3UhMTFSeD4rYFrXandi2vxpbS06hrMYGaGRooIEMGbIMJEQG4/aUGATVHcaHK97Fxo0b4efn160ZSb1YfETUJg6HA7t27VJOi1ZWViqnRTMzM9GvX78ue9OB0y1ho7kCa38oh1uSoddpYNDroNX+6++TJBl2lxsOl4TqU5WYktkfc8f9Gn46377KQf/C4iOiDqmrq1Mu0Ofn58Pf37/FadHQ0NBO+XssDXbMyduDstpzMPrrob9KkVksFjQ7XQgNj0B8RDAWjB2OqBBDp2Qhz8biI6JOI8syysrKlNXgrl27kJCQoBThsGHDoNfr2/x1LQ12PLa+CHU2B0IM+quuKB2OZpSXlyMhIQE6nR4NdhfCjf5YMSWd5UcsPiLqOg6HA7t371aeDx4/fhzp6elKEfbv3/+qJeZ0S5i5rhDltTaEBrbmOZ2MY8cqEBISjPDwXsqfWpuciI8wYvXUjKuuFsm7sfiIqNucOXOmxSV6nU6nlODIkSMvuS269sdyfPBdGUJbsdIDAKvVitraGsTHJ7T4eFmWYbU7MfNXCZh6XVxnflvkYVh8RCTEhW3RCyVYXFyM+Ph45dpESkoKmtzAPcu/R4Bee8lVmuRyonTbRzh77ABc9kYEhEXAb0AmEjN+c8n3FrrcEppdEv721PUIMfCUp69i8RGRKjgcDpSUlCjPBysqKhCRMQblwUmICAv6eb5myxWf29GMEzv+f0QPG4WAkHAcKfonTny7EVmPvQlDWMQl/x5rkxPP/XYg7kvr1w3fFakRi4+IVOns2bOYuvoHlJ9uhKOxAQAQHGyE0RgMo9H4i7uDzc3NOHbsGKzffYS46+9GxOD0S37dRocLcb2M+GDayC7/Hkid2n68ioioG4SEhqEBQegfEwKtNhrNzQ7YbDbU19fj1KlK+PsHwGg0IjjYiMDAQFRVVSEsKAA1Z2sQFHH5N08Y9DocrbHBLcnKhBfyLSw+IlKlaqsd0MjK5fSAgAAEBAQgPDwcsiyjqakJ586dQ3W1BXa7HYGGAFTn/wVRydchqNfl30Sv1Z6f8FJttaNPj8Du+nZIRVh8RKRKTrcEDS69ItNoNAgKClIOsLicTpR++SG0Oh0G3Dzpql9bq9HA6ZY6NS95Dl5mISJV8tNpIePqRxBkWcbRbR/B0WjF0LufhFZ39d/nJVnmCDMfxv/yRKRK0aEGQNZAkq5cfqV/X4/G05VIvu8Z6Pyu/mZ1STo/0Do6lBNcfBW3OolIlXRaDeIjjTh22oYg/0v/qLLX16Jq93Zo9X4oWD5b+fPE0VMQlZx16c9xuZEQaeTBFh/G6wxEpFp/3nkCS7853MpRZa3De3zErU4iUq3RSdHQaTVwddJBFJdbgk6rwS1JfBO7L2PxEZFqhRj8MG1UHGwOFzq6OSXLMmwOF6aPiuO4Mh/H4iMiVZtkikV8RDAa7K4OfZ0GuwsJEcGYZIrtpGTkqVh8RKRqfjotFowdjnCjP6xNzjav/GRZhrXJiXCjP+aPHc5XEhEPtxCRZ2jrG9iB88/0bI7zK735fAM7/YzFR0Qew+mWsMlcgTU/lMMtydDrNDDodcpYM+D8PT27yw2XW4ZOp8H06+IwyRTLlR4pWHxE5HEa7E5s21+NrSWncLTGBhkytBoNJPn85fSESCNuT4nBLUnRPMhCv8DiIyKP5pbOD5x2uiX46bSIDjXwcjpdEYuPiIh8Cje9iYjIp7D4iIjIp7D4iIjIp7D4iIjIp7D4iIjIp7D4iIjIp/w/jjaBefnS6z0AAAAASUVORK5CYII=\n",
"text/plain": [
"