tutorials/combinatorial_optimization/MAXCUT_CN.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# QAOA 求解最大割问题\n",
    "\n",
    "<em> Copyright (c) 2021 Institute for Quantum Computing, Baidu Inc. All Rights Reserved. </em>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 概览\n",
    "\n",
    "在[量子近似优化算法教程](./QAOA_CN.ipynb)中，我们介绍了如何将经典的组合优化问题编码为量子优化问题，并用量子近似优化算法 [1]（quantum approximate optimization algorithm, QAOA）求解。在本教程中，我们将以最大割问题为例来进一步阐述 QAOA。\n",
    "\n",
    "### 最大割问题\n",
    "\n",
    "最大割问题（Max-Cut Problem）是图论中常见的一个组合优化问题，在统计物理学和电路设计中都有重要应用。最大割问题是一个 NP 困难问题，因此目前并不存在一个高效的算法能完美地解决该问题。\n",
    "\n",
    "在图论中，一个图是由一对集合 $G=(V, E)$ 表示，其中集合 $V$ 中的元素为该图的顶点，集合 $E$ 中的每个元素是一对顶点，表示连接这两个顶点的一条边。例如下方图片中的图可以由 $V=\\{0,1,2,3\\}$ 和 $E=\\{(0,1),(1,2),(2,3),(3,0)\\}$ 表示。\n",
    "\n",
    "![G](figures/maxcut-fig-maxcut_g.png \"图 1：一个有四个顶点和四条边的图\")\n",
    "<div style=\"text-align:center\">图 1：一个有四个顶点和四条边的图 </div>\n",
    "\n",
    "\n",
    "一个图上的割（cut）是指将该图的顶点集 $V$ 分割成两个互不相交的集合的一种划分，每个割都对应一个边的集合，这些边的两个顶点被划分在不同的集合中。于是我们可以将这个割的大小定义为这个边的集合的大小，即被割开的边的条数。最大割问题就是要找到一个割使得被割开的边的条数最多。图 2 展示了图 1 中图的一个最大割，该最大割的大小为 $4$，即割开了图中所有的边。\n",
    "\n",
    "![Max cut on G](figures/maxcut-fig-maxcut_cut.png \"图 2：图 1 中图的一个最大割\")\n",
    "<div style=\"text-align:center\">图 2：图 1 中图的一个最大割 </div>\n",
    "\n",
    "\n",
    "假设输入的图 $G=(V, E)$ 有 $n=|V|$ 个顶点和 $m=|E|$ 条边，那么我们可以将最大割问题描述为 $n$ 个比特和 $m$ 个子句的组合优化问题。每个比特对应图 $G$ 中的一个顶点 $v$，其取值 $z_v$ 为 $0$ 或 $1$，分别对应该顶点属于集合 $S_{0}$ 或 $S_{1}$，因此这 $n$ 个比特的每种取值 $z$ 都对应一个割。每个子句则对应图 $G$ 中的一条边 $(u,v)$，一个子句要求其对应的边连接的两个顶点的取值不同，即 $z_u\\neq z_v$，表示该条边被割开。也就是说，当该条边连接这的两个顶点被割划分到不同的集合上时，我们说该子句被满足。因此，对于图 $G$ 中的每条边 $(u,v)$，我们有\n",
    "\n",
    "$$\n",
    "C_{(u,v)}(z) = z_u+z_v-2z_uz_v,\n",
    "\\tag{1}\n",
    "$$\n",
    "\n",
    "其中 $C_{(u,v)}(z) = 1$ 当且仅当该条边被割开。否则，该函数等于 $0$。整个组合优化问题的目标函数是\n",
    "\n",
    "$$\n",
    "C(z) = \\sum_{(u,v)\\in E}C_{(u,v)}(z) = \\sum_{(u,v)\\in E}z_u+z_v-2z_uz_v.\n",
    "\\tag{2}\n",
    "$$\n",
    "\n",
    "因此，解决最大割问题就是要找到一个取值 $z$ 使得公式（2）中的目标函数最大。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 编码最大割问题\n",
    "\n",
    "为了将最大割问题转化为一个量子问题，我们要用到 $n$ 个量子比特，每个量子比特对应图 $G$ 中的一个顶点。一个量子比特处于量子态 $|0\\rangle$ 或 $|1\\rangle$，表示其对应的顶点属于集合 $S_{0}$ 或 $S_{1}$。值得注意的是，$|0\\rangle$ 和 $|1\\rangle$ 是 Pauli $Z$ 门的两个本征态，并且它们的本征值分别为 $1$ 和 $-1$，即\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "Z|0\\rangle&=|0\\rangle,\\tag{3}\\\\\n",
    "Z|1\\rangle&=-|1\\rangle.\\tag{4}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "因此我们可以使用 Pauli $Z$ 门来构建该最大割问题的哈密顿量 $H_C$。因为通过映射 $f(x):x\\to(x+1)/2$ 可以将 $-1$ 映射到 $0$ 上 并且仍将 $1$ 映射到 $1$ 上，所以我们可以将式（2）中的 $z$ 替换为 $(Z+I)/2$（$I$ 是单位矩阵），得到原问题目标函数对应的哈密顿量\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "H_C &= \\sum_{(u,v)\\in E} \\frac{Z_u+I}{2} + \\frac{Z_v+I}{2} - 2\\cdot\\frac{Z_u+I}{2}\\frac{Z_v+I}{2}\\tag{5}\\\\\n",
    "&= \\sum_{(u,v)\\in E} \\frac{Z_u+Z_v+2I - (Z_uZ_v+Z_u+Z_v+I)}{2}\\tag{6}\\\\\n",
    "&= \\sum_{(u,v)\\in E} \\frac{I - Z_uZ_v}{2}.\\tag{7}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "该哈密顿量关于一个量子态 $|\\psi\\rangle$ 的期望值为\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\langle\\psi|H_C|\\psi\\rangle &= \\langle\\psi|\\sum_{(u,v)\\in E} \\frac{I - Z_uZ_v}{2}|\\psi\\rangle\\tag{8}\\\\\n",
    "&= \\langle\\psi|\\sum_{(u,v)\\in E} \\frac{I}{2}|\\psi\\rangle - \\langle\\psi|\\sum_{(u,v)\\in E} \\frac{Z_uZ_v}{2}|\\psi\\rangle\\tag{9}\\\\\n",
    "&= \\frac{|E|}{2} - \\frac{1}{2}\\langle\\psi|\\sum_{(u,v)\\in E} Z_uZ_v|\\psi\\rangle.\\tag{10}\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "如果我们记\n",
    "\n",
    "$$\n",
    "H_D = -\\sum_{(u,v)\\in E} Z_uZ_v,\n",
    "\\tag{11}\n",
    "$$\n",
    "\n",
    "那么找到量子态 $|\\psi\\rangle$ 使得 $\\langle\\psi|H_C|\\psi\\rangle$ 最大等价于找到量子态 $|\\psi\\rangle$ 使得 $\\langle\\psi|H_D|\\psi\\rangle$ 最大。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Paddle Quantum 实现\n",
    "\n",
    "接下来，我们用量桨实现 QAOA 来求解最大割问题。有许多方法可以找到参数 $\\vec{\\gamma},\\vec{\\beta}$，我们这里使用经典机器学习中的梯度下降方法。\n",
    "\n",
    "要在量桨上实现 QAOA，首先要做的便是加载需要用到的包。其中 `networkx` 包可以帮助我们方便地处理图。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-04-30T09:11:04.006915Z",
     "start_time": "2021-04-30T09:11:03.951557Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<style>pre { white-space: pre !important; }</style>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from IPython.core.display import HTML\n",
    "display(HTML(\"<style>pre { white-space: pre !important; }</style>\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-04-30T09:11:08.078417Z",
     "start_time": "2021-04-30T09:11:04.569758Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/Caskroom/miniconda/base/envs/pq_new/lib/python3.8/site-packages/paddle/tensor/creation.py:125: DeprecationWarning: `np.object` is a deprecated alias for the builtin `object`. To silence this warning, use `object` by itself. Doing this will not modify any behavior and is safe. \n",
      "Deprecated in NumPy 1.20; for more details and guidance: https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations\n",
      "  if data.dtype == np.object:\n",
      "/usr/local/Caskroom/miniconda/base/envs/pq_new/lib/python3.8/site-packages/paddle/tensor/creation.py:125: DeprecationWarning: `np.object` is a deprecated alias for the builtin `object`. To silence this warning, use `object` by itself. Doing this will not modify any behavior and is safe. \n",
      "Deprecated in NumPy 1.20; for more details and guidance: https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations\n",
      "  if data.dtype == np.object:\n"
     ]
    }
   ],
   "source": [
    "# 加载量桨、飞桨的相关模块\n",
    "import paddle\n",
    "from paddle_quantum.ansatz import Circuit\n",
    "from paddle_quantum.qinfo import pauli_str_to_matrix\n",
    "from paddle_quantum.loss import ExpecVal\n",
    "from paddle_quantum import Hamiltonian\n",
    "\n",
    "# 加载额外需要用到的包\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import networkx as nx\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "接下来，我们生成该最大割问题中的图 $G$。为了运算方便，这里的顶点从 $0$ 开始计数。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-04-30T09:11:08.411878Z",
     "start_time": "2021-04-30T09:11:08.093215Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAADnCAYAAAC9roUQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAjJUlEQVR4nO3deXSV5bX48e8OgYQwxoIMSk1FQHEoqFWrgAOt/Gpaa53Qqz8rTtU6UHpRm94Otl6N+oNqh1utUKWr1mqv97Y/Nd6iFSesI1JHhAoeFVEEGQLkBBKy7x/7TYgMQnLe8w7n7M9arCUvJ+/Zy5xnn32e8zz7EVXFOedcNEriDsA554qJJ13nnIuQJ13nnIuQJ13nnIuQJ13nnIuQJ13nnIuQJ13nnIuQJ13nnIuQJ13nnIuQJ13nnIuQJ13nnIuQJ13nnIuQJ13nnIuQJ13nnIuQJ13nnIuQJ13nnIuQJ13nnIuQJ13nnIuQJ13nnIuQJ13nnItQadwBOLerqmrqSoChwHCgO9AN2ARkgUXA4kxtdUt8ETq3c+KnAbukCpLseKAaGAvsB7QAzYAEfzT4U4p9clsAPAXUAY96EnZJ40nXJU5VTV0lMAmYCvQCemAJdlcpsAGoB6YDd2Zqq1eHHadzneFJ1yVGVU1dBXAjcAFW0VaEcNsGrAKeCVydqa1uCOGeznWaJ12XCFU1dWOBe4BKbL42bFlgNTAxU1s9Nw/3d26XeNJ1saqqqSsDbgbOJT/JdmtZYBYwJVNbvTGC53PuEzzputhU1dT1BB4GRhFNwm2VBeYDEzK11esjfF7nPOm6eAQJdy4wAiiPIYRGYCEwxhOvi5JvjnCRC6YUHia+hEvwvCOA2UE8zkXCk66Lw83YlEJcCbdVOTAa+FnMcbgi4tMLLlLBKoXZRDuHuzNZ4Hhf1eCi4EnXRSZYh/tPYHDcsWzHMmCYr+N1+ebTCy5KN2HrcJOoErgh7iBc4fNK10Ui2Nq7jPjncT9NIzDYtwy7fPJK10VlEra1N8lasE0azuWNV7ou74JuYUuBQXHHsguWAUO8O5nLF++n66IwHusW1mlj9unHsSP6c8DgPowc3Jte5V3b/u2M25/h2bdX5Rpjq97AccDfwrqhc+150nVRqMbaM3baOUfsxfH7DwwpnE9VAZyAJ12XJ550XRTG0rF+uNtQ4MO1jby2bC3rNzZz0qg9wolsWyXAuHzd3Dmf03V5FcznbiDHVQvlXUtobLJp1iM+txv3XPTFtn8LeXoBbLNEj0xttQ8OFzpfveDybSiwOdebtCbciLRgcTsXOk+6Lt+GY2eapUkzFrdzofOk6/KtOznO58ZASFZvCFdAPOm6fOtGOpOut3t0eeFJ1+XbJmzxQZoo4Ef5uLzwpOvyLUs6k2427iBcYfKk6/JtEelbD16Kxe1c6NI2GFz6LCaEN/evHjiIg/bsC8Dgvp9c8nv2EXtx3L4DAHhl6RoefPWDXJ+uBIvbudD55giXd1U1dfOAg3O5x7RTD+LUQ4bs9HH3zXuPqfe9kstTAczL1FYfmutNnNsen15wUXiK9MzrtgBPxh2EK1w+veCiUAecD/Ts7A2m3vdKGBXsrmgAHoriiVxx8krXReFRYF3cQeyiemBO3EG4wuVJ1+Vd0BB8GlZFJlkDMM0bmLt88qTronInyX+9lQCz4g7CFbakDwJXIDK11as3LV/yaEvTxkR+oaaqWWCGH0rp8s2Trss7ERkgIvd++Psrq1sa1yeyD0NLw9ouH9w5eVbccbjC50nX5Y2Ys4E3gNO1eeOGtc/86ZagqkyMlqaNuuK/r+u2afniZ0TkhyLSLe6YXOHypOvyQkQ+iy0V+z2wG/AIcMC6eQ9OEZFZJKe3QRb0txvfX3A71hHtp8ALIuKbI1xeeNJ1oRKREhG5BHgd+AqwBpgETFDVTPCwKcB8oDGOGNtpBOaXdC2/TFW/hZ0CvBg4CHhORG4SEe+r60Ll24BdaERkODATO4gS4L+BS1X1w60fW1VT1xOYC4wgx/PTOqkRWAiMydRWr2+9KCIVWLU7BStK3gIuUNUnYojRFSCvdF3ORKRURK4CXsYS7nLgVFU9ZXsJFyBIdGOwijfqqYYs8BJbJVwAVW1Q1anAF4HXgH2Ax0XkVhHpHXGcrgB50nU5EZHPA88BN2IV6++Akar6Xzv72SDhHYut4Y0q8WaD5ztu64Tbnqo+DxwCXAM0ARcDr4tIdRRBusLl0wuuU0SkHPgBcDXWw+Nd4CJVnd2Z+1XV1I0B7gUqyc/5ZFlgNTAxU1s9tyM/KCIHAL8FDgsu/QH4jqquDDdEVwy80nUdJiJHYtMC/wZ0AX4JHNDZhAsQJMJh2JxwI+FtGW4I7jcTGNbRhAugqq8BRwLfxZL3WcACETlDRBK57tgll1e6bpeJSE/gOuBy7PDGhcD5qvp0mM9TVVNXCZwLTAV6AxV0rEBowZJtPdbzYVZYO81EZCgwA5sWAbgf+Laqvh/G/V3h86TrdomIfBm4HagCNmNzuNeqat6WfVXV1JUA47GlZ+OAkVhCbW5p3NAHEUrKKuqxXr2lWGJ+A+uH+xAwJx/Na4Lq9nxgOvamUI+9QcxUH1BuJzzpuk8lIpVYcpkUXJoPnKeq/4g6liAJ7w0MX/Hn6+vo0pX+J155KvaRfxGwOFNbHdkLWkT2AG4FvhZcegy4UFX9qB+3Q5503Q6JyMnAfwADsSPJrwGmq2pTnHEBiIgCqGqsc6pB1Xs6Nq/dH3sD+AHwc1XdHGdsLpk86bptiMhA4FfAKcGludgGgYXxRfVJSUm6rUSkH3AL9iUbwPPYfPdrsQXlEslXL7g2QYOab2LzoqcA64HLgKOTlHCTSFVXqurZwFeBpdjyspdE5MfeQMe155WuA0BE9gJ+A0wILv0VuFhV34kvqh1LWqXbXrBz7QbgkuDSa1jV+3x8Ubmk8Eq3yAUNai7DGtRMAFYB5wAnJDXhJp2q1qvqt4FjsN4NBwDPiMj0oLeDK2Je6RYxEdkX2zRwVHDpP4HLVXV5fFHtmiRXuu0FXcquwZaUlQBLsPnxx+KMy8XHK90iJCJdReT7WIOao4APgZNV9fQ0JNw0UdWsql4NHA68gi15myMit4tIn3ijc3HwSrfIiMho4A5gVHDpDmCqqqbqbLC0VLrtiUhXrFfFD7GG6cuwefMHYg3MRcqTbpEIGtT8GLgS65eQwRby/y3OuDorjUm3lYiMxBroHBFcuge4QlVXxBeVi4pPLxQBERmDTSV8D/ud3wIcmNaEm3aq+gbWS/g7WI+IM7AGOmd5A53C55VuARORXkAtcGlwaQG2dOmZ+KIKR5or3fZE5HNYT4svBZfqgEtU9b34onL55JVugRKRCdj60EuBZuBaYHQhJNxCoqpvA8djDXTWAtVYs/SLRcTHZwHySrfAiMhngJ9ha20B5mHV7cvxRRW+Qql02xORwVivi5OCS09g8+7/jC0oFzp/Jy0QwRbeU7EtvOdgjbuvAo4otIRbqFR1GXAycBrwEXA08IqIXCkipbEG50LjlW4BEJFBWIX0jeDSk1iFtCi+qPKrECvd9orlE0sx8ko3xYLqdhJW3X4DWIft9z+2kBNuMVDVj1X1m1gD93exQzJfFJFrRaQs3uhcLrzSTantfOv9ELbQvii+9S70Sre9Ql6FUoy80k0ZEekiIldgKxO+BHwMnA18tVgSbrFR1XWqehl2ZNEiYD/gaRG5RUR6xBud6yivdFMk2Mk0E/hicOkeYLKqfhRfVPEopkq3vWBn4Y+wL0lTv7OwGHnSTYEd7Nm/RFXvjzWwGBVr0m0lIgdjW4lHBZfuAP5VVdfEFZPbNZ50E05EDsEG1EHBpRnAVcU+uIo96ULbm/FUrKdGGfABdhz8X+KMy306T7oJtYM+rBeq6pw440oKT7pbBH2RfwscGVxKTV/kYuRfpCWQiIzDGtRcFVz6GdagxhOu24aqvgmMBS4HNmCbKxaIyDneQCd5vNJNkO2crfU6tjToufiiSiavdLdvB2fdfUtV340vKtde0Sbdqpq6EmAoMBzojn1BtQnIYstyFmdqq1uiikdETsAGy55AE3A9cL2qbooqhjTxpLtjQXV7DnAzUImd6vw94FZVjew1nbQxlhRFk3SDF8B4rIvTWGytYwvWgUuCPxr8KcWmXhYAT2Ht9h7NxwtERPph/W3PCi69gFW3r4b9XIXEk+7OichA4JfAqcGludj5bAvz8XxJHWNJU/BJt6qmrhKYhH0h1Qvogf3yd5Vi82T1wHTgzkxtdc5H2wTVyOnYoOiPvfv/ELhFVTfnev9C50l314nIycCvgQHARuwL2umq2hTG/ZM6xpKqYJNuVU1dBXAjcAH2bhvG0dcN2LvzTODqTG11Q2duErTwuxU4Mbj0OLYy4a0QYiwKnnQ7RkQqsYQ2Kbg0H/tENb+z90zyGEuygky6VTV1Y7HdWpXYXFLYssBqYGKmtnrurv5QUN2eD0wD+mDv7FOBmVqIv4g88qTbOSJyPNazYy9gM5Y0r1XVxo7cJ6ljLA0KKulW1dSVYV8enEt+XghbywKzgCmZ2uqNn/ZAEdkb29hwXHDpAWxX2ft5jbBAedLtPBHpCVyHLTETYCFW9T69s59N8hhLi4JJulU1dT2Bh7FtkVG8GFplsY9qEzK11eu3/kcR6QJcgb3IuwMrsRf7vV7ddp4n3dyJyJHYpop9sXnVXwHfV9VtXseQ3DGWNgWRdIMXw1xgBFAeQwiNWLUwpv2LQkQOwF7UhwWX7sYa1KyMPsTC4kk3HEEDnR9gvT1KgXewdb2z2z8uqWMsjVKfdIOPO48Bo4nnxdCqEXgJOO6dG76qQA3wb0BX4H2s1+2DMcZXUDzphktERmEFwsHBpd8B31XVVUkcY2meaiiEbcA3Yx934nwxEDz/6KY1H/4eO1rlGizh3gbs7wnXJZmq/gM4HNtEsRH4JvCGiJxCwsYYti0+tVKddINvUM8l2vmlT9O9S4/K08r2HHkA8BZwjKpeoqpr4w7MuZ1R1WZVvRHraPcUMKBsz/3v0+ami0jQGAMmVdXUjYk7kM5K7fRCsEbwn8DguGPZWkvjhvX18+6vWvPkXR/HHUuh8umF/BKRkpKKvpcNOu8Xt5T23C2J/4+XAcPSuI43zZXuTdgawcQpKe/Rpe9RZ/447jic6yxVbRlyxV3Du/SoTOrcaSXWHCp1UlnpBtsOlxH/HNOnaQQGF/J2xjh5pZtfPsbyJ62V7iRs22GStWDzzc6lkY+xPEldpRt0MloKDOrsPXqXl/LlkQM4rOoz7D+4N/17lVFZ0Y2NzZvJfNzAnDc/4o6n32ZtNud+IMuAIcXQOSlqXunmTxhjDKBvRVcuHLM34/fbnSGV1pbhvdUNPLrgI25/akkY4wtSOMZK4w6gE8ZjnYw67ah9+jH9tFHbXO9WWsKBe/ThwD36cMYXhvAvM59l8YoNuTxVb2zbr5/U6tIk5zE2fEBPfn/e4Qzo/cnZiX0H9mbfgb059ZA9+b93PMei5Tnvc0jdGEvj9EI11jouZ/XZJh54eRnTH1nI7U8uYXn9lp4fA3qXc/1JB+b6FBXACbnexLmI5TTGykpLuO3sQ9oS7tpsE7c9sZjbnljcVt0O6F3OrWcdQllpzikodWMsjZXuWDrWq3Mbaxqa+MkDr/PHF96lsWnLp5LbnlzM7Mnj6N+rDIAvVO1Gj25d2LCp0+1tS4BxucTqXAxyGmMnjdqDvfv1bPv75Hvm8/iiFQA8+/bHzDrXdsUP7d+Tr4/agz+9+F4usaZujKWq0g3mmkbmep9nlnzMnX/PfCLhAqzasIkXMqva/l5SInTN/Z14ZFVNnc87ulQIY4xN2H9g23/XNza1JVyAJxatYF3jlrnc/9PusTlI1RhLVdLFzlvK66kKQ/tveYfOfLyBNQ05T/a3YHE7lwY5j7GRg3q3/ffSVZ/cu6AKS1dn2/6+36Ccpo5bpWqMpS3pDsfOW8qLyeOHMWLglhfBzx5ZFMZtm7G4nUuDnMdYZUXXtv9et3HbW61r3HJtt4puuTxVq1SNsbTN6XYnx/nc7RGBfzthPy4Ys3fbtVv+toj7X14Wyu1Jzr5153Ym1DEm27mVhD8RkKoxlrak242Qk26Pbl34xZmjGb/vAABaWpTav77JjKeWhPUUApSFdTPn8iznMba6oYmBfboA0Kt82xTTs2zLtVUNm3J5qlapGmNpS7qbsA73odijb3dmnnMo+wVzUA2bmpnyp5eZ/fqHYT0FWLxJ3b/u3NZyHmNvfFDPwD62XGzPyu6I2FwuWJU7ZLct51cu+GBdLk/VKlVjLG1zullCSrqjh/TlL98+si3hLluT5bTbngk74YLFm93po5xLhpzH2MNvbBlDvcq7cszw3dv+fszw3T9R6YY03lI1xtJW6S4ihJgP/mwld19wOOVd7SNQ8+YWHnzlA47cpx9H7tPvE4998JVlfLC2Qwelbq0Ui9u5NMh5jP15/vtcNG7vtrW6P584iruffxeAMw/7bNvjlqxcz1/+Ecq5rKkaY2lLuosJoTrfu1+PtoQLUNqlhIvG7b3dx766dE2uSbcEi9u5NMh5jG1sbuHiu+Zx13mHs3vvcnp378rFR39yRddH9Y1cfNc8NjaH0jIhVWMsVdMLQVOLBXHH0UFvZGqr09VVyBWtsMbYouXrOf7nT/Lrx99i0fJ1NGxqpmFTM4uWr+PXj7/F8T9/Moy+C61SNcbSVumCHSMymhy+Yb3vpaXc99LS8CLasRbgySieyLkQ5TzGwLbb3zR7ITfNXhhOVNuXujGWqko3UAfk1PorQg3AQ3EH4VwH+RjLozQm3UeBUNaZRKAemBN3EM51kI+xPEpd0g3mnKZh73BJ1gBMS1NzZefAx1i+pS7pBu4k+bGXALPiDsK5TvIxlidJ/5+6XcFBdDNJ6IJo3dy8SVVnpO3APOdaZWqrVzetfO+vLU0bk7oqIAukcoylMukGrgYS+T98c8Oabu/dfPr+IvK5uGNxrqNEpK+IzPhg1ndOamlcn9Q+tauA78UdRGekNulmaqsbgIkkrNrVls2bVt4/ba1uyh4HvCYik0Wky05/0LkEEJGvA28AF2jzxk1rn/rDDFVN1BjDxvzEIAekTmqTLkCmtnouNqeTlBdFVkq6zNz43mvDgHuw85tuAeaKSM4nXjiXLyKyu4jcA/wFOwX4GWDUupdnXyQis0jQGAPuzNRWPx13IJ2V6qQbmALMB3LaqxuCxiCO76rqClU9E/g6dkT0EcB8EfmhiITStdm5MIg5G9uFNhFbETAZGKuqrTvTEjfGYo4jJ6Ka1HnyXVdVU9cTmAuMAMp38vB8aAQWAmMytdWf2NsoIn2A/wdcGFx6BThfVV+MNsTCIiIKoKpJnXNMPBEZAtzGltN0HwG+papvb/3YJI+xtCmESpfglzAGexeM+mNQFniJHbwYVHWtql4EjAeWAAcBz4nITSKSmm73rnCISImIXAK8jiXcNcAkYML2Ei4ke4ylTUEkXWh7URyLrS+M6kWRDZ7vuJ29GFR1DnAgMD24dCXwiogcnd8QndtCRIYDjwG/BnoBfwZGquos3cnH3qSPsbQoiOmFrVXV1I0B7gUqyc/ZSVlsudrE4Mu8DhGRw4DfAgcEl24DrlbV+vBCLGw+vdAxIlKKzYX+BJseWA5cqqr/1Zn7JX2MJVnBVLrtBb+kYdgGikbC287YENxvJjCssy8GVX0eOAS4BmgCLgZeF5HqkOJ0ro2IfB54DrgRS7i/w6rbTiVcSP4YS7KCrHTbq6qpqwTOBaYCvbFlXB15s2nBXgj12H70WWHughGRA7Cq97Dg0h+A76jqyrCeoxB5pbtzIlIG/ADbRFAKvAtcpKqzw3yepI+xpCn4pNuqqqauBPsy6yvAOGAk9stuxvqGCnbWkmIv0BJskfiTWOu4OflqrBFsnpgM/Dv2UW0lcDlw787m2YqVJ91PJyJfxN7M98Ne0/8BfF9V89Y9LMljLEmKJuluLXiB7A0MxxJdGXaiaBY7b2lx1N3oRWQoMAP7sgLgfuDbqhrKQVKFxJPu9olIT+zN+wosyS0ELlDVyD+mJ3GMJUHRJt2kEhEBzsdWOfTGPnJNBWZ61buFJ91ticiXgduBKmAzcBPwU1WNe1ODa8eTbkKJyB7ArcDXgkuPAReqamoO4MsnT7pbiEgl9iY9Kbj0D+A8VZ0fW1Buhwpy9UIhCKYUvg6cAazAphxeFZHvegMd10pEvoHNi07CPrp/HzjME25yeaWbAiLSD2ucc1Zw6XlsK/FrsQUVs2KvdEVkIPBL4NTg0tPY3O2b8UXldoVXuimgqitV9Wzgq8BSbHnZSyLyY2+gU1yCBjXnYNXtqcB64DJgnCfcdPBKN2VEpDe2yP3i4NJrWNX7fHxRRa8YK10R2Qv4DTAhuDQba1DzTnxRuY7ySjdlVLVeVS8BjgHewrYSPyMi00SkItbgXF4EDWouwxrUTMC2x34T+Ion3PTxSjfFgi5l12BLykqwLmYXqOpjccYVhWKpdEVkBLbJ4ajg0n3AZaq6PL6oXC680k0xVc2q6tXA4cCr2EL0OSLym6CPr0spEekqIjXAy1jC/RA4WVVP84Sbbl7pFojgC7WrgB8C3bATKy5W1QdiDSxPCrnSFZHRWHU7Orh0BzBVVQu2H0Ex8aRbYIKz2H6LHREEdlbbFaq6Ir6owleISVdEyoEfYW+eXYAMtiHmb3HG5cLl0wsFRlXfwDr8fwfr3HQGsEBE/iXYYuwSSETGYDvJarBx+XPgQE+4hccr3QImIp/D9uJ/KbhUB1yiqu/FF1U4CqXSFZFeQC1waXBpAbYE8Jn4onL55JVuAQvOuzoea6CzFqjGmqV/S0T8dx8zEZmArbO+FGt/+O/AaE+4hc0r3SIhIoOxnqonBZeewOYL/xlbUDlIc6UrIrsBNwPnBJfmYdXty/FF5aLi1U6RUNVlwMnA6cBHwNHYwZhXBudnuQiIyKnYFMI52LE0VwFHeMItHl7pFiER+QzwM7ZUWi9ildYr8UXVMWmrdEVkEPAr7I0P7LSEC1V1UXxRuTh4pVuEVPVjVf0mdqzKu8ChwDwR+WlwrpYLSdCgZhLWoOZkYB1wCXCsJ9zi5JVukdvOt+dvYFXvs/FFtXNpqHS3s3rkf7AGNalfPeI6zyvdIqeq61T1MuwgwUXYYYJ/F5GbRaRHvNGlk4h0EZErsJUJXwI+Bs4Gqj3hOq90XZtgR9SPgSuxHVFvY0d2J26BflIrXRHZD9sR+MXg0j3AZFX9KL6oXJJ40nXbEJGDscQxKrh0B/Cvqromrpi2lrSkKyJdsZUIP2JL74tLVPX+WANzieNJ121XkESuxCrfbsAH2HHwf4kzrlZJSroicgj2xnRQcGkGcFWS3qRccnjSdZ9KRPbFqt4jg0v/CVwed3vBJCTdHfQzvlBV58QVk0s+T7pup4Itw98GbgB6AKuwhjp3aYQvoKqauhJgKDB8xZ9rH6RLKf1PvPI0IIt9Cbg4U1vdEkUsIjIOmAkMA1qwg0N/pKobonh+l16edN0uE5Eq7Iyu44NL/4P17H03H88XJNnxWM+IscB+WIJrbtm4oQ8IJWUV9YACpVi1uQB4Cmvu82jYSTg4o+4GbK0t2BE656vqc2E+jytcnnRdhwTtIc/BegdUYqfRXg3cpqqhJLiqmrpKYBL2sb0XVl13ZBpBgQ1APTAduDNTW51zA3AROQG4DRgCNAHXA9er6qZc7+2Khydd1ykiMhDb1npKcGkudj7bws7es6qmrgI76fgCrKIN46DNBqwCnglcnamtbujoDUSkH/Ymc3Zw6QWsun01hPhckfGk63IiIqdg3csGABuxL5amqWpzR+5TVVM3FlvTWgl0DzlMsHnf1cDETG313F35gaCqPx34JdA/uMcPgVtUdXMeYnRFwJOuy1nQqnA6cG5w6SWsEvzHzn62qqauDKsizyU/yXZrWWAWMCVTW71xRw8KWmHeCpwYXHocW5nwVr4DdIXNk64LjYgcj/Ua2AvYjE0VXKuqjdt7fFVNXU/gYWwTRhQJt1UWmA9MyNRWr2//D0F1ez4wDeiDzQtfCcwMa87aFTdPui5UItITuA64HPvy602s6v17+8cFCXcuMAIojzpOrJftQmBMa+IVkb2xjQ3HBY95ENtVtjSG+FyB8oY3LlSqul5VJ2NLvN4E9gXmisgvgoTcOqXwMPElXILnHQHM3v20a7qLyBSsQc1xwErgTOBET7gubJ50XV6o6tPAaKzqbcEq39eCKYibsSmFuBJuq3LVloNbGte/jTV17w7cDeynqvdEufHDFQ9Pui5vVLVRVX+ANUmfD+xVtuf+s3Vz00VEO4e7QyIl5RUjjhxQ9tmDVgBfU9WzVHVl3HG5wuVzui4SIlJa0r339wad96trS3vtFnc421Bt+VCkZGhn1vE61xFe6bpIqGrzkMl3D+zSo+92VzLETaSkD7a917m88krXRSLY2ruM+OdxP00jMDiMLcPO7YhXui4qk7Av1JKshS0bPJzLC690Xd4F3cKWAoM6e48pXxrG/oP7sM/uPams6EaPbl3INm3m/TVZXsys5vfPvsPC5evCCHcZMCSqFpGu+JTGHYArCuOxbmGdNnn88G2u9epSwr4Du7LvwN6cfugQLr37JR5ZkHNv9d7YWt3EnQvnCoMnXReFaqw9Y6etWLeRF99ZxburGljb0ERFWSljh/Xj83v2BaBbaQlXTRgRRtKtAE7Ak67LE0+6Lgpj6Vg/3G184fptc+D0RxbytylHM7R/TwCG7BZGJ0hKsOPoncsL/yLN5VUwnzsyzHuKQN+KrnztoMHs0XfLHos3PwxlThdgZFVNXewHXrrC5JWuy7ehWMexnO3Ztztzrz5uu/+2asMmfvLA62E8DdgqhqGAt3F0ofNK1+XbcKBDDc076p/L13HmjGeZ/96asG7ZjMXtXOg86bp8606O87mt1mSbuO6hBdz41ze54+m3yXxsB+8OG9CL/3/pUZz4+cFhPA1YvInoDeEKj08vuHzrRkhJd/3GZmY8taTt79c9tIDfTTqMMfv0o7xrF244+UCeWfwxK9bv8ECIXSVAWa43cW57vNJ1+bYJO503dJtblEfbLRGr6FbKqCF9w7i1Yue9ORc6T7ou37LkmHQPq9qNPt27bnNdBI4Z0f8T1zSc/K5Y3M6FzqcXXL4tIsfX2emH7snXPj+Y55as4vVla6lvbKayohvHjujPsAFbNrrVNzbx3JJVucYLFu+iMG7k3NY86bp8W0wIn6jKSrswbnh/xg3vv91/X9fYxBV/nM+6jaEslCjB4nYudN7wxuVdVU3dPODgzv78F6oqOeHAQYweUsmgPuX0rbCphvpsE0tWbGDu4pX88fl3Wbl+U1ghz8vUVh8a1s2ca88rXReFp7Dz0jq1iuGFzGpeyETW4rYFeDKqJ3PFx79Ic1GoAzbEHcQuagAeijsIV7g86booPAqE1hghz+qBOXEH4QqXJ12Xd0FD8GlYFZlkDcA0b2Du8smTrovKnST/9VYCzIo7CFfYkj4IXIEIDnucSXI3HWSBGX4opcs3T7ouSlcDSU1qq4DvxR2EK3yedF1kMrXVDcBEklftZoGJQXzO5ZUnXRepTG31XGzeNCmJNwvcmamtfjruQFxx8KTr4jAFmA80xhxHYxDHd2OOwxUR3wbsYlFVU9cTmAuMAMpjCKERWAiMydRWr4/h+V2R8krXxSJIdGOwSjPqqYYs8BKecF0MPOm62AQJ71hsDW9UiTcbPN9xnnBdHHx6wSVCVU3dGOBeoJL8nE+WxZarTQy+zHMuFl7pukQIEuEwbANFI+FtGW4I7jcTGOYJ18XNK12XOFU1dZXAucBUoDdQQccKhBYs2dZjPR9m+U4zlxSedF1iVdXUlQDjga8A44CRWEJtxnrzCnaemWK9oUuAN7B+uA8Bc7x5jUsaT7ouNYIkvDcwHJv3LcNO7c1iZ5otztRW+wvaJZonXeeci5B/keaccxHypOuccxHypOuccxHypOuccxHypOuccxHypOuccxHypOuccxHypOuccxHypOuccxHypOuccxHypOuccxHypOuccxHypOuccxHypOuccxHypOuccxHypOuccxHypOuccxHypOuccxHypOuccxHypOuccxH6X4ZvoaaQCRNmAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# n 是图 G 的顶点数，同时也是量子比特的个数\n",
    "n = 4\n",
    "G = nx.Graph()\n",
    "V = range(n)\n",
    "G.add_nodes_from(V)\n",
    "E = [(0, 1), (1, 2), (2, 3), (3, 0), (1, 3)]\n",
    "G.add_edges_from(E)\n",
    "\n",
    "# 将生成的图 G 打印出来\n",
    "pos = nx.circular_layout(G)\n",
    "options = {\n",
    "    \"with_labels\": True,\n",
    "    \"font_size\": 20,\n",
    "    \"font_weight\": \"bold\",\n",
    "    \"font_color\": \"white\",\n",
    "    \"node_size\": 2000,\n",
    "    \"width\": 2,\n",
    "}\n",
    "nx.draw_networkx(G, pos, **options)\n",
    "ax = plt.gca()\n",
    "ax.margins(0.20)\n",
    "plt.axis(\"off\")\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 编码哈密顿量\n",
    "\n",
    "量桨中，哈密顿量可以以 `list` 的形式输入。这里我们构建式（11）中的哈密顿量 $H_D$。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-04-30T09:11:08.426170Z",
     "start_time": "2021-04-30T09:11:08.418352Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-1.0, 'z0,z1'], [-1.0, 'z1,z2'], [-1.0, 'z2,z3'], [-1.0, 'z3,z0'], [-1.0, 'z1,z3']]\n"
     ]
    }
   ],
   "source": [
    "# 以 list 的形式构建哈密顿量 H_D\n",
    "H_D_list = []\n",
    "for (u, v) in E:\n",
    "    H_D_list.append([-1.0, 'z'+str(u) + ',z' + str(v)])\n",
    "print(H_D_list)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "可以看到，在这个例子中，哈密顿量 $H_D$ 为\n",
    "\n",
    "$$\n",
    "H_D = -Z_0Z_1 - Z_1Z_2 - Z_2Z_3 - Z_3Z_0 - Z_1Z_3.\n",
    "\\tag{12}\n",
    "$$\n",
    "\n",
    "我们可以查看哈密顿量 $H_D$ 的矩阵形式，并且获取它的本征值信息："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-04-30T09:11:08.792299Z",
     "start_time": "2021-04-30T09:11:08.777145Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-5.  1. -1.  1.  1.  3.  1. -1. -1.  1.  3.  1.  1. -1.  1. -5.]\n",
      "H_max: 3.0\n"
     ]
    }
   ],
   "source": [
    "# 将哈密顿量 H_D 从 list 形式转为矩阵形式\n",
    "H_D_matrix = pauli_str_to_matrix(H_D_list, n)\n",
    "# 取出 H_D 对角线上的元素\n",
    "H_D_diag = np.diag(H_D_matrix).real\n",
    "# 获取 H_D 的最大本征值\n",
    "H_max = np.max(H_D_diag)\n",
    "\n",
    "print(H_D_diag)\n",
    "print('H_max:', H_max)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 搭建 QAOA 电路\n",
    "\n",
    "前面我们介绍了 QAOA 需要将两个酉变换 $U_C(\\gamma)$ 和 $U_B(\\beta)$ 交替地作用在初始态 $|s\\rangle = |+\\rangle^{\\otimes n}$ 上。在这里，我们使用量桨中提供的量子门和量子电路模板搭建出一个量子电路来实现这一步骤。要注意的是，在最大割问题中，我们将最大化哈密顿量 $H_C$ 的期望值的问题简化为了最大化哈密顿量 $H_D$ 的期望值的问题，因此要用到的酉变换也就变成了 $U_D(\\gamma)$ 和 $U_B(\\beta)$。通过交替地摆放两个参数可调的电路模块，我们得以搭建QAOA电路\n",
    "\n",
    "$$\n",
    "U_B(\\beta_p)U_D(\\gamma_p)\\cdots U_B(\\beta_1)U_D(\\gamma_1),\n",
    "\\tag{13}\n",
    "$$\n",
    "\n",
    "其中，$U_D(\\gamma) = e^{-i\\gamma H_D}$ 可以由下图中的电路搭建实现。另一个酉变换 $U_B(\\beta)$ 则等价于在每个量子比特上作用一个 $R_x$ 门。\n",
    "\n",
    "![U_D circuit](figures/maxcut-fig-cir_ud.png \"图 3：酉变换 $e^{i\\gamma Z\\otimes Z}$ 的量子电路实现\")\n",
    "<div style=\"text-align:center\">图 3：酉变换 $e^{i\\gamma Z\\otimes Z}$ 的量子电路实现 </div>\n",
    "\n",
    "\n",
    "因此，实现一层酉变换 $U_B(\\beta)U_D(\\gamma)$ 的量子电路如图 4 所示。\n",
    "\n",
    "![U_BU_D circuit](figures/maxcut-fig-cir_ubud.png \"图 4：酉变换 $U_B(\\beta)U_D(\\gamma)$ 的量子电路实现\")\n",
    "<div style=\"text-align:center\">图 4：酉变换 $U_B(\\beta)U_D(\\gamma)$ 的量子电路实现 </div>\n",
    "\n",
    "量桨中，电路运行前每个量子比特默认的初始状态为 $|0\\rangle$（可以通过输入参数来自定义初始状态），我们可以通过添加一层 Hadamard 门使每个量子比特的状态由 $|0\\rangle$ 变为 $|+\\rangle$，由此得到 QAOA 要求的初始态 $|s\\rangle = |+\\rangle^{\\otimes n}$。在量桨中，可以通过调用 `superposition_layer()` 在量子电路中添加一层 Hadamard 门。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def circuit_QAOA(num_qubits, depth, edges, vertices):\n",
    "    # 初始化 n 个量子比特的量子电路\n",
    "    cir = Circuit(num_qubits)\n",
    "    # 制备量子态 |s>\n",
    "    cir.superposition_layer()\n",
    "    # 搭建p层U_D电路\n",
    "    cir.qaoa_layer(edges, vertices, depth)\n",
    "    \n",
    "    return cir\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "搭建 QAOA 量子电路的工作完成后，如果运行该量子电路，得到的输出态为\n",
    "\n",
    "$$\n",
    "|\\vec{\\gamma},\\vec{\\beta}\\rangle = U_B(\\beta_p)U_D(\\gamma_p)\\cdots U_B(\\beta_1)U_D(\\gamma_1)|s\\rangle.\n",
    "\\tag{14}\n",
    "$$\n",
    "\n",
    "### 计算损失函数\n",
    "\n",
    "由上一步搭建的电路的输出态，我们可以计算最大割问题的目标函数\n",
    "\n",
    "$$\n",
    "F_p(\\vec{\\gamma},\\vec{\\beta}) = \\langle\\vec{\\gamma},\\vec{\\beta}|H_D|\\vec{\\gamma},\\vec{\\beta}\\rangle.\n",
    "\\tag{15}\n",
    "$$\n",
    "\n",
    "要最大化该目标函数等价于最小化 $-F_p$。因此我们定义 $L(\\vec{\\gamma},\\vec{\\beta}) = -F_p(\\vec{\\gamma},\\vec{\\beta})$ 为损失函数，即要最小化的函数，然后利用经典的优化算法寻找最优参数 $\\vec{\\gamma},\\vec{\\beta}$。下面的代码给出了通过量桨构造的损失函数："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 构造损失函数\n",
    "loss_func = ExpecVal(Hamiltonian(H_D_list))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 训练量子神经网络\n",
    "定义好了用于 QAOA 的量子神经网络后，我们使用梯度下降的方法来更新其中的参数，使得式（15）的期望值最大。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-04-30T09:11:12.968989Z",
     "start_time": "2021-04-30T09:11:12.961264Z"
    }
   },
   "outputs": [],
   "source": [
    "depth = 4   # 量子电路的层数\n",
    "ITR = 120   # 训练迭代的次数\n",
    "LR = 0.1    # 基于梯度下降的优化方法的学习率\n",
    "SEED = 1024 #设置全局随机数种子"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "这里，我们在飞桨中优化上面定义的损失函数。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "iter: 10   loss: -2.5212\n",
      "iter: 20   loss: -2.7688\n",
      "iter: 30   loss: -2.9486\n",
      "iter: 40   loss: -2.9832\n",
      "iter: 50   loss: -2.9907\n",
      "iter: 60   loss: -2.9969\n",
      "iter: 70   loss: -2.9990\n",
      "iter: 80   loss: -2.9997\n",
      "iter: 90   loss: -2.9999\n",
      "iter: 100   loss: -3.0000\n",
      "iter: 110   loss: -3.0000\n",
      "iter: 120   loss: -3.0000\n"
     ]
    }
   ],
   "source": [
    "paddle.seed(SEED)\n",
    "\n",
    "cir = circuit_QAOA(n, depth, E, V)\n",
    "# 使用 Adam 优化器\n",
    "opt = paddle.optimizer.Adam(learning_rate=LR, parameters=cir.parameters())\n",
    "\n",
    "for itr in range(1, ITR + 1):\n",
    "    state = cir()\n",
    "    # 计算梯度并优化\n",
    "    loss = -loss_func(state)\n",
    "    loss.backward()\n",
    "    opt.minimize(loss)\n",
    "    opt.clear_grad()\n",
    "    if itr % 10 == 0:\n",
    "        print(\"iter:\", itr, \"  loss:\", \"%.4f\" % loss.numpy())\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 解码量子答案\n",
    "当求得损失函数的最小值以及相对应的一组参数 $\\vec{\\gamma}^*,\\vec{\\beta}^*$ 后，我们的任务还没有完成。为了进一步求得 Max-Cut 问题的近似解，需要从 QAOA 输出的量子态 $|\\vec{\\gamma}^*,\\vec{\\beta}^*\\rangle$ 中解码出经典优化问题的答案。物理上，解码量子态需要对量子态进行测量，然后统计测量结果的概率分布：\n",
    "\n",
    "$$\n",
    "p(z)=|\\langle z|\\vec{\\gamma}^*,\\vec{\\beta}^*\\rangle|^2.\n",
    "\\tag{16}\n",
    "$$\n",
    "\n",
    "通常情况下，某个比特串出现的概率越大，意味着其对应 Max-Cut 问题最优解的可能性越大。\n",
    "\n",
    "Paddle Quantum 提供了查看 QAOA 量子电路输出状态的测量结果概率分布的函数："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-04-30T09:11:55.548009Z",
     "start_time": "2021-04-30T09:11:54.556907Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAYrUlEQVR4nO3de7RkdXnm8e9Dc9EIKgiioVubFRu1vYDaoCuZeAUDEwMYMYIY8RbMCKMskyhOJhjRZFCjEx0xsSMqJsFGjdFObAUTLxmdUbtBFAFbOgjSJMZWvMZBJbzzx66WorrOOVXn1K5Tzf5+1qp1qvbe9dZbLdZTtfdv/3aqCklSd+2x3A1IkpaXQSBJHWcQSFLHGQSS1HEGgSR1nEEgSR2353I3MK4DDzywVq9evdxtSNJu5bLLLvtWVR00bN1uFwSrV69my5Yty92GJO1Wktww1zp3DUlSxxkEktRxBoEkdZxBIEkdZxBIUscZBJLUcQaBJHWcQSBJHbfbnVAmafFWn/3hRT/3+vN+dYKdaJb4i0CSOs4gkKSOMwgkqeMMAknqOINAkjqu1SBIcmySrUm2JTl7yPrnJNmR5Ire7QVt9iNJ2lVrw0eTrADOB44BtgObk2ysqqsHNr24qs5sqw9J0vza/EVwFLCtqq6rqp8AG4ATWnw9SdIitBkEhwA39j3e3ls26GlJvpTk/UlWtdiPJGmI5T5Y/HfA6qp6OPAx4MJhGyU5PcmWJFt27Ngx1QYl6c6uzSC4Cej/hr+yt+xnqurbVfXj3sO3A48aVqiq1lfVuqpad9BBQ6+9LElapDaDYDOwJsmhSfYGTgY29m+Q5L59D48HrmmxH0nSEK2NGqqqW5OcCVwCrADeUVVXJTkX2FJVG4EXJzkeuBW4GXhOW/1IkoZrdfbRqtoEbBpYdk7f/VcAr2izB0nS/Jb7YLEkaZkZBJLUcQaBJHWcQSBJHeelKjVVXipRmj3+IpCkjjMIJKnjDAJJ6jiDQJI6ziCQpI4zCCSp4wwCSeo4g0CSOs4gkKSOMwgkqeMMAknqOINAkjrOIJCkjjMIJKnjDAJJ6jiDQJI6ziCQpI4zCCSp4wwCSeo4g0CSOs4gkKSOMwgkqeMMAknqOINAkjrOIJCkjjMIJKnjWg2CJMcm2ZpkW5Kz59nuaUkqybo2+5Ek7aq1IEiyAjgfOA5YC5ySZO2Q7fYDXgJ8rq1eJElza/MXwVHAtqq6rqp+AmwAThiy3auB1wK3tNiLJGkObQbBIcCNfY+395b9TJJHAquq6sPzFUpyepItSbbs2LFj8p1KUoct28HiJHsAbwR+Z6Ftq2p9Va2rqnUHHXRQ+81JUoe0GQQ3Aav6Hq/sLdtpP+ChwCeTXA88BtjoAWNJmq42g2AzsCbJoUn2Bk4GNu5cWVXfq6oDq2p1Va0GPgscX1VbWuxJkjSgtSCoqluBM4FLgGuA91bVVUnOTXJ8W68rSRrPnm0Wr6pNwKaBZefMse3j2+xFkjScZxZLUscZBJLUcQaBJHWcQSBJHWcQSFLHGQSS1HEGgSR1nEEgSR1nEEhSxxkEktRxBoEkdZxBIEkdZxBIUscZBJLUcQaBJHXcSEGQ5JeS3K13/1lJ3pjk/u22JkmahlF/EfwZ8KMkh9NcbP6fgXe31pUkaWpGDYJbq6qAE4C3VNX5NBeflyTt5ka9VOUPkrwC+E3gl5PsAezVXluSpGkZ9RfBM4AfA8+rqm8AK4HXt9aVJGlqRgqC3of/3wD79BZ9C/jbtpqSJE3PqKOGfgt4P/C23qJDgA+21JMkaYpG3TV0BvBLwPcBqupa4N5tNSVJmp5Rg+DHVfWTnQ+S7AlUOy1JkqZp1CD4VJL/Btw1yTHA+4C/a68tSdK0jBoEZwM7gCuBFwKbgP/eVlOSpOkZ6TyCqroN+IveTZJ0JzJvECR5b1X9RpIrGXJMoKoe3lpnkqSpWOgXwUt6f5/SdiOSpOUx7zGCqvrX3t0XVdUN/TfgRe23J0lq26gHi48Zsuy4hZ6U5NgkW5NsS3L2kPW/neTKJFck+XSStSP2I0makHmDIMl/6R0feGCSL/XdvgZ8aYHnrgDOpwmMtcApQz7oL6qqh1XVEcDrgDcu9o1IkhZnoWMEFwEfAf4HzRDSnX5QVTcv8NyjgG1VdR1Akg0001hfvXODqvp+3/Z3w5PUJGnqFgqCqqrrk5wxuCLJAQuEwSHAjX2PtwOPHlLnDOClwN7AExduWZI0SQsdI7io9/cyYEvv72V9j5esqs6vql8AXs4cJ6klOT3JliRbduzYMYmXlST1zPuLoKqe0vt76CJq3wSs6nu8srdsLhtoLok5rI/1wHqAdevWuftIkiZooRPKHjnf+qq6fJ7Vm4E1SQ6lCYCTgWcO1F/Tm8kU4FeBa5EkTdVCxwjeMM+6Yp59+lV1a5IzgUuAFcA7quqqJOcCW6pqI3BmkqOBnwLfAU4bq3tJ0pIttGvoCUspXlWbaCao6192Tt/9l+zyJEnSVC20a+iJVfXxJL8+bH1VfaCdtiRJ07LQrqHHAR8Hfm3IugIMAknazS20a+iVvb/PnU47kqRpG/Xi9fdK8uYklye5LMmbktyr7eYkSe0bddK5DTRXKHsacFLv/sVtNSVJmp6RrlAG3LeqXt33+DVJntFGQ5Kk6Rr1F8GlSU5Oskfv9hs05wdIknZzCw0f/QHN6KAAZwF/1Vu1B/BD4HfbbE6S1L6FRg3tN61GJEnLY9RjBCTZH1gD3GXnsqr6pzaakiRNz0hBkOQFNBeyXwlcATwG+L94/QBJ2u2NerD4JcCRwA29+YceAXy3raYkSdMzahDcUlW3ACTZp6q+AjywvbYkSdMy6jGC7UnuCXwQ+FiS7wA3tNWUJGl6RgqCqnpq7+4fJvkEcA/go611JUmamnFGDT0S+E805xV8pqp+0lpXkqSpGXXSuXOAC4F7AQcC70wy9ELzkqTdy6i/CE4FDu87YHwezTDS17TUlyRpSkYdNfQv9J1IBuxDc0F6SdJubqG5hv4XzTGB7wFXJflY7/ExwOfbb0+S1LaFdg1t6f29DPjbvuWfbKUbSdLULTTp3IU77yfZGzis93BrVf20zcYkSdMx6lxDj6cZNXQ9zZTUq5Kc5qRzkrT7G3XU0BuAJ1fVVoAkhwHvAR7VVmOSpOkYddTQXjtDAKCqvgrs1U5LkqRpGvUXwWVJ3s7tVyg7ldsPJEuSdmOjBsFvA2cAL+49/t/AW1vpSJI0VQsGQZIVwBer6kHAG9tvSZI0TQseI6iq/wC2JrnfFPqRJE3ZqLuG9qc5s/jzwL/vXFhVx7fSlSRpakYNgj9otQtJ0rJZaK6hu9AcKH4AcCVwQVXdOmrxJMcCbwJWAG+vqvMG1r8UeAFwK7ADeF5VeeUzSZqihY4RXAisowmB42hOLBtJ7yDz+b3nrQVOSbJ2YLMvAOuq6uHA+4HXjVpfkjQZC+0aWltVDwNIcgHjzTh6FLCtqq7rPX8DcAJw9c4NquoTfdt/FnjWGPUlSROw0C+Cn00sN84uoZ5DgBv7Hm/vLZvL84GPDFuR5PQkW5Js2bFjx5htSJLms9AvgsOTfL93P8Bde48DVFXdfRJNJHkWzS6oxw1bX1XrgfUA69atq0m8piSpsdA01CuWUPsmYFXf45UMuapZkqOB3wceV1U/XsLrSZIWYdRJ5xZjM7AmyaG9axmcDGzs3yDJI4C3AcdX1Tdb7EWSNIfWgqB3TOFM4BLgGuC9VXVVknOT7DwR7fXAvsD7klyRZOMc5SRJLRn1hLJFqapNwKaBZef03T+6zdeXJC2szV1DkqTdgEEgSR1nEEhSxxkEktRxBoEkdZxBIEkdZxBIUscZBJLUcQaBJHWcQSBJHWcQSFLHGQSS1HEGgSR1nEEgSR1nEEhSxxkEktRxBoEkdZxBIEkdZxBIUscZBJLUcQaBJHWcQSBJHWcQSFLHGQSS1HEGgSR1nEEgSR1nEEhSxxkEktRxBoEkdZxBIEkd12oQJDk2ydYk25KcPWT9Y5NcnuTWJCe12YskabjWgiDJCuB84DhgLXBKkrUDm30deA5wUVt9SJLmt2eLtY8CtlXVdQBJNgAnAFfv3KCqru+tu63FPiRJ82hz19AhwI19j7f3lkmSZshucbA4yelJtiTZsmPHjuVuR5LuVNoMgpuAVX2PV/aWja2q1lfVuqpad9BBB02kOUlSo80g2AysSXJokr2Bk4GNLb6eJGkRWguCqroVOBO4BLgGeG9VXZXk3CTHAyQ5Msl24OnA25Jc1VY/kqTh2hw1RFVtAjYNLDun7/5mml1GkqRlslscLJYktccgkKSOMwgkqeMMAknqOINAkjrOIJCkjjMIJKnjDAJJ6jiDQJI6ziCQpI4zCCSp4wwCSeo4g0CSOs4gkKSOMwgkqeMMAknqOINAkjrOIJCkjjMIJKnjDAJJ6jiDQJI6ziCQpI4zCCSp4wwCSeo4g0CSOs4gkKSOMwgkqeMMAknqOINAkjrOIJCkjjMIJKnjWg2CJMcm2ZpkW5Kzh6zfJ8nFvfWfS7K6zX4kSbtqLQiSrADOB44D1gKnJFk7sNnzge9U1QOA/wm8tq1+JEnDtfmL4ChgW1VdV1U/ATYAJwxscwJwYe/++4EnJUmLPUmSBuzZYu1DgBv7Hm8HHj3XNlV1a5LvAfcCvtW/UZLTgdN7D3+YZGsrHcOBg689I7UmXW+3rJXxfy/ulu9zGWvNW89//6nUm3Rv/e4/14o2g2Biqmo9sL7t10myparWzVqtSdfrQq1J1+tCrUnX60KtSdebdG+janPX0E3Aqr7HK3vLhm6TZE/gHsC3W+xJkjSgzSDYDKxJcmiSvYGTgY0D22wETuvdPwn4eFVViz1Jkga0tmuot8//TOASYAXwjqq6Ksm5wJaq2ghcAPxlkm3AzTRhsZwmuftp0ruyZrW3Wa016XpdqDXpel2oNel6re8CHyZ+AZekbvPMYknqOINAkjrOIJCkjjMIJKnjdosTytrSm87iKJoznKE5r+HzkxzCmuRBVfWVMZ9zD+DYgb4uqarvTrCvY6rqY4t43kz2luRBNFOW9Pe1saqumWBfz62qd06qnjQrOjtqKMmTgbcC13L7iW4rgQcAL6qqSyf0Ol+vqvuNsf2zgVcClw70dQzwqqp693L0Ncu9JXk5cArNfFbb+/o6GdhQVectR199z5vJ8Ow9byYDdBp9LaG3XwFOHOjtQ1X10Qn2dU5VnTupegu+XoeD4BrguKq6fmD5ocCmqnrwGLXePNcq4LSquvsYtbYCjx78kEiyP/C5qjpsjFqDJ/D19/XEqrrbqLVmubckXwUeUlU/HVi+N3BVVa0Zo9aX5unrsKraZ9RavXozGZ6958xkgE6rr0X29qfAYcC7B3p7NnBtVb1kOfpaqi7vGtqT2/+H7HcTsNeYtZ4L/A7w4yHrThmzVoBh6Xxbb904fhl4FvDDIa9x1Ji1dj5vFnu7Dfh54IaB5fftrRvHwcCvAN8Z0tf/GbMWwO8Dj5orPGk+UEayQHjeaxG9PZ/hAfpG4Cpg5A/cBQL04OXqq4Xe/vOwLzxJLga+CowcBEm+P09fdx2zryXpchC8A9icZAO3z5K6iuZbxwVj1toMfLmqdvmgSPKHY9b6I+DyJJf29XU/mm+Qrx6z1meBH1XVp4b0tZgZXGe1t7OAf0xy7UBfDwDOHLPW3wP7VtUVQ/r65Ji1YHbDc2cPsxigk+xr0r3dkuTIqto8sPxI4JYxa30XOLKq/m1wRZIbd928PZ3dNQTQu1DO8ey6H/LqMescANxSVT+aUF/70/yHO7hPefA/5Kmb1d6S7MGuB/43V9V/LF9XkOQ04ByaXUO7hGdVvWuMWh8BXldVnxiy7p+q6rFj9nYs8Baa42S7BOg4+7yTXAC8s6o+PWTdRVX1zOXoq4XeHgn8GbAft+9RWAV8Dzijqi4bo9ZraD5vPj9k3Wur6uWj1lqqTgfBTr0Pcqrq5lmqNauSHEzfB+6wbzTLUWuO+vtW1eA36KnWmtXwhJkO0Jnsa6ck9+GO/91+Yzn7WarOBkGS+wGvA55Ik+YB7g58HDh78CDyiLWeRPNzb9G1FnidK6vqYctVK8kRwJ/TTBe+neZ9rqR5zy+qqsvHqPUImm9W9+COB1HHrrXA60zsoNtSau1O4dl7jWUN0EkP7Z7VoeLTqDWKLh8juBj4U+DUnd8y0lxn+ek0oxUesxy1kvz6XKuA+4zR00Rr9bwLeGFVfW7gdR4DvBM4fIxa75xUrSQvnWsVsO8YPU20Vq/eEQwJzyTfZULhuZhaI7iaZnfM1GvNN7Q7ydhDuyddbx6XMrl/s0nWWlCXg+DAqrq4f0HvQ3xDknEPfE6y1sXAXzP8AONdlrEWwN0GP7gBquqzScYaijrhWn8MvB64dci6cc+en2QtmNHw7D1vVgP0TcDRcw3tBkYe2j3pegsMFb/nOE1NstZSdTkILkvyVuBC7jhq6DTgC8tY60vAn1TVlwdXJDl6GWsBfCTJh2mGPPa/z2cD455MM8lalwMfHHagLskLlrEWzG54wuwG6CSHdk+63iSHik+y1pJ0OQieTTNe+VUMjBpi/OGjk6x1FjDX+OKnLmMtqurFSY5j1zM+z6+qTctVi+b/UHNd4nTc679OshbMbnjC7AboJId2T7reJIeKT7LWknT2YLE0LXME3sZFBN6kaz0QuLmqdgxZd/A4B6EnWav3nAcz/H2ONbR70vUmOVR80sPOl9RLV4MgyZ403+JPZGDOEOCCwbMal6HWU2lOqpmJWiO81vqqOt1a0u6ny0HwHpqhihdyxzlDTgMOqKpnWGuXegfMtQr4YlWttNYu9e4BvILm2+jBNAfuv0kTxucNTj0xrVoD9U4E7j2h3pZca4HX+UhVHTeJWpOuN6u1RtHlYwSPql3nDNkOfDbNJGbW2tUOmtP++6dGqN7je1trqPfSnE/yhJ0nHfVORnpOb92Tl6lWf73HD9Q7bQm9LblW7+zdoauAI8boaeL1ZrXWUnU5CG5O8nTgb6rqNvjZ2YxPZ9c5SazVuA54UlV9fXBFxp8bpQu1AFZX1Wv7F/Q+KM9L8txlrDVfvdcmed4y1toMfIo7hvFO9xyz1qTrzWqtpamqTt6A1TTj7L9JM2vgV3v3LwYOtdbQemcAh8+x7r9aa+hzLgVeBhzct+xg4OXAPyxXrVnuDfgysGaOdTcu4n1OrN6s1lrqrbPHCGDOkQQfqkVc/KILtXr1JnbBkI7U2h84u1dv566lf6MZWnxejTHf0CRrzXJvSU4CrqyqXWahTXJiVX1w1FqTrjertZaqs9csTnPxi4to9v9+rncDeE+Ss601tN7LaKbMCPD53i2L7O1OXwugqr5TVS+vqgdV1QG924OrmVnyxOWqNcu9VdX7h3049uw/Tq1J15vVWks2zZ8fs3Sj2U2y15Dle9Ncachau0lvs1prhNf6+izWmuXefJ/t3Lp8sHiSF7/oQq1Z7m1Wa5EJXh1rkrUmXW9Wa0263qzWWqouB8FZTO6qVl2oNcu9zWotmOzVsSZ9Gc1Z7c33ubj3uWidDYKq+miSw5jAxS+6UGuWe5vVWj2TvPTlpC+jOau9+T4X9z4XrdOjhiRJHR41JElqGASS1HEGge7UkqxM8qEk1ya5LslbkuwzwvOGXmM3ybnpXdQnyVlJfm6O7Z6S5AtJvpjk6iQv7C0/McnaEV5/pO2kSTAIdKeVJMAHaC6YsgZYA9wVeN1ia1bVOVX1D72HZwG7BEGSvYD1wK9V1eHAI4BP9lafCIzyAT/qdtKSebBYd1pJngS8sqoe27fs7jTnCKwCTgLWVdWZvXV/T3Npz0/2fhH8Bc2smd8ATq6qHUneRTPa4+eBPwG2At+qqif0vcYBwFeA+1fV/+tb/ou9536vd3sa8ETgdJoT1rYBv0kz8+TgdgDnAwcBPwJ+q6q+MpF/KHWevwh0Z/YQ4A6XTqyq7wPX05wXMJ+7AVuq6iE0M0S+cqDOm4F/oZkS+gkD626mmWPnhiTvSXJqkj2quSThRuD3quqIqvpn4ANVdWTvl8M1wPPn2G49zaR3jwJ+F3jr2P8a0hw6ex6BtIDbaGZpBfgrml1MI6uqFyR5GHA0zQf3MTTXDRj00CSvoZl2eF/gksENkuwL/CLwvmZvFwALHueQRmUQ6M7saprdPz/T2zV0H5pdOg/ljr+K7zJPrbH3oVbVlcCVSf4S+BrDg+BdwIlV9cUkzwEeP2SbPYDvVtUR4/YgjcJdQ7oz+0fg55I8GyDJCuANwFt6++6vB45IskeSVTRnE++0B7eHyDOBTw+p/wNgv8GFSfZN8vi+RUdw+9xFg8/ZD/jX3gHmU4fV7u3O+lqaiwyRxuHzvXFpHAaB7rSqGQnxVOCk3txB3wZuq6o/6m3yGZpv6lcDbwYu73v6vwNHJfkyzQHdc4e8xHrgo0k+MbA8wMuSbE1yBfAqbv81sAH4vd7Q0l8A/oBmevDP0BxgZo7tTgWen+SLwFU08/5LE+GoIXVGb9TOe4CnVtXlC20vdYVBIEkd564hSeo4g0CSOs4gkKSOMwgkqeMMAknqOINAkjrOIJCkjvv/uYzWl+Mg94MAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "state = cir()\n",
    "# 模拟重复测量电路输出态 1024 次\n",
    "prob_measure = state.measure(shots=1024, plot=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "通过测量，找到出现几率最高的比特串。 此时，记其在该比特串中对应的比特取值为 $0$ 的顶点属于集合 $S_0$ 以及对应比特取值为 $1$ 的顶点属于集合 $S_1$，这两个顶点集合之间存在的边就是该图的一个可能的最大割方案。\n",
    "\n",
    "下面的代码选取测量结果中出现几率最大的比特串，然后将其映射回经典解，并且画出对应的最大割方案：\n",
    "- 红色顶点属于集合 $S_0$，\n",
    "- 蓝色顶点属于集合 $S_1$，\n",
    "- 虚线表示被割的边。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-04-30T09:11:55.906817Z",
     "start_time": "2021-04-30T09:11:55.625551Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "找到的割的比特串形式： 1010\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAb4AAAEuCAYAAADx63eqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAArx0lEQVR4nO3deXQUVdoG8Kc7a3cIaXaSsIgLOgQBBUQjED4PokACzjCAyh5l0wERUZxBgRnHT0FAFkc2ZVdBGBeCII6CmASH7QPZJCECYiAKBhLJ0tn6fn9cu0lC9qruqu56fufknHQnqXrFpJ++dd+61ySEECAiIjIIs9YFEBEReRKDj4iIDIXBR0REhsLgIyIiQ2HwERGRoTD4iIjIUBh8RERkKAw+IiIyFAYfEREZCoOPiIgMhcFHRESGwuAjIiJDYfAREZGhMPiIiMhQGHxERGQoDD4iIjIUBh8RERkKg4+IiAyFwUdERIbC4CMiIkNh8BERkaEw+IiIyFAYfEREZCgMPiIiMhQGHxERGQqDj4iIDIXBR0REhsLgIyIiQ2HwERGRoTD4iIjIUBh8RERkKAw+IiIyFAYfEREZir/WBRAZgt0OHD0KHDoEXLgAXLsGFBcDISFAgwZAx45A585As2ZaV0rk8xh8RO5y5Ajw1lvArl3ATz8BVqsMu7y8st/n7y8D0G6X33PXXUB8PPDnPwNBQZqUTuTLTEIIoXURRD7Dbge2bAHmzAHOnAEKCoCSktofp149wGQCxo8H/vIXoHVr9WslMigGH5Fatm0DRo0CCguBnBx1jhkYCJjNwOjRwLx5cmRIRIow+IiUunoVGDcO2L79xsuYarFYAJsN2LQJ6NHDPecgMggGH5ESX30FDB4sA6+gwP3ns1jk/N+iRYCfn/vPR+SDGHxEdbVlCzByJJCf79nzWq1Ar17ARx+x+YWoDhh8RHXx4Ydy3s3ToedksQDR0cCOHUBAgDY1EHkpBh9RbX3xBfDII9qFnpPFAvTtK0eeJpO2tRB5Ea7cQlQbly8DQ4ZoH3qArGHnTmDlSq0rIfIqHPER1UZsLPCf/8hbFvQiJAQ4eRJo1UrrSoi8Akd8RDW1eTOwe7e+Qg+Q3aSPPQbwPSxRjTD4iGri2jVg7Fj33aenRHEx8N13wLp1WldC5BUYfEQ1sW6dDBi9ys0FXnmFoz6iGuAcH1F1hADatAF+/FHrSqoWEiLnH++7T+tKiHSNIz6i6iQmAr/+qvw4vXsDCxYAX38NZGfLQHV+xMQoP35enlzPk4iqxG2JiKqzYIE6c3tPPy3v/3MXIeR6oVeuAA0buu88RF6OIz6i6iQmqjN3JoTchDYhAXjvPeXHq0hgILB/v3uOTeQjGHxEVbl0STaOqGHYMKBFC2DAAPfddJ6XBxw44J5jE/kIBh9RVQ4dAoKD1TmWJ1Z7KS4G9uxx/3mIvBiDj6gqBw6oN+LzlCNHtK6ASNcYfERVSU3V9/17FcnMBBwOrasg0i0GH1FVvG20B8gNaj2xKS6Rl2LwEVWlqEjrCmrPbPbOuok8hMFHVBWLResKaq+42DvrJvIQBh9RVRo00LqC2jObAX+uTUFUGf51EFWlc2d5s7kaK7cMGQJ07So/b9my7NcmTpR7/QGyk/TDD+t+njZtuCM7URUYfERV6dxZvdFT377A6NEVf23o0Oufr1mjLPi4SDVRlXipk6gqd97pmRvP1WK1At27a10Fka5xWyKi6tx2G5CWpnUVNRMaKneJ79xZ60qIdIsjPqLqDB+u3rJl7hYUBHTqpHUVRLrG4COqzvjx3rGzucUCTJkib2Anokox+Iiq07w58OCD+u+UdDiAsWO1roJI9xh8RDXx/POycUSvzGagXz+gaVOtKyHSPTa3ENWEEBA9esDx7bfw0+MC0BYLcPgwcPvtWldCpHsc8RHVQOrp0/iT3Y58PYZeSAiKZ8zA5qNHwfexRNVj8BFVoaSkBAsWLEDHjh3xyaFDmFW/Por11OFpNgO33ILZOTkYMmQIYmNjceHCBa2rItI1Bh9RJVJSUtCzZ08899xzsNvtGDVqFF46exb+99wjbxvQg+Bg4MMP8Yf27WGz2bB9+3ZERUVhzZo1HP0RVYJzfEQVKCgoQJs2bZCRkYGIiAisWLEC/fv3l1/MzgbuuQc4e1bb7X8sFuCjj4CHHwYAXLx4EePHj8e2bdsAAP369cOKFSsQGRmpXY1EOsTgI6rE2rVrsXv3brz55ptoUH6XhkuXgHvvBdLTtQk/iwVYuxYYPLjM00IIbNiwAZMnT0ZWVhbCwsJw8OBB3HrrrZ6vkUinGHxEkHN5ixYtQkBAACZNmlSzH7pyBYiJAX74wXPreZpMMvQ2b5a3L1TCOforKSnBZ599BpPe70Ek8iAGHxleamoqxowZg7179yI4OBhnzpxBeHh4zX44Px944QXg3XfdH35Wq9zOaPNmuXh2NYQQyMvLQ0hICAA5Z7l3716MHj2aQUiGxuYWMqzSHZt79+5FeHg4Nm/eXPPQA+Toa8kS4MsvgchI9+x87hzlvfACcOxYjUJP/pjJFXolJSWIj49HfHw8Oz/J8Bh8ZEipqallOjZHjhyJEydOINa5GWxtRUcDp08DkyfLHRLq1VNeZGCg7Nrs1QvYtw+YNQsICKjTocxmMyZMmMDOTyLwUicZVO/evfHVV18hPDwcK1asqHvgVaSgANiyBZgzR87/ATXfwd3fX17SFAKYMAF4+mmgdWvVSmPnJxGDjwxECOGa20pNTcXcuXPxxhtv3NixqaaTJ4HkZCApCfj2W+DcOcDfHzn5+RAAQq1W2RUaFgbcdZcc3XXpIptm3HSvYPnOz+bNm+PMmTOwuOMyLZEOMfjI5zk7NpOTk7FlyxZtGzsKC4HMTNwcEQE/AKfPnwfq15fB52HO0V+vXr3w3HPPefz8RFph8JFPK92xCQBff/01YmJiNK4KrvDV+s9PCAEhBMxmOd3/wQcfwG63s/OTfBqbW8gnVdSxmZCQoIvQ0xOTyeQKvcuXL2PixIns/CSfx+Ajn6N6x6ZBNG7cGEuWLCnT+bl69WrNR6VEamPwkc/ZsGFDmVHe2rVr3dvA4iNMJhNGjBiBkydPIi4uDtnZ2Rz9kU/iHB/5hMLCQgQGBgKQC0z/4x//wLRp03QbeHqZ46tM+c7P+++/H0lJSVqXRaQKBh95NWfH5ltvvYWDBw+iYcOGWpdUI3oPPqeMjAw89dRTmDFjBrp06aJ1OUSqYPCR1yrfsbl8+XKMGzdO46pqxluCryKTJk1C586dMWrUKHZ+kldi8JHXcY7yZsyYAbvd7p7VV9zMW4Pv22+/RXR0NACu+kLei80t5FXYsamte++9F+vWreOan+TVGHzkVdLT09mxqSFn56fzzUZ2djbGjBnDzk/yKrzUSbqXmZmJRo0auR6vX78esbGxXh143nqps7TynZ+TJ0/GokWLtC6LqFoMPtIt51zeyy+/jB07dqBnz55al6QaXwg+p4yMDMyaNQvz5s1D/fr1AZRdEJxIb3ipk3Sp9FxeXl4edu7cqXVJVAlnc5Ez9PLy8tCtWzfO/ZFuMfhIVypbY/PVV1/VujSqoffeew8HDhzg3B/pFi91km6cO3cOw4YNc92XN3LkSCxcuNCr5/Iq40uXOssrP/cXFhaGhQsX8r4/0g0GH+nG5cuXERUVBX9/f6+7L6+2fDn4nLjbO+kVg480lZaWhtatWyMgIAAAsG/fPrRt29YnR3mlGSH4gBtHf5s2bcKQIUO0LosMjnN8pAnnXN6dd96J1157zfV8t27dfD70jKT0fX9z587F4MGDXV/Lz8/XsDIyMgYfeVz51Vd++uknnx/5GF1ERASef/5510j36NGjaN26Nff7I00w+MhjKuvYXLlyJZseDOa9997D5cuXud8faYJzfOQRV69eRWxsrCE6NmvCKHN8lamo8/PNN9/E6NGj+SaI3I4jPvIIm80Gi8XCNTYJQNW7vV+8eFHr8sjHccRHbpOamorAwEDcdNNNAGR7u8ViYeCBI77SSo/+ioqKcOzYMbRp00brssiHMfhIdaX3y+vWrRt27doFs5kXF0pj8N0oIyMDR44cQd++fQEADocDly5dQvPmzTWujHwNX41IVeU7Nlu3bg273a51WeQFwsPDXaEHAEuXLsUdd9zBzk9SHYOPVFFZx+batWthtVq1Lo+8UHJycpm5P3Z+klp4qZMUczgc6N27N3bv3g2AHZs1wUud1WPnJ7kLR3ykmNlsRu/evdmxSaqqrPOzf//+uHLlitblkRfjiI/qJDU1FefOnUOfPn0AAMXFxbh27RoDr4Y44qud0qO/m2++Gf/9739d67sS1RaDj2qldMem1WrFiRMn2HVXBwy+usnIyEBOTg5uu+02AMClS5dQVFTEHR+oVnipk2qsfMdmbGwsgoKCtC6LDCQ8PNwVekIITJw4EVFRUdztnWqFwUfVqqpjk5c2SSt2ux2FhYXIzs7mbu9UKww+qlZ8fLxrlDdy5EicOHHCpzeJJe9gsViwdetWrFu3DjabDdu3b0dUVBTv+6NqcY6PqrV3714MHToUS5cuZeCphHN86iq/23v//v3x6aefws/PT+PKSI844qMbpKSkYN68ea7H0dHRSEtLY+iRbkVERJQZ/bVt25ahR5XiiI9cSkpKsHDhQrz00kuw2+3YuXOn63YFUhdHfO6TkZGBsLAw14pBBw4cQEREBDs/yYUjPgIgR3k9e/bEtGnTYLfbMWrUKHTt2lXrsohqLTw83BV6v/32GwYNGsTOTyqDwWdwJSUlmD9/Pjp16oS9e/ciIiIC27Ztw5o1a9ixSV6voKAAHTt2ZOcnlcHgM7jXX3+9zCjv+PHj6N+/v9ZlEamiSZMmFXZ+cvRnbJzjM7isrCz06dMHs2bNYuB5EOf4PK985+cTTzyBd955R+OqSAsc8RlMamoqRo8e7dojz2azYd++fQw98nnlOz//+Mc/al0SaYTBZxAlJSV488030bFjR6xduxZz5sxxfY1bvJBROHd8OHv2bJk3e+vXr+fcn4Ew+AwgNTUVMTExmDp1qmv1lcmTJ2tdFpFmbDab6/P9+/dj9OjRnPszEAafDys9yktOTkZ4eDi2bt3KNTaJSmnZsiX69+/Pzk8DYXOLD9u2bRvi4uIAcFd0vWFzi75UtNv7woULMWrUKE4F+CAGnw8TQmD8+PGIi4tzBSDpA4NPn8p3fr7yyit46aWXNK6K1MbgU0FxMZCfD5hMgMUCaLVEYGpqKiZOnIglS5agXbt22hRBNcLg0y/n6G/27NlITExERESExvUAdjtQWAgEBckPDkKVYfDV0unTwP79wLffAklJwKlT8hfSGXbFxYDVCkRFAT17At26yY9WrdxXU+ld0e12OwYOHIhPPvnEfSckxRh8+ldcXAx/f38A8m/sxRdfxDPPPIMWLVq47Zy5ucCBA8DBg8DXXwOHDgG//AKYzfLD4ZAfkZFA165ATAzQuTPQpQsQHOy2snyPoGrl5QmxerUQd9whhMUiRL16Qsj3YVV/mM1ChIYKERwsRNeuQmzZIkRhobq1paSkiOjoaAFAABAjR44UV65cUfckpDrn/y/yDgsWLBAARFhYmFi1apVwOByqHv/ECSHGjpWvL/XrCxEYWLPXmKAg+f316gnxzDNC/PCDqmX5LP7lVeHKFSEmT5a/VDUNu6o+QkOFsNmEmDlThqkSxcXFYv78+SI4OFgAEOHh4SIhIUGd/3ByOwafd7l48aKIi4tz/X/r16+fSE9PV3zcnTuF6NxZBp6/v7LXl4AA+Sa7Rw8hkpJU+I/2YfzLq8TWrUI0aCDfUSkNvPIfFosQkZFCJCfXvb6zZ88Ki8XCUZ6XYvB5H4fDIdatWydsNpvi0d/Vq0IMHSqE1ar+64vzNWbCBCFyc9X/d/AFnOMr5+pVYNw4YPt2IC/PveeyWID4eOCNN+Tn1SkpKYHZbHbND61evRpNmjThBrFeiHN83isjIwPjx49HQkICAGDLli0YNGhQjX9++3ZgxAg5n1dQ4K4q5WuKzQZs2gT06OG+83gjBl8paWnA/fcD2dnu/YUszWIBWrYEEhOBpk0r/77U1FSMGTMGY8aMwZNPPumZ4shtGHzeTfze+bl582Z8/PHHNdrtXQjg5ZeBN990/5vq0iwW4PXXAS7WdB2D73fHj8t3RdnZ8hfUkwICgGbNgP/+V3ZrlVa+Y/OOO+7A8ePHa/SHRvrF4PM96enpmDZtGubPn3/Dbu9CAH/5C7BmjWdDz8lqBaZPB2bO9Py59YhLlgH44Qd560FWludDDwCKioCMDCA6Grh8+frzqamp6NmzJ5577jnXGpt79+5l6BHp0LRp07Bp06YK1/ycPl270APkeefMAebN0+b8uqPR3KJuXLsmRESEvPXAHZPMte3K6tBBCLudHZu+Dmxu8TkXLlwQsbGxN3R+rljhviaW2n5YrUJ8/LHW/1LaM/ylzvh44IMP5MoIehASAkyfXowtWzrj6NGjXGPTR/FSp28Souyan6Gh7VFQ8H8oLAzQujQXm032MzRqpHUl2jF08O3aBcTFaXf5oTIWC7Bhw/cIDPyBHZs+isHn2zIyMjBu3Hhs2/YigK4A9BN8gYFA376AkRd3MmzwXbsG3HorcOmS1pXcyGQC2rUDjhwBfl8xiXwMg8/3LV0q8OyzxSgo0E/oOVmtwIYNgFE3oTdsc8uSJTL89EgI4Nw5YMsWrSshorqw24Hp0026DD1AXuV6+mm57qcRGTL4SkqAhQvljgp6lZsLzJ2rdRVEVBdbtmjTIV4b164BX36pdRXaMGTwff65fppZqnLqlLy/kIi8y5w5QE6O1lVULSfHuG+uDRl8c+cqv8xpswEjRwLvvCO3DrlwQa72kp0ttxT5+98BpY2YhYVylQci8h6HDwNnzig/TsOGwKuvAkePytera9fk56++qvy1xSk5GTh/Xp1jeRPDNbdcuybbeIuKlB1n0KDq5+AuXgQeeABISan7eerVA377jRtP+ho2t/iu6dOB+fPllEpdRUUBX3wBVLYH7sWLQJ8+wIkTdT8HIPfw+9//BZ59VtlxvI3hRnyHD8uOJrVkZQEbN8o1+ObNk7+QThERwPLlyo5fUmLMd2RE3mrPHmWhFxwMfPTR9dC7elVepZozR34OyK/9+99yN3Yl7HZZr9EYrln+0CF15veuXAGeeQZYubJsk8ycOXJerlkz+bhHDzlqq+v1/oAAWXPr1sprJiL3EkL5vPywYUDbttcfP/647EsAZEht3y4/v/12+b2rVik738GDyn7eGxluxLdnjzo7L+zeDSxefGNn6K+/yp0WnMxmecNoXeXkyMWriUj/zp9XfotA6XvrsrOvhx4gP//tt+uP//QnZecC5L3Mer21y10MF3yHD7v/HHfccf3ztDQ5OqwrhwP49lvlNRGR+333nbxKo0SnTtc/P3u27NeEKPtcx47KzgXIqR+jdY8bLvjc3WI8cybQvn3Zx0plZys/BhG5X3a28hFf6TU0S4/uKnqucWNl56rqPL7McHN8hYXuOa7JJJtbpk69/tzs2XIBbKU8tSkuESljt6u7GkpF3dxqd3gL4R33NavJcMFndsMYt149GXDO9aQdDuCFF2RLsxq4/R6Rd/DzUx5MmZnXN6SuX//Gr5d+7tdflZ0LkPUa7TXGcMGntP23vFatgIQEoEMH+Tg3FxgxAvj4Y/XOERys3rGIyH0sFuVvro8cuR58bdrIYHLe7mkyyeecvvtO2bmcLBZ1juMtDDfH17y5esfq1g3Yt+966P30E9C9u7qhB8hwJSL9i4xUPuIrvV1Q/fpyCyGnvn2B0NDrj9V4rSkpqfxGeV9luJVbnnoKWLpU+XHuuw/46qvr75SKi+XC1z//fOP3btoEpKfX7TxBQcBrrxlvZQVfx5VbfJMaK0MFB8uRnPNevqys6wthjBt3fbmy1FTZ1al0fi4oSF6pMtLlTsNd6rzvPmD9euXdnW3blr084O8PTJtW8fcePFj34AsOBjp3rtvPEpFnhYYCTZqUXcGptux2eX/ef/4DhIfLdYGnTy/7PRkZ8nvUaEpp29ZYoQcY8FKnt4VIXh5w111aV0FENdWli/JjnDghb4t67TX5eW6u/DhxQj7Xvr3ydTqdundX5zjexHCXOktK5DsovW8Z4nTLLfImePItvNTpu5YulVd/8vK0rqR6oaFyybM//1nrSjzLcCM+Pz/giSeUr67gCSEhnNsj8jaPPeY9O5ubTEBcnNZVeJ7hgg8AJk3yjmvaDoe8NYKIvIfNJrct0/trTGAgMH68+rd4eQNDBt8ttwBdu2pdRdX8/YFHH634BlYi0repU/UfKGYz8PTTWlehDUMGHwDMmCEvJepVQADw3HNaV0FEdXH33XLbIJNJn3O4/v5Ar17G3e7MsMH30ENA797KtgxynzzExJxEu3b6/KMhouq9+uqPMJn0udBucLDcS9SoDBt8APDOO3pcqscB4BI+//wuDBgwABeV3BBERB7ncDiwaNEiDBr0Bzgc/wCQq3VJZYSEyMU2WrTQuhLtGDr4GjcGVq+W+1HphcViwssvH0dYmAXbtm1DVFQU1q9fz7Z3Ii+QlpaGXr16YcqUKcjPz8djj6XjjjuCVN9Roa78/eW9zPHxWleiMUFizBghrFYh5FKw2n1YrUK88oqsKT09XfTr108AEABEbGysyMzM1PYfilTj/P9KvuPtt98WFotFABDNmjUTn3zyiRBCiJQUIerX1/71xWQSonFjIdLTNf6H0gFDj/icVq4EHnxQ28ueVqu8v3DGDPk4MjIS27Ztw+rVqxEWFoZz584hRM/dOEQGl5eXh/z8fAwbNgwnT57EwIEDAcglwXbt0raZzmQCwsKA5OTrOz8YmeFWbqlMUZG8feDzzz2/4oLVCkyYIDeyreiSyIULF5CdnY127doBAK5cuYKCggKEh4d7tlBSDVdu8X4OhwOnTp1y/V2WlJRgz549eOCBByr8/v375RtsT+92bjbLewsTE4HfSzU8jvh+FxAAbN4MPPmk50Z+JpM81+zZctPayuYBIiMjXX9cADBp0iS0a9eOc39EGklLS0NMTAyio6Nx4cIFAICfn1+loQcA99xzfcTlqdcYqxW4+Wbg0CGGXmkMvlLMZmDRIjnqCw937y+n1SovgezbBzz/fM1/rqCgAFevXkVWVhZGjhyJgQMHsvOTyEOcHZsdOnRAUlISgoODce7cuRr/fPv2cjuh+Hj3vr4431RPnQqcPAncdJP7zuWVtJ1i1K/cXCEmTBDCYhEiIEC9CeagIHnM2bOFKCysW20Oh0OsWrVKhIWFCQDCZrOJdevWCYfDoe4/ArkN2NzidU6fPi26d+/u+n83bNgwRQ1nyclCREYKUa+euk0s9eoJcdttQhw5ouJ/vI/hX141UlKEeOopIUJC5EddfxlDQ4Ww2YT429+E+OkndWr76aefRN++fV1/iHFxcaKoqEidg5NbMfi8y/vvv19hx6ZSBQVCbNwoxF13yTfEfn51e30JCBAiOFiI6GghPv1UCL4MVI2XOqvRti3wr38Bly7Jy6B33y3X4AsJkVt6VDQvZzbLNTatVvlx//3yZvlLl4BXX1XvxtEWLVrgs88+w6pVqxAWFoYWLVrA399wewsTud1tt92GwsLCGzo2lQoMBIYOBf7v/2Tzy6hRcpolIEB2YVa23mdwsHyN8fcHWraUi00fPSrnEAcMkM9T5djVWQcOB3D6tJww3r8f+OUXuUmkySQDsUULOZHdubO8tu6Jm1fT09Nhs9lQr149AMCRI0fQtGlTREREuP/kVGvs6tQ3h8OBL774Ag8//LDruZSUFNx+++0eOX92NnD4sHyNOXpU7h9aUHA98Dp2lK8vnToBv//JUy0w+HxQbm4uOnTogCtXrmDx4sUYPny464WW9IHBp19paWmIj49HYmIitm7dijgjbljn43ip0wfl5eXh9ttvL9P5mZGRoXVZRLpWumMzMTERzZo149SBj2Lw+aAmTZqUmftLSEjgfX9EVSi/xqZzLq9v375al0ZuwEudPi49PR3jxo3Djh07AADx8fF49913Na6KeKlTP3bt2oXY2Fjk5+ejWbNmWL58uWrNK6RPHPH5uPKdn/369dO6JCJd6dKlC5o0aaJ6xybpF0d8BpKZmYlGjRq5Hn/44Yfo3r07Oz81wBGfdhwOB1atWoXHH38c1t/3JCv/t0G+jSM+Ayn9h33kyBEMGzaM+/2RoTjn8saOHYsZzq1QAIaewTD4DKpJkyZ48MEH2flJhlBRx2avXr20Los0wuAzqMjISHZ+kiFU1rHJuTzj4hwf3dD5+fe//x0zZ87UuCrfxjk+z/jhhx9w5513smOTymDwEQD5ArxmzRrMmjULiYmJaN26tdYl+TQGn+cMHjwYQUFBWLx4MRo2bKh1OaQDDD4qo7CwEIGBgQDkvMisWbMwceJEdn6qjMHnHg6HA2+99RZiYmLQsWNHAGV/p4kAzvFROaVfIJYtW4Z//vOf7Pwkr+Ccy3vmmWcwZswYlJSUAABDj27A4KNKDRw4EP369XN1fg4YMIC7vZPuVNSxOWvWLPj5+WldGukUL3VSlYQQWLt2LaZMmYLs7GzYbDbu+KACXupUR+mdFABg2LBhnMujajH4qEYuXLiAcePGYfv27QCAjRs3YujQoRpX5b0YfMrl5+fjpptuwqVLl9ixSbXC4KMac47+3n//fWzfvp1btijA4FPH0qVLkZyczFEe1QqDj2pNCOF64f75558xbdo0zJ07l52ftcDgqz2Hw4ElS5bAarVi7NixAMr+LhLVFIOPFBkxYgQ2bNjAub9aYvDVTlpaGsaMGYOkpCSEhITgzJkzaNq0qdZlkZdiVycp8vrrr7Pzk9ymdMdmUlISmjdvjvfff5+hR4pwxEeKsfOz9jjiq17pUR4ADB8+HIsWLeJcHinG4CPVlO78NJvNOHbsGNq1a6d1WbrE4Ktez549kZiYiObNm2P58uUYMGCA1iWRj2Dwkaqco7/z589zoesqMPiqd+zYMcyfPx8LFizgKI9UxeAjt/v888/x9ttvY9myZez8/B2DryznGpsHDhzAunXreImc3IrBR24lhECHDh1w/Phxzv2VwuC7rvzqK8nJyYiOjta4KvJl7OoktzKZTNixYwf69u3Lzk8qo6I1Nj/++GOGHrkdR3zkEc79/p599ll2foIjvvKjvMcffxyLFy9Go0aNNK6MjIDBRx5Verf3yMhIfP/99wgNDdW6LI8zevC98MILeOONN9CsWTMsW7YMjzzyiNYlkYEw+MjjnJ2f4eHheOihhwAAxcXF8PPzM8zoz4jBV1xc7FrfNS8vDzNnzsRf//pXjvLI4xh8pAsvvPACTp06heXLlyM8PFzrctzOSMHnXGNz6dKl2LdvH8LCwrQuiQyOzS2kuaysLLz77rtISEhAu3btuNu7D3Huij5lyhSkpKRgy5YtWpdExOAj7dlsNnz33XdlOj8HDhyIjIwMrUujOqqoY/OTTz7BE088oXVpRLzUSfphpM5PX77UyV3RSe844iPdMJlMGDNmDI4fP+4a/fHSmPdJSUkpM8rbsGEDQ490hSM+0iUhBNatW4c+ffq4ml2ysrIQFhbmE6M/XxvxZWVlwWazuR6vXLkSgwYNYuCRLjH4yCsUFxcjOjratVK/t3d++krwOTs2X375ZXz55Ze45557tC6JqFq81Ele4eTJk0hNTWXnp46U7ti8du0aEhIStC6JqEYYfOQVnAtds/NTe5V1bL7yyital0ZUI7zUSV6los7Pt99+G4899pjWpdWKt17q/PHHHzFixAh2bJJX44iPvEpFnZ+XL1/WuizDCAgIwLFjx9ixSV6NIz7yWkIIJCQkIDY2FmazfA+XlpaGW265Rfedn9404jt79ixatWoFPz8/AMA333yD9u3bM/DIa3HER17LZDJhwIABrtA7f/487r77bu73pxLnXF5UVBQWLFjger5nz54MPfJqDD7yGadOnYLZbMa2bdsQFRXFzk8F0tLSEBMTgylTpiA/Px+pqalal0SkGgYf+Yw+ffrgxIkT6NevH3d7r6PSHZtJSUlo3rw5Pv30U6xcuVLr0ohUwzk+8jnO/f6mTJni6vxcv349YmNjtS7NRY9zfJmZmXjkkUeQlJQEABg+fDgWLVrEy5rkczjiI59jMpkwevRo1+gvJycHLVq00Los3bPZbBBCuEZ569evZ+iRT+KIj3yaEAJHjx5Fx44dXc9988036NGjh6adn3oZ8aWlpcFisSAyMhKAbBCqV68eA498Gkd85NNMJlOZ0Pvoo48QExPj+VVfCgqAQ4eAFSuA2bPxBoCFAPC3vwFz5wI7dwK//uqxckrP5T355JOuAG7VqhVDj3yev9YFEHlSUVERwsLCkJCQgMTERPfu93fsGPCvfwFffQX8+CMQHAyUlAB5eZjm/J7XXgMCAgCLBbDbgdBQ4O67gSeeAP74RyAwUPWyyu+X16hRI9jtdlgsFtXPRaRHvNRJhpOeno5x48Zhx44dAIDY2FgsX74cERERyg9eWAh89BEwZw6QkiIfl5TU/jihoYDZDEyYADz9NNCypeLSnDsp/PWvf0V+fj6aNWuG5cuXY+DAgYqPTeRNGHxkSBV1fu7YsQP33ntv3Q+6fTswcqQMu2vX1Ck0KAgwmeQIcO5cwGqt02FKSkrw4IMPYvfu3QDYsUnGxjk+MqTynZ8NGzbEnXfeWbeDZWUBjz4KDB4MZGaqF3qAnBu024FVq4C2bYHk5Dodxs/PD/fffz87NonAER8RhBD45Zdf0Lx5cwBATk4OPvvsMwwZMqT6ub+vvwYGDQJyc2VIuZvFAowdCyxYAPy+dmZl0tLScOHCBcTExAAACgsLkZOTw8Ajw+OIjwzPZDK5Qg8AXnzxRTz66KPVd35+/DHQrx9w5YpnQg8A8vOBd96RjS+FhRV+S+mOzaFDhyIzMxMAEBgYyNAjAoOP6AZdu3Z1dX5Wutv7v/8NDBsmg8jT8vKAL78E4uKA4uIyXyq9K3p+fj569+6t+50qiDyNlzqJKlC+8zMuLg7Lly9HeHi4vD1hwAAZQFqyWoH+/YFNm+AQgh2bRDXE4COqRPnd3hs0aIATe/YgvEcPIDtb6/KkkBBg4UKM2LMHGzZsAMBd0Ymqw+AjqoZz9Gez2fB+bi7w+eeVzq9pIiQEye+8g8FTp2Lp0qUc5RFVg8FHVANCCBRu3IigJ5/U/hJnef7+wL33In/nTljqeJ8fkZEw+Ihq4to1oFUrec+eHoWEAMuWAcOHa10Jke6xq5OoJjZsAIqKtK6icrm5wOzZAN/HElWLIz6i6ggB3HwzcO6c1pVULSRE3uagZNk1IgPgiI+oOklJwOXLyo8zezawdStw+rS86b2oSHaHHj0KvP020L69suPn5wPz5imvk8jHccRHVJ0//Qn45BPllxGr+/nCQrne59atdT9HcDBw8SLQoEHdj0Hk4xh8RNVp1EiO0JT6+Wc5ejxzRh6vXj2gTx+ga9fr33PyJBAVVfdzhIUBmzYBDz2kvF4iH8XgI6rK5ctAixbuu2/PZAK+/x64/Xb5OD+/zlsPAZCb2s6cCbz0kjr1EfkgzvERVeXQIbkjgtpMJqBhQ7mdUatW158/dkzZcYuK5I4RRFQpf60LINK1/fvlrQJqad268u7QX38FJk9Wfo4jR5Qfg8iHccRHVJXU1Bt2QHCLkyeB//kfYN8+5cfKzAQcDuXHIfJRHPERVUXN0R4gm1qmTZPLjDVvDsTGArfeCrRrJ0eX8fHAxo3KzuHnJ/cHdMclWiIfwOYWoqrExgKffea+4/v5yUWve/eWj3NzgVtuAX75pe7HDAiQl03r11enRiIfw0udRFUJDnbv8UtKgISE649DQoBu3ZQds7jY/XUTeTEGH1FV1NrTrkePim8qN5mAvn3LPqf0IozZLEd9RFQhzvERVeXuu+V9dUq3IoqPBx57TN5qcPiw3OWhcWOgXz85v+eUnQ3s2aPsXDfdJAOViCrE4COqSufOshFFDUFBckWVylZV+e03eV/fb78pO8999yn7eSIfx+AjqkqHDnI1FaXefVcG2r33ypVgGjWSz1+9CqSkyF0VVqwALl1Sdh6rFejeXXm9RD6MXZ1E1bn1VuCHH7SuomZCQ4Fdu4AuXbSuhEi32NxCVJ3HH/eeLsmgIOCuu7SugkjXGHxE1Zk40Tt2Ng8OBp55Rt4bSESVYvARVSc8XN5g7g2dkuPGaV0Bke4x+Ihq4vnnlW0X5G5mM/Dww0DTplpXQqR7bG4hqgkhZLfk/v2eWbS6tiwWuYXSH/6gdSVEuscRH1FNmEzABx/I5hG9CQmRG88y9IhqhMFHVFOtWgHz5smg0QuTCWjTBpg+XetKiLwGL3US1YYQct3NAweAwkKtq5HzjgcPcrRHVAsc8RHVhskkd1No2VK9pczqymIBPvyQoUdUSww+otpq0ABITgYiI7XbBcFikcug9e+vzfmJvBgvdRLV1a+/AjExwJkzgN3umXOaTDL0Nm4E4uI8c04iH8MRH1FdNW4sb28YPVqGkbtZrXLd0L17GXpECnDER6SGb76RWwplZamzm0NpJpNcjmzqVGDWLG4yS6QQg49ILXl5MpiWLZOPc3KUHS8wUK7I0rUrsHgx0KmT4hKJiMFHpL78fGDzZmDOHODcOXkLRE1Hgf7+8pKmwwGMHQtMmiTv0yMi1TD4iNzp6FHZAZqUBOzbB5w/Ly9Vms3yEqbDARQVyZviO3WSzTJduwIPPOA9WyEReRkGH5EnFRQAly/LEWBJiQy3+vWBhg21rozIMBh8RERkKLydgYiIDIXBR0REhsLgIyIiQ2HwERGRoTD4iIjIUBh8RERkKAw+IiIyFAYfEREZCoOPiIgMhcFHRESGwuAjIiJDYfAREZGhMPiIiMhQGHxERGQoDD4iIjIUBh8RERkKg4+IiAyFwUdERIbC4CMiIkNh8BERkaEw+IiIyFAYfEREZCgMPiIiMhQGHxERGQqDj4iIDIXBR0REhsLgIyIiQ2HwERGRoTD4iIjIUBh8RERkKAw+IiIyFAYfEREZCoOPiIgM5f8BNPMiUWo/SsQAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 找到测量结果中出现几率最大的比特串\n",
    "cut_bitstring = max(prob_measure, key=prob_measure.get)\n",
    "print(\"找到的割的比特串形式：\", cut_bitstring)\n",
    "\n",
    "# 在图上画出上面得到的比特串对应的割\n",
    "node_cut = [\"blue\" if cut_bitstring[v] == \"1\" else \"red\" for v in V]\n",
    "\n",
    "edge_cut = []\n",
    "for u in range(n):\n",
    "    for v in range(u + 1, n):\n",
    "        if (u, v) in E or (v, u) in E:\n",
    "            if cut_bitstring[u] == cut_bitstring[v]:\n",
    "                edge_cut.append(\"solid\")\n",
    "            else:\n",
    "                edge_cut.append(\"dashed\")\n",
    "\n",
    "nx.draw(G, pos, node_color=node_cut, style=edge_cut, **options)\n",
    "ax = plt.gca()\n",
    "ax.margins(0.20)\n",
    "plt.axis(\"off\")\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "可以看到，在这个例子中 QAOA 找到了图上的一个最大割。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_______\n",
    "\n",
    "## 参考文献\n",
    "\n",
    "[1] Farhi, E., Goldstone, J. & Gutmann, S. A Quantum Approximate Optimization Algorithm. [arXiv:1411.4028 (2014).](https://arxiv.org/abs/1411.4028)"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "3b61f83e8397e1c9fcea57a3d9915794102e67724879b24295f8014f41a14d85"
  },
  "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.13"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}