docs/source/tutorials/qaoa.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "# Quantum Approximation Optimization Algorithm (QAOA)"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Overview"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "QAOA is a hybrid classical-quantum algorithm that combines quantum circuits, and classical optimization of those circuits. In this tutorial, we utilize QAOA to solve the maximum cut (MAX CUT) combinatorial optimization problem: Given a graph $G=(V, E)$ with nodes $V$ and edges $E$, find a subset $S \\in V$ such that the number of edges between $S$ and $S \\backslash V$ is maximized. And this problem can be reduced to that of finding the ground state of an antiferromagnetic Ising model whose Hamiltonian reads:\n",
    "\n",
    "$$H_C = \\frac{1}{2} \\sum_{i,j\\in E} C_{ij} \\sigma^{z}_{i} \\sigma^{z}_{j},$$\n",
    "\n",
    "where $\\sigma^{z}_{i}$ is the Pauli-z matrix on $i$th qubit and $C_{ij}$ is the weight of the edge between nodes $i$ and $j$. We set $C_{ij}=1$ for simplicity. If $\\sigma^{z}_{i}=\\sigma^{z}_{j}$, $i,j\\in S$ or $i,j\\in S \\backslash V$; if $\\sigma^{z}_{i}= -\\sigma^{z}_{j}$, $i\\in S, j\\in S \\backslash V$ or $i\\in S \\backslash V, j \\in S$. Obviously, the number of edges between $S$ and $S \\backslash V$ is maximized with the graph structure decoded from the ground state."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Setup"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "source": [
    "import tensorcircuit as tc\n",
    "import tensorflow as tf\n",
    "import networkx as nx\n",
    "\n",
    "K = tc.set_backend(\"tensorflow\")\n",
    "\n",
    "nlayers = 3  # the number of layers\n",
    "ncircuits = 2  # the number of circuits with different initial parameters"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Define the Graph"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "source": [
    "def dict2graph(d):\n",
    "    g = nx.to_networkx_graph(d)\n",
    "    for e in g.edges:\n",
    "        if not g[e[0]][e[1]].get(\"weight\"):\n",
    "            g[e[0]][e[1]][\"weight\"] = 1.0\n",
    "    nx.draw(g, with_labels=True)\n",
    "    return g\n",
    "\n",
    "\n",
    "# a graph instance\n",
    "# each node is connected to three nodes\n",
    "# for example, node 0 is connected to nodes 1,7,3\n",
    "example_graph_dict = {\n",
    "    0: {1: {\"weight\": 1.0}, 7: {\"weight\": 1.0}, 3: {\"weight\": 1.0}},\n",
    "    1: {0: {\"weight\": 1.0}, 2: {\"weight\": 1.0}, 3: {\"weight\": 1.0}},\n",
    "    2: {1: {\"weight\": 1.0}, 3: {\"weight\": 1.0}, 5: {\"weight\": 1.0}},\n",
    "    4: {7: {\"weight\": 1.0}, 6: {\"weight\": 1.0}, 5: {\"weight\": 1.0}},\n",
    "    7: {4: {\"weight\": 1.0}, 6: {\"weight\": 1.0}, 0: {\"weight\": 1.0}},\n",
    "    3: {1: {\"weight\": 1.0}, 2: {\"weight\": 1.0}, 0: {\"weight\": 1.0}},\n",
    "    6: {7: {\"weight\": 1.0}, 4: {\"weight\": 1.0}, 5: {\"weight\": 1.0}},\n",
    "    5: {6: {\"weight\": 1.0}, 4: {\"weight\": 1.0}, 2: {\"weight\": 1.0}},\n",
    "}\n",
    "\n",
    "example_graph = dict2graph(example_graph_dict)"
   ],
   "outputs": [
    {
     "output_type": "display_data",
     "data": {
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ],
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAb4AAAEuCAYAAADx63eqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABKxElEQVR4nO3dd1yV9f//8ceBAxwQAcU9crEUJ87EElducecgLU0t1DS1FE1zj9TcZmU2hH7ZF7WszDJnOdKcqCCIouJkiKABMs7vDz+QyIZzuM7hvO63W7dbnHOd9/WilCfveam0Wq0WIYQQwkSYKV2AEEIIUZIk+IQQQpgUCT4hhBAmRYJPCCGESZHgE0IIYVIk+IQQQpgUCT4hhBAmRYJPCCGESZHgE0IIYVIk+IQQQpgUCT4hhBAmRYJPCCGESZHgE0IIYVIk+IQQQpgUCT4hhBAmRYJPCCGESZHgE0IIYVIk+IQQQpgUCT4hhBAmRYJPCCGESZHgE0IIYVIk+IQQQpgUtdIFFFX0o2QCT0UScjee+KRU7DRq3KrYMah5DRxtrZQuTwghhIFSabVardJFFMa5m3FsOHiFQ6FRACSnpme+p1GboQW8XCvi296JJjUdlClSCCGEwTKq4PM/HsGi3SEkpaaRV9UqFWjU5szq4YZPm9olVp8QQgjDZzRDnU9DL5jElPR8r9VqITEljUW7gwEk/IQQQmQyih7fuZtxDPn8OIkpadnee3zpEHFH/h9p8VGYlymHY8/JaGo2zHzf2sKcbWPb0LiGQwlWLIQQwlAZRY9vw8ErJKVmD73Ea2d4cPArKnpPx7KaC2mPYrNdk5SaxsaDV9jk06IkShVCCGHgDD74oh8lcyg0Ksc5vYd/BWDvORSr6m4AqMtWyHaNVgsHLkcR8yhZVnsKIYQw/H18gacic3xdm55G8p0rpP/7kFubxhC5YSSxv39CekpytmtVQODpnNsRQghhWgy+xxdyNz7LloUMaY/jID2Vfy8fobLPMlRm5kRtX8jDo9so135ElmuTUtP5J+wOgxuVx87ODnNz8xKqXuiD7OEUQhSHwQdffFJqjq+rLJ7+gCvbvDdq2/JP/71l3xyDD2Dvob+oO603CQkJ2NjY4ODgkPmPvb19gb+2t7fH0tJSf9+wyFXeezjvsuqPUNnDKYTIl8EHn50m5xLNNbaYPzenp1Kpcm2nX89urPpmBunp6SQkJBAXF8fDhw+Ji4vL/Cfj6zt37hAcHJzr+xYWFnkGZX7hqdFo8qxVZJffHs6k/4Xg75fucTg0WvZwCiFyZfDB51bFDiv13RyHO20bdSbh1M9Y120O5mriT/6AjVPLbNdp1Ga4VS0LgJmZGfb29tjb2xepHq1Wy7///ptjcGa8Fhsby9WrV3MN1/T09EL1Mp//2tbW1qSCU/ZwCiF0yeD38UU/SsZz2f4cg0+blkrsH5/x+NIhVGoLyri9RLkOb6BSZx2KtFKbcXR6R4OZ/0lKSsoMxNx6lXm9lpSUhJ2dXZGD05jmOXPbw3k3YAbJty+jMnv6fZiXdaT62E+zXCN7OIUQOTH44AMYu/Uf9gbfy/OYstyoVNC1QeVStY8vJSWF+Pj4XIMxv68TEhIoU6ZMoYLz2ddKcp4zt//3dwNmUKZhB8o26ZrrZ0vj/3shRPEZ/FAnwHgvJ/4Mi87x5Jb8aNTm+Ho56aEq5VhYWODo6Iijo2ORPv/sPGdu4Zgxz5nb+1ZWVsUari3IPGdeezgLQvZwCiFyYhQ9PijcPE8GawszZvWoL/M8OqbVann8+HGhe5oZXz948AAg36C8bFaTownlSdVmD8i7ATNIib4BgEX56ji8/BqaWo2zXadRm/FuFxfGvVxPv/9RhBBGwyh6fPDfIgV5OoPyVCoVtra22NraUr169SK18ew8Z25BeSU5kVTLnHuF5Tq8gYVjTVTmFjwOPsz97Quo+sZaLMpVzXqf1HRC7iQUqUYhROlkNMEHT8OvcQ0HNh68woHL2fdykfoEC0tLOtWvjK+XkyxqMGAajQaNRkPlypVzvWbU1ye5HXI/x/esqrlm/rtto048vnSIxPB/sGjRO9u18UkpxS9YCFFqGFXwATSu4cAmnxbEPErmiwMXWffNDl7p5Y2dxoJrZ49QJ/0uq3wWKF2m0IHc9nDmSKUCch4GsNNY6KYgIUSpYPBndebG0daKHnUsKXtxB1+MbMmqV5vi17cFv/24XenShI7UdrDEnOxzuulJj0i8egpt6hO06Wk8uniA5JsXnu7nfM6zeziFEAKMsMf3rLi4OMqVK5f5dcuWLXnw4AFhYWE4OzsrWJkojhs3brBmzRq+/m47ZX3WZvv1TJueRtxhf1JiI0FlhoVjDSr2/wCL8tnnG5+kpNK5ngSfEOI/RtvjA3jw4AEODg6ZX5uZmdG7d2927dqlXFGiyE6dOsWwYcNo2rTp06+PHqJLw2o8v+vB3Maeqq+v4oUp/8cL726j6oiVWNdplq09FeD45A5tmrrz0UcfkZiYqP9vQghh8Iw6+J7v8QF4e3vz448/KlSRKKz09HR+/vlnOnToQN++ffHw8ODatWusXLmSWrVqMd7LCY26aKfMaCzM+WLKYA4fPszx48dxcXFh8+bNpKbmfPC5EMI0GHXwPd/jA+jUqRPnzp0jOjpamaJEgSQmJvLZZ5/h7u7OnDlzGDNmDFevXmXatGlZzlFtUtOBWT3csLYo3B/Vp3s43Whcw4H69euzY8cOAgMDCQgIoGHDhmzfvh0j2cIqhNAxow++53t8Go2GTp068csvvyhUlchLVFQU8+bNo3bt2uzatYuNGzdmDnFaWOS8+tKnTW1m9aiPtYV5tmHP56lUT8/ozOnggtatW7N//35Wr17NggULaNOmDQcOHNDRdyaEMBZGHXxxcXHZenzwdLhT5vkMy+XLlxk3bhwuLi5ERkZy8ODBzCHOgjxpwqdNbbaNbUPXBpWxUpuhUWf9o6tRm2GlNqNrg8psG9sm14MLVCoV3bp14/Tp00yaNIk333yTbt26cebMGV18m0III2A0R5bl5PXXX+fll19m1KhRWV6Pjo6mXr163Lt3D41Go1B1QqvVcvjwYVauXMnx48d5++238fX1zXPTekHEPEom8HQkIXcSiE9KwU5jgVvVsgz0KPwT2J88ecLnn3/OwoUL8fLyYuHChdSrJ8ebCVGalartDBkqVKhAkyZN2L9/Pz169FCgMtOWmppKYGAgK1eu5OHDh0yZMoXvvvsOGxsbnbTvaGuls7M3LS0tGT9+PCNHjmTVqlW0bt2aV199ldmzZ1OlShWd3EMIYViMeqgzp8UtGfr06SOrO0tYfHw8q1atol69emzcuJEPPviAkJAQ3nrrLZ2Fnr7Y2toye/ZsgoODsbS0xN3dndmzZ/Pw4UOlSxNC6JhRB19uPT74b54vPb3gT3MQRXPz5k3ee+896tSpw99//01gYCCHDx/G29sbMzPj+iNWsWJFVq1axenTp7l58yYuLi58/PHHJCUlKV2aEEJHjOun0nPy6vE5OztTrlw5/vnnn5ItyoScOXMGHx8fmjRpQmpqKqdOneK7776jZcuWSpdWbLVq1eKrr75i3759HDx4EFdXV7766ivS0gr/TEghhGEx6uDLq8cHMtypD+np6ezevZuOHTvSu3dvmjRpwtWrV1m1ahW1a9dWujyda9iwIbt27SIgIIDNmzfTpEkTdu3aJXsAhTBiRruqMzU1FY1Gw5MnT3IdTjt27Bhjx44lKCiohKsrfZKSkvD39+fjjz/GysqKqVOnMnjwYCwtLZUurcRotVp++eUX/Pz8sLOzY+nSpbz00ktKlyWEKCSjDb7o6GhcXFyIjY3N9Zq0tDSqVavGsWPHqFu3bglWV3pER0fzySefsGHDBjw8PJg6dSodO3Ys0N670iotLY2AgADmzJlDw4YNWbJkCY0aNVK6LCFEARntUGd+w5wA5ubmcmh1EYWGhvL222/j7OxMREQE+/btY/fu3XTq1MmkQw+e/rkaMWIEly9fpnPnznTu3JkRI0YQERGhdGlCiAIw2uDLa2HLs2Ser+C0Wi1//vknffv2xdPTkwoVKhAcHMwXX3yBu7u70uUZHCsrKyZPnkxYWBh16tShefPmTJ48maioKKVLE0LkwWiDryA9PoDOnTtz6tSpPIdETV1qairff/89rVu3ZtSoUbzyyitERESwYMEC2cRdAHZ2dsybN49Lly6RlpaGm5sb8+bNIyEhQenShBA5MNrgK2iPz8bGho4dO7J79279F2VkEhISWL16Nc7Ozqxdu5aZM2cSEhKCr68vZcqUUbo8o1O5cmXWrVvHyZMnCQ0NxdnZmXXr1vHkyROlSxNCPMNog6+gPT6Q4c7nRUZGMn36dOrUqcPRo0f57rvv+Ouvv+jbty/m5kV79p34T926dQkICGDPnj3s3r0bNzc3AgIC5DAFIQyE0QZfQXt8AL169WLv3r0kJyfrtygDd/bsWV577TUaN25MUlISJ0+ezBziFLrXtGlTfv31V7Zs2cK6deto1qwZu3fvlj2AQijMqIOvoD2+SpUq4e7ubpLPXtNqtezZs4fOnTvTs2dPGjZsSHh4OGvWrKFOnTpKl2cSvLy8OHbsGHPnzmXq1Kl4eXlx/PhxpcsSwmQZbfDl9iy+3Hh7e5vUcGdycjJbtmyhUaNGTJ8+nREjRnDt2jWmT59e4F8YhO6oVCr69etHUFAQI0aMYNCgQfTr14/g4GClSxPC5Bht8BWmxwf/HVpd2oeZYmJiWLRoEbVr1+b7779n9erVnD17lhEjRpjUKSuGSq1WM3r0aEJDQ/H09KR9+/aMGjWKmzdvKl2aECbDaIOvMItbAFxdXbG1teXUqVN6rEo5V65cYfz48Tg7OxMeHs7evXszhzhNfcO5IbK2tmbatGmEhoZSpUoVmjZtyrRp04iJiVG6NCFKPaMNvsIsbsmQ0esrLbRaLUeOHKF///68+OKLODg4cPHiRbZs2ULDhg2VLk8UgIODA4sXLyYoKIhHjx7h6urK4sWLefz4sdKlCVFqGW3wFbbHB6VnW0Nqair/93//x4svvsiIESPo1KkTERERLFq0iKpVqypdniiCatWqsWnTJo4ePcq5c+dwdnbmk08+ISUlRenShCh1jPaQ6ooVK3LhwgUqV65c4M+kpaVRtWpVTpw4YZSP0Hn06BFbtmxh9erVVK1alalTp+Lt7S1770qhf/75Bz8/PyIiIli4cCGDBg0yuof6CmGojPJvklarLfSqTnh6uHDPnj2Nbrjz9u3b+Pn5Ubt2bQ4fPkxAQEDmEKeEXunUokUL9u7dyyeffMLy5ctp2bIle/fuVbosIUoFowy+x48fY2FhgZWVVaE/a0zbGs6fP8/IkSNxd3fn8ePHnDhxgsDAQF588UWlSxMlpHPnzpw4cYIZM2Ywfvx4OnfuzMmTJ5UuSwijZpTBV9itDM/q0qULJ0+e5MGDBzquSje0Wi2//fYbr7zyCt26dcPNzY3w8HDWrl0rzxQ0UWZmZgwaNIiLFy8yaNAg+vbty6BBgwgNDVW6NCGMklEGX1EWtmQoU6YM7du359dff9VxVcWTnJzMV199RZMmTZg2bRrDhg3j2rVr+Pn5Ub58eaXLEwbAwsKCcePGERYWhoeHB23btmXcuHHcvn1b6dKEMCpGGXxF2crwLEPa1hAbG8vixYupU6cO3377LStWrOD8+fO8/vrrRRrKFaWfjY0Nfn5+hIaGYm9vT6NGjfDz8yMuLk7p0oQwCkYZfMXp8cHTQ6t/++03RR8XEx4ezsSJE3FyciI0NJQ9e/bw+++/88orr8iGc1Eg5cuX56OPPuLcuXNERUXh4uLC8uXLSUxMVLo0IQyaUQZfcXt8VapUwc3NjUOHDumuqAI6duwYAwYMoHXr1tja2nLhwgW++uorGjduXOK1iNKhRo0abN68mUOHDnHs2DFcXFzYvHkzqampSpcmhEEyyuArbo8PSnYze1paGtu3b6dt27YMHz4cLy8vIiIiWLJkCdWqVSuRGkTpV79+fXbs2EFgYCABAQE0atSIHTt2lPrzaYUoLKMMvuL2+KBkDq1+/Pgx69evx8XFhRUrVjBlyhTCwsKYOHEitra2eruvMG2tW7dm//79rFq1ivnz59OmTRuTfCSXELkx2uArbo+vfv36WFlZcfbsWd0U9Yw7d+4wc+ZMateuzf79+9m6dSvHjh1j4MCBsuFclAiVSkW3bt04ffo0kyZNYvTo0XTr1o0zZ84oXZoQijPK4NPFUKdKpdL5cOeFCxd44403cHd3Jz4+nmPHjrFjxw7atm2rs3sIURhmZmYMGzaMkJAQevXqRY8ePRg2bBjh4eFKlyaEYowy+HQx1Am62dag1WrZu3cv3bp1o0uXLjg5OREWFsb69etxcnIqdo1C6IKlpSUTJkwgLCwMNzc3WrVqxYQJE7h7967SpQlR4owy+HTR4wNo27YtN27c4MaNG4X+7JMnT/jmm29o2rQp7777Lq+++ioRERHMmjULR0fHYtcmhD7Y2toyZ84cQkJCsLCwwN3dndmzZxMfH690aUKUGKMMPl31+NRqdaEPrX7w4AFLly6lTp06bN26lWXLlhEUFMQbb7whG86F0ahYsSKrVq3i1KlT3LhxA2dnZ1atWkVSUpLSpQmhd0YZfLrq8UHBhzuvXbvGpEmTqFevHpcuXWL37t2ZQ5yy4VwYq9q1a/P111/zxx9/cODAAVxdXfnqq69IS0tTujQh9MYog09XPT6AV155hePHj/Pw4cMc3//7778ZNGgQLVu2xNramqCgIL755huaNGmik/sLYQgaNWrErl27CAgIYPPmzTRp0kTv232EUIrRPYg2JSUFa2trUlJSdNbT6tmzJyNGjODVV18Fnm4437VrFytXruTWrVtMnjyZUaNGUbZsWZ3cTwhDptVq+fnnn/Hz88PBwYGlS5fSrl07pcsSQmeMrscXFxeHvb29TocXM7Y1/Pvvv2zcuBE3NzeWLl3KpEmTCAsLY9KkSRJ6wmSoVCp69+7NuXPnGDNmDMOHD6d3794EBQUpXZoQOmF0waeLzevPa926NTt37uSFF15g7969fPnllxw/fpxBgwahVqt1ei8hjIW5uTkjR44kNDSUTp060blzZ0aMGEFERITSpQlRLEYXfLpc2HLx4kVGjx5Nhw4dsLOzY+XKlezcuZN27drJghUh/sfKyorJkycTFhZG7dq1ad68OZMnTyYqKkrp0oQoEqMJvuhHyWw6FM6KI9E8bjGCydvOsOlQODGPkgvVjlarZd++ffTo0YNOnTpRu3ZtwsLCeOeddzh16pSeqhfC+NnZ2TF//nwuXbpEWloabm5uzJs3j4SEBKVLE6JQDH5xy7mbcWw4eIVDoU9/u0xOTc98T6M2Qwt4uVbEt70TTWo65NrOkydP2LZtGx9//DHJyclMnTqV4cOHo9FogKfHjfXq1Ytr165Jb0+IArh69SqzZ89m3759zJo1i3HjxmFpaal0WULky6CDz/94BIt2h5CUmkZeVapUoFGbM6uHGz5tamd5Ly4ujs8++4y1a9fi6urK1KlT6datG2ZmWTu7Wq2WevXqsXPnTtmqIEQhnD17Fj8/Py5fvsyCBQsYOnRotr9fQhgSg/3T+TT0gklMyTv0ALRaSExJY9HuYPyPRwAQERHB5MmTqVu3LkFBQfz000+ZQ5w5/aVUqVQ6ObtTCFPTtGlTfv31V7Zs2cK6devw8PDg119/lT2AwmAZZI/v3M04hnx+nMSUnE+PSIm9xe0vJlDGzZMKvadlec/SHJyv/8Lx3dsYPXo077zzDjVq1CjQfQ8cOMB7773HP//8U+zvQQhTpNVq+eGHH5g5cyaVK1dm6dKltGnTRumyhMjCIHt8Gw5eISk19yOTYn/fhFVV5xzfe5KSTlK99kRERPDRRx8VOPQA2rVrx7Vr14iMjCx0zUKIpyMn/fr1IygoiNdee41BgwbRr18/goODlS5NiEwGF3zRj5I5FBqV6/Dm40uHMNOUQVMrl3k4MzMi0+x4oir8JLuFhQXdu3fnp59+KvRnhRD/UavVjB49mtDQUNq2bcvLL7/M6NGjuXnzptKlCWF4wRd4KvfeVnryv8T9GUC5jm/m2YYKCDxdtF6bt7e3Th9OK4Qps7a25r333iMsLIxKlSrRtGlT3nvvPWJiYpQuTZgwgwu+kLvxWbYsPCvu8FZsm7yC2q5Cnm0kpaYTcqdoe4u6devGkSNH5PlkQuiQg4MDS5YsISgoiISEBFxdXVm8eDGPHz9WujRhggwu+OKTUnN8/cm9qyRdP4ddS+8CtpNSpPuXLVsWT09PfvvttyJ9XgiRu2rVqrFp0yaOHj3K2bNncXFxYdOmTaSkFO3vqxBFYXDBZ6fJ+WzMpBtBpD68R+TGN7i5zof4Ezv59/JR7nw5KZd2LIpcg2xrEEK/XFxc+P777/nxxx8JDAzE3d2d77//nvT0nEd7hNAlg9vOsOlQOKv+CM023JmekoQ2OTHz6/gTO0h9eI/yXcdjbmOf5VqN2ox3u7gw7uV6RaohMjKSJk2acPfuXSwsih6gQoiC+eOPP5gxYwYAS5YsoUuXLgpXJEozg+vxDWye8/YDMwsN5rblMv9RWWhQqS2zhR6AFhjoUfBtDM+rUaMGderU4a+//ipyG0KIguvcuTMnTpzg/fffx9fXl86dO8t+WqE3Bhd8FWytaO9SkfyOy3R4aXi2zevw9PiyDq4VcbS1KlYdMtwpRMkyMzNj8ODBXLp0iYEDB9KnTx8GDx5MaGio0qWJUsbggg9gvJcTGrV5kT6rUZvj6+VU7BoyHk5rYCPBQpR6FhYWvPXWW4SFhdGsWTPatm3LuHHjuH37ttKliVLCIIOvSU0HZvVww9qicOVZW5gxq4cbjWs4FLuGxo0bk56ezsWLF4vdlhCi8MqUKZN5+LWdnR2NGjXCz8+PuLg4pUsTRs4ggw/Ap01tZvWoj7WFeb7DnqSno7EwY1aP+tmezlBUKpUqs9cnhFCOo6Mjy5cv5+zZs9y/fx9nZ2eWL19OYmJi/h8WIgcGG3zwNPy2jW1D1waVsVKboVFnLVejNsNKbYb9o+u8YnZRZ6GXQeb5hDAcNWvW5IsvvuDw4cMcO3YMFxcXNm/eTGpqznt/hciNwW1nyE3Mo2QCT0cScieB+KQU7DQWuFUty0CPGjyKuYuHhwcnT56kbt26OrtnSkoKlStX5sKFC1SrVk1n7Qohiu/48ePMmDGDe/fusWjRIvr16ycPkRYFYjTBl58lS5Zw9OhRnR8wPWzYMNq3b8+4ceN02q4Qovi0Wi179uzBz88PKysrli5dSocOHZQuSxg4gx7qLIwpU6YQGhqq86FJGe4UwnCpVCq6d+/O6dOnmTRpEqNHj6Zbt26cOXNG6dKEASs1PT54evrDmDFjuHjxIjY2Njpp8+HDh9SsWZPbt29ja2urkzaFEPrx5MkTPvvsMxYtWkSHDh1YsGAB9eoV7QQnUXqVmh4fPD39oVWrVixZskRnbdrb29OmTRt+//13nbUphNAPS0tLJkyYQFhYGG5ubrRq1YoJEyZw7949pUsTBqRUBR/AypUr+eSTTwgLC9NZm7KtQQjjYmtry5w5cwgJCcHCwoIGDRowZ84cedyYAErZUGeGFStWsHfvXvbs2aOTVV43btzAw8ODu3fvolbn/PQIIYThioiI4MMPP2TPnj3MmDEDX19frKwKdqxh9KNkAk9FEnI3nvikVOw0atyq2DGoeY1iH40olFEqgy8lJYWmTZsyf/58BgwYoJM2mzVrxpo1a3j55Zd10p4QouQFBQUxc+ZMgoKCmDdvHj4+Ppib53w84rmbcWw4eIVDoVEAWZ4Yo1GboQW8XCvi296JJjUdSqB6oSulMvgADh48yIgRI7h06ZJOFqXMnTuXhIQEVq5cqYPqhBBK+uuvv5g+fTrx8fEsXryYXr16ZRkd8j8ewaLdISSlppHXT0iV6un5wLN6uOn8AA2hP6U2+AB8fHyoXr06y5YtK3ZbZ86cYdCgQYSFhckmWSFKAa1Wy88//4yfnx8ODg4sXbqUdu3a/S/0gklMKfhDca11fGSi0K9SHXx3796lYcOG/Pnnn9SvX79YbWm1WmrVqsWePXto0KCBjioUQigtLS0Nf39/5syZg1PrTkS6DCQ57b8fi9rUFGJ+30hSxFnSkx6hdqhCufYjsa7XIks71hbmbBvbRieH5Av9KnWrOp9VpUoVZs+ezfjx44v9eKGMQ6tlM7sQpYu5uTkjR47k8uXLUP8VklLSsryvTU9DXbYCVYYtpea723B4+TWiflxGalzWLRJJqWlsPHilJEsXRVSqgw9g/PjxxMbG8t133xW7LdnWIETp9ShVRWS6PSqzrD8WzSw1OLw0HLVDZVQqM2ycWqG2r0zy3awhp9XCgctRxDxKLsmyRRGU+uBTq9Vs3LiRadOmFXsPj5eXFyEhIdy9e1dH1QkhDEXgqcgCXZf2+AEpsbewrPhCtvdUQODpgrUjlFPqgw+gbdu2dO3alblz5xarHUtLS7p27crPP/+sm8KEEAYj5G58li0LOdGmpRK9awW2jTph4Vgz2/tJqemE3EnQV4lCR0wi+ACWLl2Kv78/QUFBxWpHhjuFKJ3ik/J+rp9Wm070zyvBXE35Lm/l0U6KrksTOmYywVepUiXmzZuHr69vsRa6dO/enUOHDvH48WMdVieEUJqdJvdTmbRaLTG715L2OI6K/WaiMs/9WjuNhT7KEzpkMsEHMHbsWBITE9m6dWuR2yhXrhwtW7Zk7969OqxMCKE0typ2WKlz/pEY+9sGUmJuUmngHMwscj+mTKM2w61qWX2VKHSkVO/jy8mJEyfw9vYmODgYBweHIrWxdu1azp49y5YtW3RbnBBCMdGPkvFctj/bPF/qw/vc+mQUmFugMvvveLPy3cZj6571obdWajOOTu8oZ3gaOJMLPoBx48ZhaWnJunXrivT5iIgIWrVqxZ07d3I9508IYXzGbv2HvcH38jymLDcqFXRtUJlNPi3yv1goyqSGOjMsXryY77//ntOnTxfp87Vr16Zq1aocO3ZMx5UJIZQ03ssJjbpov8xq1Ob4ejnpuCKhDyYZfI6OjixatAhfX1/S0wt+Ht+zvL295RQXIUqZJjUdmNXDDWuLwv1otFKrmNXDTY4rMxImGXwAo0aNQqVS8eWXXxbp87KtQYjSyadNbWb1qI+1hTn5nUevUoFalY729Ha83R1LpkBRbCY5x5fh9OnTdO/enUuXLuHoWLg/tFqtlho1arB//35cXV31VKEQQinnI+PYePAKBy5HoeLp5vQM5qShVlvQwbUib7evx4YF07l//z47duzAzMxk+xNGw6SDD2DChAmkpKTw6aefFvqzb7/9NnXq1OH999/XQ2VCCEMQ8yiZwNORhNxJID4phccPogj5ez/7Pl+UuXrzyZMndOrUiQ4dOjB//nyFKxb5Mfngi4uLo379+vz444+0atWqUJ/ds2cPCxYs4MiRI3qqTghhaOLj46lWrRqxsbFYWlpmvn7//n1atmzJihUrGDRokIIVivyYfJ/cwcGBZcuW4evrS1paWv4feEaHDh24ePEi9+/f11N1QghDY2dnh7Ozc7ZV4ZUqVeKHH37A19eXs2fPKlOcKBCTDz6A1157DWtraz777LNCfc7KyoouXbrIodVCmJh27drx119/ZXu9WbNmbNiwgb59+8ovxAZMgo+nD5nduHEjH374IVFRUYX6rGxrEML0tGvXLtcpjsGDB+Pj48OAAQN48uRJCVcmCsLk5/ieNWXKFOLi4gp1FFlsbCy1a9fm7t272NjY6LE6IYShiIyMpFmzZty/fx9VDnse0tPT6d+/P5UrV2bTpk05XiOUIz2+Z8ydO5fffvutUItVypcvT/Pmzdm3b58eKxNCGJIaNWpQpkwZQkNDc3zfzMyMrVu3cuTIET755JMSrk7kR4LvGXZ2dqxYsYLx48eTmpr3s7meJZvZhTA9eQ13ApQtW5Yff/yR+fPnc+DAgRKsTORHgu85Q4YMoXz58mzcuLHAn/H29uann34q9KpQIYTx8vT0zHGBy7Pq1avHt99+y9ChQ7l27VoJVSbyI8H3HJVKxYYNG5g/fz537twp0Gfq1q1LpUqVOHHihJ6rE0IYitxWdj6vY8eOzJo1iz59+pCQkFAClYn8SPDloH79+owePbpQJ7LIcKcQpsXd3Z2oqKgCbVuYMGECbdq0YcSIEUU+GF/ojgRfLmbPns2hQ4c4dOhQga6XbQ1CmBYzMzNefPHFAi2GyxhJioqKYt68eSVQnciLBF8ubG1tWbVqFePHjyclJSXf61u0aEFcXBxhYWElUJ0QwhAUdLgTwNLSku3bt/PVV18RGBio58pEXiT48tC/f3+qV6/O2rVr873WzMyM3r17S69PCBOS38rO51WuXJmdO3fy9ttvy7FmCpIN7PkICwvjxRdf5Ny5c1SvXj3Pa3/55ReWLVvG4cOHS6g6IYSSEhMTqVChAlFRUYU6wGLbtm3MmDGDEydOULFiRT1WKHIiPb58ODs78/bbbzNlypR8r+3UqRPnzp0jOjq6BCoTQijN2tqaxo0bF3pF96uvvsqwYcMYOHCgHGumAAm+AvDz8+PEiRP88ccfeV6n0Wjo3Lkzv/zySwlVJoRQWmHm+Z61YMEC7O3teeedd/RQlciLBF8B2NjYsGbNGiZMmEBycnKe18q2BiFMi6enZ5GeyWlmZoa/vz9//vmnHGtWwmSOrxB69+5N27Zt8fPzy/Wa6Oho6tWrx71799BoNCVYnRBCCVFRUTg7OxMTE4O5uXmhPx8eHk7btm3Ztm0bXl5eui9QZCM9vkJYs2YNK1as4MaNG7leU6FCBZo0aSKHVgthIipWrEiVKlW4cOFCkT6fcazZkCFD5FizEiLBVwh169Zl0qRJTJ48Oc/rvL29ZbhTCBNS1OHODJ06dWLmzJl4e3vz6NEjHVYmciLBV0jvv/8+58+f59dff831mj59+vDTTz/J0URCmIiiLnB51sSJE2nZsqUca1YCJPgKSaPRsG7dOiZOnEhSUlKO1zg7O1OuXDlOnjxZwtUJIZSgi+BTqVRs3LiRe/fuMX/+fB1VJnIiwVcE3bt3p3Hjxnz00Ue5XiNndwphOpycnEhKSspz/r8grKys2LFjB1u2bGH79u06qk48T4KviFavXs3atWu5evVqju/LtgYhTIdKpSr08WW5yTjW7K233uLcuXM6qE48T4KviF544QWmTZvGpEmTcny/devWREdHEx4eXsKVCSGUoIvhzgzNmzdn3bp19O3bl6ioKJ20Kf4jwVcMU6ZMISwsLMchTTMzM3r16iXDnUKYCF31+DIMGTKEoUOHyrFmeiDBVwyWlpZs2LCBSZMm8e+//2Z7X7Y1CGE6mjVrRnh4OA8fPtRZmwsXLsTe3j7XkSVRNBJ8xdSpUydat27NkiVLsr3XuXNnTp8+TUxMjAKVCSFKkoWFBS1atODYsWM6azPjWLPDhw+zadMmnbVr6iT4dGDlypV88skn2R5Ca21tTceOHdm9e7dClQkhSpKuhzsB7Ozs+PHHH/nwww85dOiQTts2VRJ8OlC9enX8/PyYMGECzx99KtsahDAdnp6eOlvg8iwnJycCAgIYMmQIEREROm/f1Mgh1TqSkpJCs2bNmDdvHgMGDMh8/f79+7i4uHDv3j2srKwUrFAIoW8PHz6kRo0axMbGYmFhofP216xZw5YtWzhy5Ai2trY6b99USI9PRywsLNiwYQPvvvtulrP2KlWqRMOGDTlw4ICC1QkhSoK9vT316tXjzJkzemn/nXfeoXnz5rz++utyrFkxSPDpUPv27Wnfvj0LFizI8rpsZhfCdOhruBOebpT/5JNPuH37drafM6LgJPh0bPny5WzZsoXg4ODM1zLm+eQ3NCFKP11uZM9JxrFmX3zxBTt27NDbfUozCT4dq1KlCrNnz2b8+PGZC11cXV0pW7Ysp0+fVrg6IYS+Zazs1OfyiSpVqrBz507GjRvH+fPn9Xaf0kqCTw98fX158OAB3333XeZrMtwphGmoWbMmVlZWXLlyRa/3ad68OWvXrsXb25vo6Gi93qu0keDTA7VazYYNG5g2bRrx8fGAbGsQwpToe7gzw9ChQxkyZAiDBg0iJSVF7/crLST49KRt27Z07dqVuXPnAtCmTRvu3LnDtWvXlC1MCKF3JRV88PRYszJlyjB58uQSuV9pIMGnR8uWLcPf35+goCDMzc3l0GohTISnp6fOT3DJjbm5Od9++y0HDhzg008/LZF7GjsJPj2qWLEi8+fPx9fXF61WS58+fST4hDABDRs25O7duyX2SKGMY83mzJnD4cOHS+SexkyCT8/GjBlDUlISW7dupUuXLpw8eZIHDx4oXZYQQo/Mzc158cUXS6zXB+Ds7MzWrVt59dVX5VizfMiRZSXg5MmT9OnTh+DgYIa9MZYXvAajqVKP+KRU7DRq3KrYMah5DRxt5UgzIUqLhQsX8vDhQ5YvX16i9129ejVffvklR48epUyZMiV6b2MhwVdCho6fzjUbV6IsKpGelobWTJ35nkZthhbwcq2Ib3snmtR0UKxOIYRuHDx4ED8/P50+pqggtFoto0aNIiEhge+//x4zMxnYe54EXwnwPx7Bwl+CSXySiiqPP4QqFWjU5szq4YZPm9olV6AQQuf+/fdfKlasSHR0NNbW1iV67+TkZLy8vOjevTtz5swp0XsbA/lVQM/8j0ewaHcwSanpeYYegFYLiSlpLNodjP/xiJIpUAihFzY2NjRs2JCTJ0+W+L0zjjX7/PPP2blzZ4nf39Cp879EFNW5m3Es2h1CYkrWMzpvrByY5Wtt6hPKNutB+VfeAiAxJZ1Fu0NoXMOBxjUcSqpcIYSOZezne/nll0v83lWrVmXnzp10794dJycnGjVqVOI1GCrp8enRhoNXSEpNy/b6C1MDM/+pMXErKrUlNm7tslyTlJrGxoP6PfJICKFfJbmRPSctWrRgzZo1cqzZcyT49CT6UTKHQqPIbwb138tHMbexx6qme5bXtVo4cDmKmEfJeqxSCKFPbdu25dixY4o+mWXYsGEMHjxYjjV7hgSfngSeiizQdY+C9lGmYUdUKlW291RA4OmCtSOEMDyVK1emYsWKXLx4UdE6Fi1aRJkyZXj33XcVrcNQSPDpScjdeJJT8/4tL/XhfZJvXqBMo045vp+Umk7InQR9lCeEKCFKD3fC0w31AQEB7Nu3j88++0zRWgyBBJ+exCel5nvNowv7sarRAAuHKnm0I0MTQhizjOfzKc3e3p5du3Yxe/Zs/vzzT6XLUZQEn57YafJfMPv4wn5sG3bMpx0LXZUkhFCAp6en4j2+DBnHmg0ePJjr168rXY5iJPj0xK2KHVbq3P/zJkUGk/YoJttqzmdp1Ga4VS2rj/KEECXExcWFx48fExlpGPP1r7zyCu+//z7e3t48fvxY6XIUIcGnJwOb18jz/ccX9mHj0hYzK5tcr0l68oTK/0Ygh+sIYbxUKpXBDHdmmDx5Mk2bNuWNN94wyZ8vEnx6UsHWivYuFclhsSYAjt0mUKH31Fw/r1JBA/s03n/nbby8vDhw4ICeKhVC6JshDXfC0zDetGkTN27cYNGiRUqXU+Ik+PRovJcTGrV5kT6rUZuzdERHLl26xOjRoxk7dixeXl4cPHhQt0UKIfTOEFZ2Pk+j0bBz504+/fRTfvjhB6XLKVFySLWeZZzV+fyxZXmxtjBjVo/6WQ6qTk1NJSAggAULFlCzZk3mzZunyDFIQojCe/LkCY6Ojty6dQs7Ozuly8ni5MmT9OjRgwMHDtCwYUOlyykR0uPTM582tZnVoz7WFua5DntmUKnA2sI8W+gBqNVqRo4cSUhICCNHjuSNN96gY8eOJr8sWQhjYGlpiYeHB8ePH1e6lGxatmzJ6tWr8fb2JiYmRulySoQEXwnwaVObbWPb0LVBZazUZmieW+2pUZthpTaja4PKbBvbJs9HEqnVal5//XVCQkJ47bXXGDlyJJ06dTK4YRQhRFaGONyZYfjw4QwcONBkjjWToc4SFvMomcDTkYTcSSA+KQU7jQVuVcsy0KNoT2BPSUlh69atLFy4kLp16zJv3jw8PT31ULkQojh+/fVXli9fzv79+5UuJUdpaWn06dOHunXrsm7dOqXL0SsJvlIiJSWFb775hoULF+Lk5MS8efNo27at0mUJIf4nLi6OmjVrEhsbi4WFYR5M8fDhQ1q3bs3UqVMZM2aM0uXojQx1lhIWFhaMHj2ay5cvM3jwYIYNG0bXrl05duyY0qUJIQAHBwfq1KnD2bNnlS4lVxnHmn3wwQcGOyyrCxJ8pYylpSVjxowhNDSUAQMGMGTIELp162aQk+pCmBpDnufL4OLiwjfffMPgwYO5ceOG0uXohQRfKWVpacnYsWMJCwujX79+vPrqq3Tv3p2///5b6dKEMFmenp4GdYJLbrp27crUqVNL7bFmMsdnIpKTk/nyyy9ZvHgxDRs2ZO7cubRq1UrpsoQwKdevX6d169bcuXMnx2dwGhKtVsvrr79OYmIi27ZtM/h6C0N6fCbCysqKt956i7CwMHr37s2AAQPo2bMnJ0+eVLo0IUzGCy+8gIWFBeHh4UqXki+VSsWnn37KjRs3WLx4sdLl6JT0+ExUcnIyX3zxBUuWLKFJkybMnTuXFi1aKF2WEKXekCFD6N69OyNHjlS6lAK5ffs2rVq1YsOGDXh7e2e+Hv0omcBTkYTcjSc+KRU7jRq3KnYMal60rVklSYLPxCUlJWUGYLNmzZg7dy7NmzdXuiwhSq3169dz7tw5Pv/8c6VLKbATJ07Qq1cv9u/fT5p9DTYcvMKh0CgAklP/O45RozZDC3i5VsS3vRNNajooU3A+JPgE8DQAP//8c5YuXUrz5s2ZO3cuHh4eSpclRKlz9uxZhg4dSnBwsNKlFIq/vz+zv/4dqzZDSU5LJ6/kUKmeHrQ/q4dbnidRKUXm+ATw9KT2iRMnEh4eTufOnenduzd9+/blzJkzSpcmRKnSqFEjbt++TXR0tNKlFI5TO8yaDyQpNe/QA9BqITEljUW7g/E/HlEi5RWGBJ/IQqPR8M4773DlyhU6dOhAz5496devn0FvuhXCmJibm9OmTRuOHj2qdCkFdu5mHIt2h5Cm+u8xa/GnfuLOV5O5vrwv0T+vyvFziSnpLNodwvnIuBKqtGAk+ESOrK2tmTRpEuHh4bRv354ePXrQv39/zp07p3RpQhg9Y9jI/qwNB6+QlJqW5TW1rSP2bV/FtnGXPD+blJrGxoNX9FleoUnwiTxZW1szefJkrly5wksvvUS3bt0YMGAA58+fV7o0IYxWu3btjGIjOzxdvXkoNCrb8KaNa1tsXF7EzDrv5wtqtXDgchQxj5L1WGXhSPCJArGxseHdd98lPDwcT09PunbtysCBAwkKClK6NCGMTqtWrTh79iyJiYlKl5KvwFORxW5DBQSeLn47uiLBJwrFxsaGKVOmcOXKFV588UW6dOnCoEGDuHDhgtKlCWE0ypQpg7u7O//884/SpeQr5G58li0LRZGUmk7InQQdVVR8EnyiSMqUKcPUqVMJDw+ndevWdO7cmcGDB3Px4kWlSxPCKBjLcGd8UqqO2jGcB9xK8IliKVOmDNOmTSM8PJyWLVvSsWNHXn31VQlAIfLh6elpFAtc7DRqHbVjOM8glOATOlGmTBnee+89wsPD8fDwoEOHDgwZMoRLly4pXZoQBsnT05OjR4+Snl68YUR9un79OneC/4HU7L01bXoa2tQnkJ4G2nS0qU/Qpqfl0MrTE13cqpbVd7kFJsEndMrW1pbp06dz9epVmjZtipeXF8OGDTO6UyqE0LcqVarg6OhocH83YmNj+fTTT3n55Zdp3rw5Ze5fwMIye2/t4ZHvuLGiP/HHA3l88QA3VvTn4ZHvcmxTCwz0qKHnygtOjiwTepWQkMD69etZtWoVXbp0Yfbs2bi5uSldlhAG4fXXX+fFF19k3LhxitaRlJTEzz//jL+/PwcOHKBr1674+PjQrVu3p8/23PoPe4Pv5XtiS05UKujaoDKbfAznEHzp8Qm9Klu2LH5+foSHh+Pu7s7LL7+Mj48Ply9fVro0IRSn5Eb2tLQ09u/fz+jRo6lWrRqbNm3C29ubGzdu8P3339OnTx8sLS0BGO/lhEZtnk+LOdOozfH1ctJl6cUmwSdKRNmyZZk5cyZXrlyhfv36tGvXjtdee43Q0FClSxNCMSUdfFqtlrNnz/Lee+9Rq1Ytpk2bRoMGDQgKCuKPP/7gjTfewN7ePtvnmtR0YFYPN6wtChcZ1hZmzOrhRuMaDjr6DnRDgk+UKDs7O2bNmsWVK1dwdXXF09OTESNGEBYWpnRpQpQ4V1dXEhISuHXrll7vc/36dZYsWUKjRo3o27cvlpaW/P7775w+fZqpU6dSvXr1fNvwaVObsS0rok1JJr9nsatUYG1hzqwe9eXpDEJksLe354MPPuDKlSs4OzvTtm1bRo4cyZUrhnWmnxD6pFKp8PT01Mt+vtjYWD777LPMRSrXr19n06ZNXL16lUWLFtGgQYNCt3l4yyL621+nq3tlrNRmaNRZI0SjNsNKbUbXBpXZNraNQYYeyOIWYSDi4uJYu3Yta9eupVevXnzwwQc4ORnWvIAQ+rB8+XJu3rzJ2rVri91WxiKVgIAA9u/fT9euXRk+fDjdu3fPnK8rqoz5wODgYDQaDTGPkgk8HUnInQTik1Kw01jgVrUsAz3kCexCFEpcXByrV69m/fr19O7dmw8++IB69eopXZYQenP06FEmTpzIqVOnivT59PR0Dh06hL+/Pzt37qRZs2b4+PjQv3//HOfriiI1NRUPDw8+/PBDBgwYoJM2lSTBJwzSgwcPWL16NRs2bMDb25tZs2ZRt25dpcsSQueSk5NxdHTkzp07lC1bsE3eWq2W8+fPExAQwLfffkvFihUZPnw4Q4cOLdB8XWFt2rSJ7777jgMHDqBS5TfDZ/hkjk8YpHLlyjFv3jzCwsKoXr06LVu25M033+TatWtKlyaETllZWeHh4cHx48fzvfbGjRssXbqUxo0b4+3tjVqt5rfffuPMmTNMmzZNL6H34MEDPvzwQ1avXl0qQg8k+ISBK1euHPPnzycsLIyqVavSokULxowZQ0REhNKlCaEzeS1wefDgAZ9//jnt27fHw8ODiIgINm7cyNWrV1m8eDHu7u56rW3+/Pl4e3vTtGlTvd6nJEnwCaNQvnx5FixYQFhYGJUrV6Z58+aMHTtWAlCUCs/v50tKSmL79u3079+f2rVr8/vvvzNlyhRu3brFpk2beOmllzAz0/+P75CQELZu3crChQv1fq+SJHN8wijFxMTw8ccfs2nTJgYOHMjMmTOpVauW0mUJUSQPHjygVq1a7Nixg++++44dO3bQrFkzhg8fzoABA3S2SKWwevbsSceOHZk6daoi99cX6fEJo+To6MiiRYu4fPkyjo6OeHh48NZbb3Hjxg2lSxOiUM6fP8+SJUtITEzE19cXNzc3zp8/z759+xg1apRiobdnzx7CwsKYOHGiIvfXJwk+YdQqVKjA4sWLuXz5MuXKlaNZs2a8/fbb3Lx5U+nShMjVzZs3WbZsGY0aNaJ3796Ym5vTr18/xo8fz7Rp06hRQ9knGaSkpPDuu++ycuXKYu//M0QSfKJUqFChAkuWLCEkJAR7e3uaNm2Kr6+vBKAwGBmLVLy8vGjatClXr15lw4YNXLt2jSVLltCnTx+DeTDtxo0bqVmzJr169VK6FL2QOT5RKkVFRbF8+XI2b97M0KFD8fPzU/y3aGF6kpKS2L17N/7+/uzbt48uXbrg4+ND9+7dsbLKerpJREQEbdu25datW4puG4iOjqZ+/focPHhQ7ytGlSI9PlEqVaxYkY8++oiQkBBsbGxo3LgxEydO1PthwEKkp6dz8OBBxowZQ/Xq1Vm/fj29evXi+vXrBAYG0rdv32yhB1CrVi1UKpXie1U//PBDhgwZUmpDDyT4RClXqVIlli9fTnBwMFZWVjRq1Ih33nmH27dvK12aKGWCgoKYPn06tWrVYvLkybi4uHDu3Dn279/PqFGjcHBwyPPzKpVK0efzAVy4cIH/+7//Y+7cuYrVUBIk+IRJqFy5MitWrCA4OBgLCwsaNmzIpEmTJABFsdy8eZOPPvqIxo0b07NnT8zMzPj1118zn3lX2OH1du3a6eVJDQWh1Wp59913mT17No6OjorUUFIk+IRJqVy5MitXruTSpUuYm5vTsGFDJk+ezJ07d5QuTRiJuLg4Nm/eTIcOHWjatClXrlxh/fr1REREsGTJEho2bFjktj09PRXr8f3000/cunWLt956S5H7lyRZ3CJM2p07d/joo4/4+uuvGTlyJNOnT6dKlSpKlyUMTHJyMr/88gsBAQH88ccfdOnSheHDh9OjR48c5+uKKjU1lfLlyxMREUH58uV11m5+kpOTcXd3Z8OGDXTt2rXE7qsU6fEJk1a1alVWrVrFxYsX0Wq1NGjQgClTpnD37l2lSxMKy3jcz5gxY6hWrRrr1q2jR48emYtU+vXrp9PQA1Cr1bRp04ajR4/qtN38rF27Fjc3N5MIPZDgEwJ4GoCrV6/mwoULpKWl0aBBA6ZOncq9e/eULk2UsKCgIGbMmEHt2rV55513cHZ25uzZsxw4cIDRo0fnu0iluEp6uPPevXssW7aMjz/+uMTuqTQJPiGeUa1aNdasWUNQUBBPnjyhfv36TJs2jfv37ytdmtCjjEUqTZo0oWfPngD88ssvnDt3jvfff5+aNWuWWC0lvbLzgw8+YOTIkbi4uJTYPZUmc3xC5CEyMpKlS5fy7bff8uabbzJt2jQqVaqkdFlCB+Li4ti+fTv+/v6cP3+eAQMGMHz48BJ78kFuHj16RJUqVYiOjkaj0ej1XmfOnKF79+6EhITovSdrSKTHJ0QeatSowfr16zl//jyPHz/Gzc2N6dOnExUVpXRpogiSk5PZuXMnAwcOpFatWuzevTtzX+dnn31G+/btFQ09AFtbW9zc3Dh16pRe76PVapk0aRLz5s0zqdADCT4hCqRGjRps2LCBc+fOkZCQgJubGzNmzCA6Olrp0kQ+MhapjB07lmrVqrF27Vq6detGREQE27dv18sileIqieHOwMBAHj58yJtvvqnX+xgiCT4hCqFmzZps3LiRM2fO8PDhQ1xdXfHz85MANEAXLlzIXKQyceJEnJycMhepvPnmm5QrV07pEnOl7+BLTEzkvffeY/Xq1Zibm+vtPoZKgk+IInjhhRf45JNPOH36NA8ePMDV1ZWZM2cSExOjdGkmLTIykuXLl9OkSRO6d++OVqvl559/5vz58yW+SKU4PD09OXr0KOnp6Xpp/+OPP6Z58+Z06NBBL+0bOlncIoQOXL9+ncWLFxMYGMhbb73F1KlTS3QDsinLWKQSEBDA2bNnMxepvPzyy4rP1xVHvXr1+Omnn2jQoIFO2719+zaNGzfmxIkT1K1bV6dtGwvj/VMhhAGpVasWn376KadOnSIqKgpnZ2c++OADYmNjlS6tVEpOTuaHH37IXKTyyy+/MGHCBG7fvp35zDtjDj3Q33Cnn58fY8aMMdnQAwk+IXSqdu3afPbZZ/zzzz/cu3cPZ2dnZs+ezYMHD5Quzeilp6dz+PBhxo0bR7Vq1Vi9enXmIpUdO3bQv39/vS//L0menp46P7D677//Zu/evcycOVOn7RobCT4h9KBOnTp8/vnnnDx5ktu3b+Pk5MScOXMkAIvg4sWL+Pn5UadOHcaPH0/dunU5c+YMBw8eNPhFKsWh6x6fVqtl8uTJLFq0iLJly+qsXWMkwSeEHtWtW5cvvviCkydPEhkZibOzM3PnziUuLk7p0gzarVu3WLFiBU2bNqVr166kp6fz008/ZT7z7oUXXlC6RL1zc3MjLi5OZ08O+fbbb0lJSWHkyJE6ac+YSfAJUQLq1q3Lli1b+Pvvv7l+/TpOTk7MmzdPAvAZDx8+ZMuWLXTs2JFGjRoREhLC6tWruXHjBsuWLaNx48ZKl1iizMzMaNu2rU6GOx8/fsyMGTNYs2aN0c996oL8FxCiBNWrV48vv/yS48ePc+3aNZycnJg/fz4PHz5UujRFZCxSGTRoEC+88AI///wz48eP5/bt22zevLlULFIpDl0Nd3700Ue0a9cOT09PHVRl/GQ7gxAKCgsLY+HChZlHZ02aNAk7Ozuly9Kr9PR0jhw5gr+/P9u3b8fd3R0fHx8GDhxYaufriurIkSNMmjSJf/75p8htXL9+HQ8PD86cOWMSQ8QFIcEnhAEIDQ1l4cKF/Prrr0yaNIl33nmn1AXgxYsXCQgIICAggLJly+Lj48OwYcPkh3EekpOTcXR05O7du9ja2hapjSFDhuDq6sq8efN0XJ3xMt0xBCEMiIuLC9988w1//fUXly9fpl69eixatIj4+HilSyuWjEUqzZo1o2vXrqSmprJr167MZ95J6OXNysqKpk2b8vfffxfp83/99RdHjhzh/fff13Flxk2CTwgD4urqytatW/nzzz8JDg7GycmJxYsXk5CQoHRpBfbw4UO+/PJLOnXqRKNGjQgODubjjz/m+vXrmc+8U6lUSpdpNIo6z5eens7kyZNZtmwZZcqU0UNlxkuCTwgD5Obmhr+/P4cPH+bixYvUq1ePJUuWGGwAPnnyhB9//JHBgwfzwgsvsGvXLnx9fbl9+zZffPEFHTp0MMnDkHWhXbt2RVrZ+fXXX2NpacnQoUP1UJVxkzk+IYxAcHAw8+fPZ9++fUyZMoUJEyYUeM4n+lEygaciCbkbT3xSKnYaNW5V7BjUvAaOtkV/HE96ejpHjx7F39+fwMBA3N3dGT58OAMHDpRzSnUoNjaW2rVrExsbi1qtLtBn4uPjcXNz48cff6Rly5Z6rtD4SPAJYUQuXbrE/PnzOXDgAFOnTsXX1zfXADx3M44NB69wKPTpQ3OTU/876V+jNkMLeLlWxLe9E01qOhSqBn9/f7799ltsbW3x8fFh6NCh1KpVqzjfmsiDu7s7W7duxcPDo0DXz5gxgzt37vD111/ruTLjJMEnhBG6ePEi8+fP5+DBg0ybNg1fX98s8zj+xyNYtDuEpNQ08vobrlKBRm3OrB5u+LSpnet1t2/f5v/9v/+Hv78/UVFRDB06FB8fHxo3bizzdSVg3LhxNGzYkIkTJ+Z7bXh4OK1atSIoKIhq1aqVQHXGR4JPCCN24cIF5s+fz+HDh5k2bRpvv/02O4OiWLQ7mMSU/3p4aYkJxOxeQ1LEGcys7SjXfiRl3L0y37e2MGNWj/pZwi8+Pj7zcT+nT5+mX79+DB8+nPbt28t8XQn75ptv+OWXX9i2bVu+1/bv358WLVqY/EHUeZHgE6IUCAoKYv78+fwVfJMyvWeR+ty6tagfPwKtFsce7/Dk3lXuB86jis9yLCv+NzxpbWGO/xvNuRV0DH9/f3777Tc6dOiAj48PPXv2xNrauqS/LfE/V69e5aWXXiIyMjLPHvaBAwcYNWoUwcHBpepJFbomwSdEKTJk/T6OR/4Lqv+CL/1JEjdXD6HamxuwKF8dgOifVmJe1pFyXq//92FtOqkRp6h7ez/Dhw9n0KBBskjFQGi1WqpXr86RI0eoU6dOjtekpaXh4eHB7NmzGThwYAlXaFwKtkRICGHwoh8lc+bekyyhB5AaewuVmXlm6AFYVKpD8o2grA2ozLBxasXOzX7FWu0pdE+lUmU+ny+34Nu8eTMODg4MGDCghKszPrKPT4hSIvBUZI6vp6ckorLKOkxpZmVD+pPEbNeaqVQEns65HaGsvDayx8XFMWfOHFavXi2LjQpAenxClBIhd+OzbFnIYGZhjTY5a8hpk//FzDL7nF1Sajp7T16kgeo25cuXx9HRkfLly2Npaam3ukXBNGrxIhv2XWbytjPZ9mMumj+fPn360KxZM6XLNAoSfEKUEvFJqTm+ri5fHW16GimxtzKHO5/cv4ZFxZz33YVei8Rv50fExsYSExNDbGwsGo0mMwQdHR2z/Htur5UrV05Wf+rAf/sxY0h27cwPZ29nvqdR32Xl3ssk3rZnq98EBas0LhJ8QpQSdpqc/zqbWWqwcX2RuD8DcOz+Dk/uX+XfK39TxWd5jtd3eqktq9b7Zn6t1WpJSEjIDMGYmJgs/x4REcHp06ezvR8fH4+dnV2+Afn8a3Z2djJc9z/P78dUqbPOvSb9r4evrtMc3+2hzEo0y3M/pnhKgk+IUsKtih1W6rs5DneWf8WXmN1riFw3HDNrOxxf8c2ylSGDRm2GW9WyWV5TqVTY2dlhZ2eX68KKnKSlpREXF5djWMbExBAcHJxjmCYmJlKuXLlChaWjoyPW1talKjCfhl7W/Zi5U5GYksai3cEAEn75kO0MQpQS0Y+S8Vy2P8fgKygrtRlHp3dUdFXnkydPiI2NzTLU+nxw5vSaVqst9HCso6OjQc5fnrsZx5DPj5OYkpbl9eifVpAUcY70lCTMy5TDrs0AyjbpmuUaawtzto1tQ+MaDiVYsXGRHp8QpUQFWyvau1Rkb/C9PI8py41KBR1cKyq+lcHS0pIqVapQpUqVQn3u33//zTUY79+/T3BwcLb3M+Yv8+tNPv++vucvNxy8QlJqWrbX7doMwrH7JFRqC1JibnL3Wz8sK9fDqopT5jVJqWlsPHiFTT4t9FafsZPgE6IUGe/lxJ9h0dl6CgWhUZvj6+WU/4UGysbGBhsbG2rUqFHgzzw7f5lTjzKv+cuyZcsWKiwLOn8Z/SiZQ6FROf7yknV4WoUKFakP7mQJPq0WDlyOIuZRsuK/xBgqCT4hSpEmNR2Y1cOtEHNDTz09q9PN5IbHijt/mdvwa8b85fPvZ8xf5hWQ51Mqk56e+/BrzG8beRy0D21qMpaV62FdL3vPTgUEno5k3Mv1ivKfpdSTOT4hSiFdP51B6EZB5i//sXAn1i7vnrc2PY3kWyEk3QjCvs1AVObZ+zD9mlZn1atN9fSdGDcJPiFKqfORcWw8eIUDl6NQ8d/Sd/jveXwdXCvi6+Vkcj09Qzbq65PsD7lfoGtj9qzHosIL2LXok+29Tm6V+GKkPIQ2JzLUKUQp1biGA5t8WhDzKJnA05GE3EkgPikFO40FblXLMtCjeE9gF/qR237MHKWnk/rgTi7tWOiootJHgk+IUs7R1krmeoxIbvsx0x7HkXT9HNZOrVCpLUmKOMvj4ENU6PN+tjZy2o8p/iPBJ4QQBmRg8xqs+iM0+xsqFQlnfiXmt42gTUdtX4lyncZg49w626VaYKBHwVe3mhoJPiGEMCC57cc0t7GnyvCl+X7eUPZjGjJ5LJEQQhiY8V5OaNRF2yBv7PsxS4IEnxBCGJiM/ZjWFoX7EW2q+zELS4Y6hRDCAGXsq5T9mLon+/iEEMKAyX5M3ZPgE0IIIyD7MXVHgk8IIYRJkcUtQgghTIoEnxBCCJMiwSeEEMKkSPAJIYQwKRJ8QgghTIoEnxBCCJMiwSeEEMKkSPAJIYQwKRJ8QgghTIoEnxBCCJMiwSeEEMKkSPAJIYQwKRJ8QgghTIoEnxBCCJMiwSeEEMKkSPAJIYQwKRJ8QgghTIoEnxBCCJMiwSeEEMKkSPAJIYQwKRJ8QgghTMr/B9gJF4Lp6pNNAAAAAElFTkSuQmCC"
     },
     "metadata": {}
    }
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Parameterized Quantum Circuit (PQC)"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "The PQC with $p$ layers can be written as:\n",
    "$$\n",
    "U(\\vec{\\beta}, \\vec{\\gamma}) = V_{p}U_{p} \\cdots V_{1}U_{1},\n",
    "$$\n",
    "where $U_{j}= e^{-i\\gamma_{j}H_{C}}$ and $V_{j}= e^{-i \\beta_{j} \\sum_{k} \\sigma^{x}_{k}}$."
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "source": [
    "def QAOAansatz(params, g=example_graph):\n",
    "    n = len(g.nodes)  # the number of nodes\n",
    "    c = tc.Circuit(n)\n",
    "    for i in range(n):\n",
    "        c.H(i)\n",
    "    # PQC\n",
    "    for j in range(nlayers):\n",
    "        # U_j\n",
    "        for e in g.edges:\n",
    "            c.exp1(\n",
    "                e[0],\n",
    "                e[1],\n",
    "                unitary=tc.gates._zz_matrix,\n",
    "                theta=g[e[0]][e[1]].get(\"weight\", 1.0) * params[2 * j],\n",
    "            )\n",
    "        # V_j\n",
    "        for i in range(n):\n",
    "            c.rx(i, theta=params[2 * j + 1])\n",
    "\n",
    "    # calculate the loss function\n",
    "    loss = 0.0\n",
    "    for e in g.edges:\n",
    "        loss += c.expectation_ps(z=[e[0], e[1]])\n",
    "\n",
    "    return K.real(loss)"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Main Optimization Loop"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "source": [
    "# use vvag to get the losses and gradients with different random circuit instances\n",
    "QAOA_vvag = K.jit(tc.backend.vvag(QAOAansatz, argnums=0, vectorized_argnums=0))"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "source": [
    "params = K.implicit_randn(\n",
    "    shape=[ncircuits, 2 * nlayers], stddev=0.1\n",
    ")  # initial parameters\n",
    "opt = K.optimizer(tf.keras.optimizers.Adam(1e-2))\n",
    "\n",
    "for i in range(50):\n",
    "    loss, grads = QAOA_vvag(params, example_graph)\n",
    "    print(K.numpy(loss))\n",
    "    params = opt.update(grads, params)  # gradient descent"
   ],
   "outputs": [
    {
     "output_type": "stream",
     "name": "stdout",
     "text": [
      "[-0.23837963 -1.1651934 ]\n",
      "[-0.5175445 -1.4539642]\n",
      "[-0.7306818 -1.6646069]\n",
      "[-0.91530037 -1.8384367 ]\n",
      "[-1.0832287 -1.9884492]\n",
      "[-1.2398103 -2.120449 ]\n",
      "[-1.3878661 -2.2374902]\n",
      "[-1.5290209 -2.341291 ]\n",
      "[-1.6642232 -2.4328852]\n",
      "[-1.7940071 -2.5128942]\n",
      "[-1.9186544 -2.5888019]\n",
      "[-2.0382538 -2.6627793]\n",
      "[-2.152771 -2.735217]\n",
      "[-2.2620971 -2.8060198]\n",
      "[-2.3657765 -2.8749723]\n",
      "[-2.4635859 -2.942443 ]\n",
      "[-2.5571456 -3.0074604]\n",
      "[-2.6474872 -3.071116 ]\n",
      "[-2.7343643 -3.1320357]\n",
      "[-2.8174913 -3.1904984]\n",
      "[-2.896546  -3.2464304]\n",
      "[-2.971222 -3.298626]\n",
      "[-3.0411685 -3.3485155]\n",
      "[-3.1060221 -3.3945203]\n",
      "[-3.1671162 -3.4365993]\n",
      "[-3.2244647 -3.47741  ]\n",
      "[-3.2800133 -3.51378  ]\n",
      "[-3.328074  -3.5467668]\n",
      "[-3.3779154 -3.5716858]\n",
      "[-3.42378   -3.6026983]\n",
      "[-3.4665916 -3.6264663]\n",
      "[-3.5065007 -3.6452012]\n",
      "[-3.5436964 -3.6676104]\n",
      "[-3.5783873 -3.6827888]\n",
      "[-3.6107998 -3.696251 ]\n",
      "[-3.6411772 -3.710956 ]\n",
      "[-3.6697989 -3.725151 ]\n",
      "[-3.6969085 -3.739223 ]\n",
      "[-3.7227716 -3.753837 ]\n",
      "[-3.747642  -3.7637105]\n",
      "[-3.7717733 -3.7778597]\n",
      "[-3.7953677 -3.7897499]\n",
      "[-3.8185773 -3.8026254]\n",
      "[-3.8415692 -3.8186839]\n",
      "[-3.864397  -3.8288355]\n",
      "[-3.887118  -3.8470592]\n",
      "[-3.9089546 -3.8578553]\n",
      "[-3.9298224 -3.8789082]\n",
      "[-3.9531326 -3.898365 ]\n",
      "[-3.9759274 -3.9132624]\n"
     ]
    }
   ],
   "metadata": {}
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}