DC-QAOA_EN.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Large-scale QAOA via Divide-and-Conquer\n",
    "\n",
    "<em> Copyright (c) 2021 Institute for Quantum Computing, Baidu Inc. All Rights Reserved. </em>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Overview\n",
    "\n",
    "[Quantum Approximation Optimization Algorithm](./QAOA_EN.ipynb) (QAOA) is a promising hybrid quantum-classical algorithm to solve combinatorial optimization problems approximately. However, its applicability is restricted by the qubit limitation for large-scale problems and its execution time scales exponentially with the problem size. Detailed and other limitations can be found in [1].\n",
    "\n",
    "Divide-and-Conquer (DC) is a largely used technique to address similar challenges as those listed above, such as quicksort and fast Fourier transform (FFT). It recursively breaks down a problem into several subproblems of the same type until these subproblems can be solved directly.  \n",
    "\n",
    "DC-QAOA proposed by Junde Li et al. in 2021 [2] applied such technique on QAOA to solve the Max-Cut problem. The methodology presented below adopts this DC-QAOA scheme with some modifications."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Methodology"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Large Graph Partitioning (LGP)\n",
    "\n",
    "Let $G(V,E)$ be an undirected graph with $n$ vertices and $k$ be the qubit (vertex) limit, i.e. available qubit size in NISQ computers. LGP partitions $G$ into exactly two subgraphs $G_0$ and $G_1$, and at least one node is shared between these two subgraphs so that no edges of $G$ are missed. Those shared nodes that enables this separation are separation nodes. Note that number of separation nodes should be smaller than the max qubit limit. \n",
    "\n",
    "The following code shown is to define such partition."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-11T06:51:46.909218Z",
     "start_time": "2021-05-11T06:51:45.642375Z"
    }
   },
   "outputs": [],
   "source": [
    "import networkx as nx # version of networkx >= 2.5\n",
    "import matplotlib.pyplot as plt\n",
    "from itertools import combinations\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "# Partition a graph into two subgraphs\n",
    "def NaiveLGP(g, k):\n",
    "    E = list(g.edges)\n",
    "    E = [(min([u, v]),max([u, v])) for (u, v) in E]\n",
    "    E.sort(key = lambda tup:tup[0])\n",
    "    V = list(g.nodes)\n",
    "    V.sort()\n",
    "\n",
    "    counter = 1\n",
    "    while counter < k:\n",
    "        num_nodes = counter\n",
    "        nodes_combo = list(combinations(V, num_nodes))\n",
    "        # Check suitable separation path\n",
    "        for p in nodes_combo:\n",
    "            V1 = [x for x in V if x not in p]\n",
    "            E1 = [e for e in g.edges if e not in g.edges(p)]\n",
    "            G1 = nx.Graph()\n",
    "            G1.add_nodes_from(V1)\n",
    "            G1.add_edges_from(E1)\n",
    "            S = [G1.subgraph(c).copy() for c in nx.connected_components(G1)]\n",
    "            if len(S) == 2:\n",
    "                # Add the seperate paths to two subgraphs\n",
    "                V_S0 = list(S[0].nodes)\n",
    "                E_S0 = list(S[0].edges)   \n",
    "                V_S1 = list(S[1].nodes)\n",
    "                E_S1 = list(S[1].edges)\n",
    "                for (u,v) in g.edges(p):\n",
    "                    if u in V_S0 or v in V_S0:\n",
    "                        S[0].add_edges_from([(u,v)])\n",
    "                    if u in V_S1 or v in V_S1:\n",
    "                        S[1].add_edges_from([(u,v)])\n",
    "                    if u in p and v in p:\n",
    "                        S[0].add_edges_from([(u,v)])\n",
    "                        S[1].add_edges_from([(u,v)])\n",
    "\n",
    "                return S\n",
    "\n",
    "        counter += 1\n",
    "    print(\"G has connectivity above k\")\n",
    "\n",
    "    return {}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One example for illustration is given below. The graph is a random with $10$ vertices and the qubit (vertex) limit is $9$ vertices. Then red vertices are those selected to be removed so that the rest is disconnected and can be partitioned. We then add separation nodes with incident edges back to these two subgraphs to avoid missing information of the graph. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-11T06:51:47.579620Z",
     "start_time": "2021-05-11T06:51:46.912332Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1MAAADnCAYAAAD7CwxiAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAABKwElEQVR4nO3dd5hcdfXH8fcnvRBCr0rviFQLvWNCl04ggFKkKr0oSO8ISBeB0Kt0aSKCVGkiICDlh9RQAyFl03N+f3xvJNmd2TKZmTvl83qefZSdO3fOZmfv3PMt5ygiMDMzMzMzs67plncAZmZmZmZm9cjJlJmZmZmZWQmcTJmZmZmZmZXAyZSZmZmZmVkJnEyZmZmZmZmVwMmUmZmZmZlZCZxMmZmZmZmZlcDJlJmZmZmZWQmcTJmZmZmZmZXAyZSZmZmZmVkJnEyZmZmZmZmVwMmUmZmZmZlZCZxMmZmZmZmZlcDJlJmZmZmZWQmcTJmZmZmZmZXAyZSZmZmZmVkJnEyZmZmZmZmVwMmUmZmZmZlZCZxMmZmZmZmZlcDJlJmZmZmZWQmcTJmZmZmZmZXAyZSZmZmZmVkJnEyZmZmZmZmVwMmUmZmZmZlZCZxMmZmZmZmZlcDJlJmZmZmZWQmcTJmZmZmZmZXAyZSZmZmZmVkJnEyZmZmZmZmVwMmUmZmZmZlZCZxMmZmZmZmZlcDJlJmZmZmZWQmcTJmZmZmZmZXAyZSZmZmZmVkJeuQdgBUmMSewO7AaMDswGngDuCqC9/OMrbMkugE/AbYB5gMC+AS4A3g4gqk5hmdmNUpiADAEWAuYE2gB3gWGRfBGnrGZWfOSELAmsDOwIOk++nPgfuDuCCblGJ7lRBGRdww2HYnvA78GtgKmAv2me3hi9r2ngdMieKT6EXYsuxHaHzgU6AsMmO7hAMaSksPfAZdFMLbqQZpZzZFYHDgK2JV0res/3cOTgMnAv4HTgbsi8AeYmVWcRC/g58CRwNyke7PpV3eNBqYAFwHnRzCi6kFabpxM1RCJnYArgd5A9w4ObwHOBX5bSzcUEgsCj5FGbPp2cPg44D1ggwg+rWxkZlbLJDYC7gT60PGqibHArcA+EUyudGxm1rwkBgIPACsy4wB3IeOBkcB6EbxZ4dCsRjiZqhES2wPX0HECMr2xwO8j+E1louoaibmAfwHz0vklpJOA4cDKEXxdodDMrIZJrEO6WenoRmV6LaTka2gtDSiZWeOQ6As8AyxDGujujKnAN8AqEbxXodCshrgARQ2QWBK4mq4lUpCWwBwssWnZgyrN7aTp767sxesJzE8aZTazJiMxB/BnupZIkR2/NbBPuWMyM8tcAixF5xMpSPfWswJ/yfZYWYNzMlUbDiYlFW1cdx0MHw7ffANvvgl77tnmkH7A8ZUNr2MSywE/AHpN//1eveCKK+C992DUKHjpJRg0qM3TewFrSSxRlWDNrJbsSZFlzZ24/vUHjvMNi5mVW1YIbCcKDHQvsww88giMHAlvvw1bb93m6d1JA8XrVzhMqwFOpnIm0R/YgyLJ1OmnwyKLwMCBsOWWcMopsMoqbQ77vsQyFQ20YwdT4Gfo0QM+/BDWXTf9DMceC7feCgsv3Ob53YGDKh+mmdWKrOLnoRSZlerk9W8gvmExs/L7ObStOty9O9x9N/z5zzDHHLDPPnD99bDkkm2e3x84ogpxWs6cTOVvWwr8sU7z+uswcWL6/xHpa/HF2xzWHdivQvF1SKI3sAsFlve1tMCJJ8L776fY77sP/vtfWHXVNqfpCewpuVy/WRNZlxkr9s2gk9e//qTBHDOzcvoVBQZ6llkGFlgAzjsPpk6FRx+Fp56CoUPbPF/A+hLzVCFWy5GTqfwtSTs3EwAXXwxjx6ZlLp98Avff3+aQnsD3KhRfZ3T6QjHPPLDUUvDaawUf7k7qqWVmzWFJOqhc2onrn4BlKxSfmTWhbOnw/F04nu8VvgsbD7Rdi2MNxclU/uaA9tf7H3AADBgAa60Fd9wBEyYUOuqFDSRFHl+w7AcwqsPN4z16wA03wDXXpBujAqYwY08qM2tsA+igYE3nrn/tD0iZmXVRHyhcJfTNN+Hzz+GII9J9zcYbp60M/YrfBfm+psE5mcrfCIr8wU5v6tQ0jfyd78B+BRf0rfZIRCiPL3jjuzBrS3vxS2kz+cSJcOCBRQ/rDozq6N/CzBrGKOi4T1TH1z/GlDswM2tq4yky0D15cio4sdlm8OmncNhhaS/4Rx8VPZfvaxqc96fk703SjUCnRi569Ci4Z2AS8HJ5w+qSzzs64MorYd55YdNN04WoiEngXlNmTeQtmBqdHdcrfP2LAF7vYILfzKzTIgiJj4CFCj3+6quw3nrf/vdTT6VVNwX0BveaanSemcrfHRS5C5h7bthxR+jfH7p1g002gZ13TuU4W5kCXFrhOIuKYCKpT9akQo9feiksuyxssQWMH1/0NBOBP0YwpSJBmlnNkNRb0rbQ7RD4tODimM5f/1qAQStKOlrSdysevJk1i/OAsYUeWGEF6N0b+vZNM1Pzzw9XX93msAAejuDLyoZpeXMylbMIxkFcCVPazNdEpCUtH30EX38N55wDBx8M997b5jQvRvBONeJtxwUUWK6z0EKw776w0kppOnz06PQ1ZEjrI6cGcFEV4jSzHChZXdKlwHDgQIi7od+vSRnRDDp//ev3AfxtKLAY8LKkv0raXZL3KZjZzLiaIgVyhg5NBXE+/xw23DDtm5pWeXQ6Y4FzKhui1QJFdLhdxypI0tyw4jXw9CDoV8o6lRbgpxH8pdyxdZXEX4G16FqncGDSZPjbVBi0L3B1+E1p1jAkLQbsCgwlzaJfC9wQEe+nx5mNtAxmYAmnHwscFMGw7LX6AFsAuwFrA/dmr/e3iPCst5l1icQfYPIe0KNX1545ZSp0+w/oexEd74u3+uaZqRxJ2gZ4BV5+DSZsD4zr4ilagDNqIZHKbA98QpHlfkVMhJ7vwanrk3o63CNpgUoEZ2bVIWk2SXtLegJ4FpgbGAIsGxGnTUukACIYCQymwOxUB1qAm0mjx9m5YnxE3BYRWwBLAc8DpwMfSDpL0goz8WOZWRORJBj4PrwKTC1cR7SgmAKjp8Jyb4O6mIRZPXIylQNJc0i6ATgT2C4ijoiY/XZS49sWOq5uFdlxJwOnVDTYLojga+DHwNt07saoBXgDWD3i8aeBHwIvAf+StEu6kJlZPZDUU9Lmkm4F3gd+ApwNLBgRB0XE88VmnSN4BhgEjCbtn+zIWGAYsE+xUd+I+DwiLoiI1YCNSbNi90t6SdIhkubr8g9pZk1BUnfgIhi1I5z8Pej2DEX2T7UyDjQcPloB/jMZeEjSbBUN1nLnZX5VJmkz4HLgT8AxEdEy4+MsAxwB7AxMZcb+KdPKN/yVNCP1VOUj7jqJfsCepJ9jdtLPMC0xmkpKor4EzgKGRTB+xudrNeAaUqXDfSOiw2qBZlZ92YDHKqRldTsB75CW1d0aEV2uzCmxEHAIsFf2rVmme3gi6frxLHBmBA+UEG93YN0s3q2AZ7J4746Irq4MMLMGJKkfcCPp3mXbiBgl0YM04H0ksAhpO8P0+6lGk1YXnQdcFsFISd2A35EGcwZHxIfV+ymsmpxMVYmkgaQ/svWBn0XEY+0fz6ykZTErA3MC35CSi+si+KSy0ZZH1kF8XWBLYNoo8KfAncCT7a0jzvY+nADsARwUEbdVNFgz67Ssat4upKSkDykhuT4iylIIR6Ivadnwj0lLBMcC75Kuf/8tz2uoP7A16Wf4Aem6dC3wRERMLcdrmFl9kTQXaa/l28BeEdFmplxiVWAHYEGgJ+m+5i/Ag4UqEks6BDgU2CwiXqlg+JYTJ1NVIGkT4ArgPuDIiBidc0h1Q9KPSbNULwEHRMSInEMya0pZdbxtSMnHSqTZ9euAp+q9aEy2T3MIsDup59/1wHUR8WaugZlZ1UhaHHgAuBU4rpzXNUnbAxcDO0dE2wY3VtecTFVQdvNxNrApsGdEPJxzSHVJUl/gVNIyon0j4p6cQzJrCtmyuA1Jlfi2AB4nJVD3RkTxrnF1Klu2uCLp5x0CfECarbrZAzlmjUvSD4C7gRMj4g8Veo21SYNQh0fEdZV4DcuHC1BUiKT1gVdIU8ArOJEqXUSMi4hDgR2B8yRdK2n2vOMyy4PEIhLnSHwuMVFiisRoiQclNsyW187ka2gFSWeRkolTgeeAJSNiy6xaXsMlUgCR/CsiDgO+CxwPrAn8n6Q7JW0jqYutH9qS+K7E6RKfSkzIfodjJB6R+Inkz2azapG0OWnl0C8qlUgBRMQTpK0eJ0s6xkW2GodnpsosW4d/Omk5zC8i4r6cQ2ookmYhVUHcEtg7Ih7MOSSzqpBYGLgKWIM0ENa65G6Q9hZ9A/wygju6dn7NRyp8sxswF98udXt9JkOve5JmBbYl/dusANxGmrH6R1eWAkksSPodrkMqylMoMRtN+j0eGsFNMxm6mbVD0i9I+7O3iojnqvSaC5CSt3+Q9oR3VMHZapyTqTKStCap58mzwC8j4qt8I2pckjYErgQeBg6LiFE5h2RWMRLfBx4lNbbt3sHhkCpmnhDB2e2fV31JVe12A1YnLXO5FnjMRRgKk7Qw3xbf6E5a9nhdRLRbGENiWdIyydmAHp14qRbgrAhOnKmAzayNbFboZNKKl8HlKp7ThdcfQFryN4G0j6ozZdetRjmZKoPshuRk0gfs/hFxZ84hNYVstHha2dE9vanTGlFWLvwlUpuBriwLaSHNUF054/nUDViblAxsQxr8uQ64yx/onZfdjK3Gt2Xh3yD9O94WESNnPJYFgH+RZvy6+js8KoKLyhCymQGSegF/BJYGtoiIL3KKoyepVc5yWRxuA1OnnEzNJEk/JFWbe5VUbS6XP8pmJmkQ6cJ4D3BURIzJOSSzspH4C7ABnZuRam08sFAEX0hamlRYYVdgFGkG6saIGF62YJtUdnM2mPTvuzHwICmxeigiJkncBWxG52akWhsPLBHBx2UK16xpZYOwfyL1hNq5da/PHOIRaZnhLsCgas+QWXl4k2uJJPWWdCqpH8EJEbGDE6l8ZPumViA1+HxZ0jo5h2RWFhLfJc0iFU2kllgCxo2D6wrWhpoacMcfJT0LPAb0Je0N+H5EnONEqjwiYmJE3B0R2wGLkpZkHgN8JC19OUwdTKtEqlcvuOIKeO89GDUKXnoJBg0q+hL7VjB8s6aQ7VV6nNRcfJu8Eyn4X9Gb40l7wZ+Q9KO8Y7Kua9qZKYl5gJ+TmkLOTtrw+x/gjxG80f5ztTJpNupdUqnuTyscrnWSpC2BS0l9In4dEeOKH0svUtPOnwLzAlOBT4CbKdJ8z6yaJE4HDqFwoQIAHnoI+vaF99+HoUMLHTFyAiy4DbT8xRudq0vSEnDT5bDVetB3huV9/frBEUfA1VfDBx/AppvCTTfBCiuk32UrXwPzRjCpOpGbNRZJy5OKPlwGnFmLvfEkbQYMIzULdguYOtJ0yZTEysBxpCUZQRqpnWYyMAl4DTg1grtmfK56Ar8GDgAOA66vxT/IZidpTuAiYBVg94j4x4yPMwdwOLAfacR/QKtTjCYtrTkPuCAC7yOxXEgMB+Yv9viOO8I228Drr6cZqsLJFKOBjSN4tkJhWjsk3iXNVnXo5ZfhxBPhjrZ1GEcBW0XwWHmjM2t8ktYlDbAeFhHX5x1PeyStRtqycHJEXJp3PNY5TbXMT2Jn4ElS9ao+zJhIQVqG0Ze0qfh6iYultLxG0vdIZSx/BKwcEdc5kapNETEiInYGfgPcJelMSX0AJJYg9f86lFRVq3UiRfa9uUlJ9wsS81UlcLO2Ziv2wIABcNJJcOihHZ5jKjBPGWOyrpmjMwfNMw8stRS89lrRQ+YuW0RmTULSjqRWBkNqPZECiIgXgLWAgyWd7l5U9aFpkimJ7UiltPvRuZ+7P7AHTLlE6n4MaQ38JcBmEeGNwHUgIv4EfB9YAnhR+vlg4BnSSH9nGm/2zZ77jFT8ptasgorulTr5ZLjySvi446uRSM3DLR8dFp3o0QNuuAGuuQbefLPgIaJtXzEzK0LJocDZwEb1VO03It4lNQtfF7g2K3BjNawpkimJRUl7nFrPRHWkH0zcC/YdAqwWEVd6Nqq+ZKVGt4Nup8Jh98CUOeja+74HKfm6oSIBmrWv4BLTFVeEjTaC887r1DmCtOfG8jG6vQelVDxk4kQ48MCih03Fv0OzTpHUnbRM/2fAGhHxSs4hdVlEfAlsSBrYv1/SwJxDsnY0RTIF/JJ2Rgfbr4bVtxtc1D0i2m4JtrqQEuAp/wfLToTuM7znDzgAnn8exo+HYcOKnqI3sIHEIhUO1ay1v5NupGew3nqwyCKpcMEnn8Dhh8O228KLLxY8Ry/gnxWN0trzCBQvZnPllTDvvOn3N7l4eZA+wHMViM2soWR9P28BVgTWjoiPcg6pZFkBre1JPeyekPSdnEOyIho+mZLoC+xFO0skLr443VC3c5aFJVYrd2xWVYdBtz6tvzl8OJxyClx1VYfPF7B/JQIza8c5pH4oM7j8clh8cVhppfR12WVw333wk5+0ef5k4NYIvql4pFbMucCEQg9ceiksuyxssUUa0CliKnBPBF9WKD6zhpAVn3oYmEjq2TQy34hmXkRMIU0IXAc8ne3ftxrT8MkUsCVpmUtBO+4II0fCI+2vpu1NquBndUhiAOl90Ob9fuedcPfdMGJEh6fpDewr4c2gVk1PA5+1/ua4cfDZZ99+jRmTbsa/bHu7PZG03MVyEsE/gfdaf3+hhWDffVMy/OmnMHp0+hoypM0pxgG/q3igZnVM0qLAU6QiY7tGRMEBjHqU9aI6GzgK+Juk9fOOyWZUSjf2erM4qehEG9OqYW2wAey1V7vn6A4sW4HYrDrmJ91UdqboRHv6kNYvj5npiMw6IYKQOAj4E+3s+TzxxILfHkfql/ZyZaKzLvglqcH7/36HH3yQ9kt1YBypyaiX+JkVIWlVUjnx0yLi4rzjqZSIuEnSp8DNkg6OiJvyjsmSZpiZGkCRilhdqIY17TxWn2ahndnJLpiM3wdWZRHcT+qL1tKFp40DXgbaznNY1UXwCGl1Q9Em4gWMIzWS3y6iLNcvs4YjaTDwAHBAIydS00TEo6TCFGdIOsKl02tDM8xMfUO6CZ7hZ51WDWvllTt9nlFljsuqZzTlGTjoid8HloMILpEYATEMxvdOhXEKmkyahb0f2DWi8F4dq74IhkmMBK4nDe70L3zkJKDbeOj+CLBDRJcSMLOmIWlP4BRgq4h4Ju94qiUi/i1pDdJ1fqFslqpokRurvGZIpt4gjfDNMKMwfTUsgFlmge7dYbnlYNVV25xjMvBSpQO1ivmYdvr1dMEoujY7YFY2EdwirTIZNr0QTpkKGsiMVeJ6kUr4nx9B8davlpsI7pSYHxgKHAHMRfp8maYXvPwqHD464rEtcgnSrMZlszHHk/6O1o2It3IOqeoi4mNJ6wC3A3+SNCSr/mc5UKO3TZLoSdrAPfv03+/bF2ad9dv/PvzwlFztt1/BTdzjgFUjeKOiwVrFSFwJ7EarAYTu3VPDzOOPh+98B/beO5UnntJ2jGcccGYEhXenmFWBpMeAP0LcSCr9Oy8pifoaeDmi/Z5GVjuyYjYrAPOR9nOOBF4BTQL+S7pJ/E9+EZrVHkk9gT+Q/nY2j4g2BXqaSdbQ9ypgMWDLrD+VVVnDJ1MAEieRRgHblMae5vjjU7+poUMLPvx8BD+sUHhWBRIrAM/SahP/8cfDCSfMeOwJJxTc0D8eWDSCTysVo1l7JK0O3AgsGRHFOxJZ3ZN0HLBIROyZdyxmtULSAOA20mzujhFRsKl5s8lm6k4FtiOVhH8355CaTrMkU/MDb5EKEXRVC2kD8APljcqqTeIx4Md0varfOOCOCHYte1BmnSTpLuDhZthk3ewkzQG8DaxYz01HzcpF0vzAfcALwP4eUGpL0n7AcaQ9ZO12T7XyaoZqfkTwCanPUFf3u7SQlnY5kWoMPwU+LbiIr7jxwJukxs9muZC0HGkgoOP20lb3IuIrYBhwaN6xmOVN0jKknnu3A79wIlVYRFwK7AvcJ2mzvONpJk2RTAFE8CgpoRpDqnbVnqmkROok4OQKh2ZVEsHXsOUB8HbA1M5s1BwLPA+sG8H4Codn1p4jgQu8wbipnAfsIWnOvAMxy4uktYDHgBMi4tRohuVUMyEi7gG2AK6QtHfe8TSLpljmNz2JhYFfQcsvoeck6Dn9PqpxgEjlJs+K4NlcgrSKkDQf8AJ85yD4cB7SDeq8pKbO03o1TIHxggkfwsDjgJsjmJRTyGZIWohUTXTxiBiZczhWRZKuAD6IiJPyjsWs2iRtC1wK7BoRf8k7nnoiaUlS/60bgeOdhFZW0yVTAJLmhT5vwvBDYfblgLlJ/aj+D7gxgi/yjdDKLasA9FfgsYg4Pn0PAWsAg4D5STOSw2FoL7h+4YjYJbeAzTKSfg9MiIgj847FqkvS0sATwKLebG/NRNKvSIXDNo+If+UcTl2SNA9wL6lF0N4R4YHhCmnWZOpA4EcRUbh2nzUcSecAy5MuzO3umZI0N2nz94K+gbE8SZqLVDznexExPO94rPok/Ql4PCIuyDsWs0qT1A04GxgMDI6I93MOqa5J6g/cRKpmvV1EjMo5pIbUNHumWhlCenNZE5C0A7ANsEtnuoRHxBfAM6R1x2Z5Ogj4kxOppnYmcFg2u27WsCT1Id2brQas6URq5mUDwtuQVl79XdICOYfUkJoumZK0KLAk8HDesVjlZVXQLga2zSpkddaNpKTbLBeSZgH2J43SWpPKShy/DeycdyxmlSJpduCh7D9/EhFf5xlPI8mqH+4P3AI8nd0XWRk1XTIF7EQa6fXa0QYnaVbgDuCIiHipi0+/C1g36/diloe9gUcj4u28A7HcnQEclS2BMmsokhYGniL1kNo5Ilw9t8wiOYPUh+pRSWvnHVMjacYL8xDSrIM1sKwj+DBSwYmru/r8iBhNGiXbtsyhmXVIUi9Sj6Ez847FasIjpGqzm+cdiFk5SVqJlEhdHhGHRcTUnENqaBFxHbALcLuk7fOOp1E0VTIlaQVgIOkP1xrbEcB3gF/NxDluwkv9LB+7AG9ExIt5B2L5y8oanwEckw0UmdU9SZsAfwEOjojzcw6naUTEX4GNgXMlHZJ3PI2gqZIp0przmz3y0dgkbQAcQqpcM2EmTvUAsKKkBcsTmVnHsqVcR5Funs2muROYE/DyHKt7knYHrgW2iYg/5R1Ps4mIl0mtYfaUdL6k7nnHVM+aJpnKRvN2xkv8Gpqk7wI3kCr3fTgz58rWbd8F7FiG0Mw6aytgFPBo3oFY7cgqkZ4FHJ13LGalUnIscAKwXkQ8mXNITSu7R1oLWBG4RVLfnEOqW02TTAE/BsYDL+cdiFWGpN7AbcB5EfG3Mp32RlxFy6okG/Q5GjjDHeutgOtIs+Ur5R2IWVdJ6gH8AfgpsHpE/CfnkJpeRIwEBgGTgIdddKs0zZRM7Qzc6BuUhnY+MJzylpJ+FPiupCXLeE6zYtYj7eu8K98wrBZly5bPIy0DNasbWfPYu4CFSDNSn+YbkU2TXVd2AZ4GnpK0SL4R1Z+mSKay0ZAdcaPehiVpD2ADYI9yJszZ0ppb8OyUVcfRwFne12ntuBzYWNLieQdi1hmS5gUeAz4Htsiq5VoNiYipEXEkcAkpoVol75jqSVMkU6Sb7Pcj4p28A7Hyk7QyaTbqpxExqgIvcSMwxFW0rJKyD6/lSXv+zArKrnGXAYfnHYtZRyQtRZrxuA/Y0z0+a1tEXAgcCDwoaVDe8dSLZkmmXHiiQWXre28HDoiI1yv0Ms8BPYGVK3R+M0hLt86dyQqU1hwuAHaUNF/egZgVI2l14O/AaRFxgrdZ1IeIuBPYGrha0s9zDqcuqNHf25L6AJ8Ay0fE8LzjsfLJSkjfB7weEYdV+LVOAXpHxBGVfB1rTtmevKeBxbwExjpD0kXA6Ig4Ju9YzFqTtDVpSeruEfFAzuFYCSQtDdxPKmF/kpPh4pphZmpT4J9OpBrSb4F+VKdU8E3ATlkCZ1ZuhwOXOpGyLvgdsI+kgXkHYjY9SQcAFwODnUjVr4h4k9SLagvgj5J65hxSzWqGG8MhuPBEw5G0GbAXsGM11mBHxGvAV6SeDGZlI2l+YHvgwrxjsfoREf8lNRbfN+9YzCCtFpF0BnAQsFZEvJh3TDZzIuIzUpXZBYC7Jc2Sb0S1qaGX+UmaFfgQWCQivs47HiuPrIrV06SCE09X8XWPJr2XfPNiZSPpTKBvRPwy71isvkhaAfgLsGjWZNwsF1mfx6uARYAtI2JEvhFZOWVVsS8FVgE2c2n7GTX6zNRPgcecSDUOSf1IBSdOrmYilbkZ2FZSryq/rjUoSbORZlh/l3MoVoci4lXgBWD3vGOx5pVdxx4A+gAbOZFqPBExGdiH1Cvs6Ww/lWUaPZlyFb8GkpUmvwz4N2k9dlVFxHvAW8DG1X5ta1j7A3+OiPfzDsTq1hnAkdnIsVlVSfou8ATwKrBDRIzLOSSrkEhOBk4C/i5pzbxjqhUNm0xJmgf4MXBv3rFY2ewLrAjsk2NVmRtJ+/DMZoqkvsAvgbPyjsXqV0Q8BQwHts07Fmsukr5PWnI/DDg4a3JvDS4iribNht8pydcdGjiZIm3o/nNEtOQdiM28rF/FicC2Of9ObwM2k9Q/xxisMfwMeDYrbmI2M84AjnZjcasWSRsAfwUOj4hzXTa7uUTEQ8BPgAskNf1+30ZOpobgJX4NQdK8wK2k7unv5BlLRHwO/INUKtSsJNmSrCNIN8FmM+t+oAfp5sasoiTtQqqSvH1E3JJ3PJaPiHgJWBPYT9LZzdw6piF/cEmLAEsBD+ccis2k7KbzZuDqiKiVJZs3kfbjmZVqB+CDiHgm70Cs/mWzAmdQnZ571qSUHA2cBmwQEX/POybLV7aXfE3Stpobs6qOTachkylgJ+BP1eg/ZBV3OjABOCHnOKZ3J7CepDnyDsTqT7YU62g8K2XldQuwcLYk2qysJHUnFX7aCVjdy5Ntmoj4ilSYqzvwkKTZcw6p6hoimZLoLTG7RPfsW27U2wAkbQdsB+xSSxtbI2IUqbfLNgBIfZFmo4mnuK0wiW4Ss0n0ne7bg4EAHswpLGtAWenic4Cjpn1Ponv22diUo8XWDklI/ZAGdvTZlbUkuQNYElgnIoZXJUarG1mfux2BfwJPSlqoveMlJDGLxKwSdb/Xs25v/iQWlThX4hugBfgUmCRN+BgO+y6893LOIdpMkLQsqUHcdrXYs+IX8PRdcDLSOGA08BkwCelVpF1p0qluA4k+EkMl/g1MIr03RkuMk7gKfnwScIY3bFsFDIPZ1pBeOUXiPdL771OgReIbid9JLJJrhJYvaWWk64FxwCjgc9Jn10tIO9Gqj6KkuYG/ASNJzVpHVTtkqw8RMTUiDgWuIPWiWmn6x7MEan2J+4CJwNfAl8AkiUclBkv1mZeo3j7PJeYBbgDWIiWDBRqoTpgEvScDlwBHRVAzsxrWMUkDgOeAsyPiqrzjmYG0FHBLwNKToW/PwkeNzv73OOAC6u2PzEqSja4dwrdLUge0PWrqFJjQDXr/C7rtEEGuBVWscUj0AM6GSQfCFKBPob5TE0izok8AQyL4spoxWo6k5Uj7jxcHesP/VvJMbzTp/XE0EZdKWoLUjPcW4DgPAFlnSdqetCx0l4h4WGI94BpgDqA/FJyNGgOMBfaK4M/VirUc6iqZykbUngbmAorcx86ghfShsUUE3j9VB7L9JLcBIyLiF3nHMwPpB6RSsLPQuVndFtLF4wAnVI0tS6QuA3YF+nXiKVNJNy4bRvBiJWOzxifRi1TRb3U69/6bSBoRXj2CDyoZm9UAaQ3gIYrfxLbW8hbctQxsEHBCRPyhsgFaI5K0NvAnuOg2OODnMMNy9/aMAw6N4LLKRVdedZNMScwBvAQsSOERlWJagHtIo3D18cM2MUmHkyqdrR0RE/KO53+kxYEXgYFdfOZY4CwiTip/UFYrJE4mzUp1tf/YSGCVCP5b9qCsKWSJ/C3A5nT+ZgXS9NWHwMoRjKxAaFYLpGVIKz0KzJQXNxbiWbh6g4ifVyYwawbSHbvD4GHQt6v7osYBu0RwZyXiKrd6Wpv4W2BeWiVSjz4K48bB6NHp6z//afO8fqQPmQ2qEqWVTNL6wOGkfVK1k0gll1How2j22eGOO2DMGHjvPdi5TcX0/sAxpHL91oAkFgcOo0Ai1Ynr0wDS3kCzUm0MbEqrROqAA+D552H8eBg2rODzugPzk5YjW+O6gkKDPB18dvUHbQA7Iy1YnTCt0aT9T9ucXiiRWnhhuO8++Oor+OQTuPBC6D7jNElf4Op6KZ5TF8lUVglrTyj8j3rggTBgQPpaZpmCp+hPukm3GiXpO6S9cLtGRG0tO0lVadam0N/LxRfDxIkw77ywyy5w6aWw3HJtzgAcWPlALScH0c5seQfXp+7AehK+YbFSHUGBpX3Dh8Mpp8BV7e867Q3sXS83LNZFac/TqpT+2QWwb4WjtMa1MWlbRBuXXAKffw7zzw8rrQTrrgv779/mMAHbVjbE8qiLZIq07GtmlugJ37DUrKzJ223ABRHx17zjKaDtnzhAv36w7bZw3HEwdiw89RTccw8MHdr6yN7APq7w13imG+gpUAin0wLYrzwRWTORWIhUjKnNyO+dd8Ldd8OIztVC3a7MoVltKDzQ0/nPrj7AAUid2aNu1tqRFEmmFl0Ubr0VJkyAzz6DBx+E5Zdvc9gApmv1UMvqJZkaSjvrfU8/Hb74Ap58MmW3RUwFtqhAbDbzziWV7z0z70CK2JlCs6JLLQWTJ8Pbb3/7vZdfLnhFIN0wr1Wh+Cw/a5OuLUV14vrUB9ilArFZ49uCmRtohPTZulsZYrHasyOFinV17bOrO/DDCsVnDSorirMuRQqenH8+7LQT9O0LCywAgwenhKqApbMq3jWtXpKpeYs9cNRRsNhisOCCcPnlcO+96b8L6EOqAmg1RNJuwEbAHjVcdrVwN+9ZZoFRrVpufPNNWs9VmN9/jWdu2qmO1YXr02wVis8a29ykz7aZVfM3K1aSwgWTuvbZFfizy7pudiheRfvxx1PuPmoUfPwxvPAC3HVXwUMnUgfvv3pJporG+dxzaf/kxIlw7bVptnrTTQseqvbOY9WXNXT7HbBNRHyTczjtKfy+GTMGZp11xu/NOmuqNNCW33+NqRvtJFNduD75vWGlKNf7xu+/xuTPLstLN4rMmktpFuqOO6B/f5hzzlQP5czCa5OCOnj/1XyAmc6t+gYi0i+qgAmkbstWAyTNDtwOHBQRr+UdTwcKfsLw1lvQowcsscS331txRXit4I8zFb//GtFX0Pmm4O1cn0YV/K5Z+74ijdyW4zzWeMYU/G7XPrsCf3ZZ131NkaJxc8yRqvlddFEaaPzqq1RxtMhAYy/q4P1XL8nUnaR+PTMYOBA22QR6904lFYcMgXXWKbruMoCHKxyndYKkbsD1wD0RcXPe8XTCfcDkNt9taUlDKyedlDb0rrEGbLUVXHddoXP0IjWctsbyFEU+MLpwfZoE9dXt3WrGwxRJ5rt3//a9N/3/L2AscEcFY7T8PEih90fXPrt6As9XOE5rMBGMB14p9NiIEfDuu7DffumaNHAg7L47vFLwaD4Dhlcw1LKol2TqagpUpOnZM5V+/eIL+PJLOOgg2HrrGfdUTufVCNp2ebE8HEva9Hxk3oF00vkUW/u7//5pB+Xnn8NNN6Wrw+uvtz5qMnArEZ59aDBZs9PbKXDD0oXr02Tg95WO1RpPBK8BbxZ67NhjU4+pY45JRdrGj0/fK6AbcG0Fw7T8nEtaldNW5z67JgHXENFmMNusE86kyMqebbaBQYPS5+M778CkSXDIIW0OGwucEzHTRXYqTrW7539GEtcDO9FOP5d2jAH2iOD28kZlXSVpMPBHYLWI+DTveDpN+iewconPbgHWIOLlMkZkNUJiFeAJCvT66aTnIvhRGUOyJiKxI6kxa8ESxB2YDNwQwR5lDcpqh/Q6sGyJzx4HrExEwYTdrD1ZRb8vgFk7OraIccB8EbW/DL5eZqYgzWYUXv/bvgnAa8Dd5Q3HukrSYqRZxp3qKpFK9iclRV3VAtzpRKpxRfBP4F5Kf3+4obPNjDvgy0+LTUB0YAxwfJnjsdqyL+mmtKtagBucSFmpIpgIHxwD49ttH1LEWODYekikoI6SqQjeAzYhTRl2djptPPBfYFBEgT0vVjWS+pKWQ50aEU/mHU+XRfwD2JWufSi1AP8AflaRmKyW7EbaV9CVhKoF2DnC+xFsZmgvWHoWmPoh6TOvM4L0WbpxBO9XLjbLXcTjpM+grn52PY6bidtMkLQ4LHww3P4ERFc/G68CzqtMZOVXN8kUQATPAasDH9JuUjVlKunC8Xfgh9m+BsuJJAGXAq8DF+YcTuki7gQ2A76OYhX+kgmkm5qbgZ8QUbTXgjWGNALHxsCtpN99O9MEY6dCfA0MjuCeqgRoDUdSN0mnAwfDV2tC3+8BT6bdB5OLDThOS6I+AH4UwQtVCtfyFHELsDXwDe1/do3Pvq4BNifCg9BWEkk/IC1/Pzdil/VAvyQlSe3tv2shvf9OAH5VD3ulpqmbPVPTk+hGunE5EliTdOOS1aKf0gtuboEd1oro2WY3pVWfpH2BA4AfRyNsZJV63QzHfR8OXS4NSEwkvf+6k0qgXwpcTMSHeYZp+ZBYmLQsdF/S+2MKqVdLL4jXYI954N5fRHz1QJ5xWv2S1Is0crsosGVEjMi+3wNWehfueBkW3ZC0J2oq6f3XmzTbcDbwSASlLL2xeib1BrYBjgaWYsbPrsnAxcClRHycW4xW9yRtTro+7RkR9377fQaQVvgcRWo4Pm2guQepPcg5wLCI2i+F3lpdJlPTk5gbWIC0+fsb+OJ9mOffwNbhfSq5k/Qj0n6SNSOicJ3FOiTpSuD1gJuAeUk3KiOBd4koR98Xq3PZ5tvFgNlIAz6fRTBc0lDgZxGxQZ7xWX2SNJBUyvwbYJeIGDfdYzsBB0TE2hL9gYWBgaQR348j+DKPmK0GSd8B5uHbPj7vehWFzSxJvyDNLG0VEc8VPgYBiwBzkgYcvwLerecBnrpPpgqRdBrQPSKOyjuWZiZpHuAF4MCIaJjlTEqje58A34+Ij/KOx+qLpJ7AO8AOEfFs3vFY/VC6Ab6ftIT94IiYMt1jAl4CfhMR9+UUopk1oez6czKwIzA4It7JOaSqqqs9U11wE7BT1hzWcpCWm3AzcG0jJVKZwcArTqSsFJFGf88hLXUw6xRJK5Aaf18L/HL6RCrzE9JyrfurHZuZNa9s2fHVpO03azRbIgUNmkxFxKuk9Zdr5B1LEzuVtAa7EcvuDgFuzDsIq2tXAmtKKrX/izURSRsAjwBHRsQ5UXhJydHAGUUeMzMrO0mzAvcBswPrR8QXOYeUi4ZMpjI3kW56rcokbUua6h1SYPS0rkkaQBoBdgNoK1lEtAAXAUfkHYvVNklDSJ9nO0TEzUWOWZ20P+qWasZmZs1L0gKkojZvA9tkn2tNqSH3TAFIWhR4FlgwvKmyaiQtQ/rj2jQiGq7sblY8YIeI2CLvWKy+SZqDtHdqxXDlR2sl24NwBKkS6qYR8Vo7x94FPBwRF1cpPDNrYpKWJ81IXQac2ewz4g07MxUR/wX+D9go71iaRTZrcwdwTCMmUpmd8RI/K4OI+AoYBhySdyxWWyR1J/Xk24W0B6G9RGo54Mek95KZWUVJWhf4G3BsRHhpMQ08MwUg6SDgBxGxW96xNLpsFPUW4JuI2DvveCpB0tyk6ewFG6JfluUuq872CrDktF5B1twk9SMN2MwCbBsR33Rw/NXAWxFxWhXCM7MmJmlH0kDPzhHxSN7x1IqGnZnK3ApsmX04WWUdQuqpc1DegVTQdsD9TqSsXLKKkHeSlnJZk5M0F6nQxGjS0r6OEqmFgC2AS6oQnpk1KSWHkSrRbuREakYNnUxFxGfAc8BmecfSyLIp3yNJo6jj846ngqZtBDcrp7OBAyX1zzsQy4+kxUmlz/8G7Bada/59KHBlRIysZGxm1ryyZcfnAz8jLTt+Jd+Iak9DJ1MZV/WrIEkLkv6Nh0bE+3nHUynZCPCywEN5x2KNJSL+AzwB7Jl3LJYPST8gvQfOjYjfdGYPQjaLtRvpJsfMrOwk9SVt4fg+sJaLJRXWDMnUHcAGkmbLO5BGkzVquw24KCIezjueCtsJuL2To8VmXXUGcJiknnkHYtUlaXNSVaxfRMRlXXjqgcCfImJ4ZSIzs2YmaU7gYWAiMMgz4MU1fDKVrTn/K7BN3rE0oN8BX5BuBBudl/hZxUTE86TiJjvnHYtVj6R9gD8Cm0fEvV143izA/qQlomZmZZW1F3oKeBLYNSIm5BxSTWv4ZCrjpX5lJmlXYBCwe0RMzTueSspKD89FWoZjVilnAEdLapbrctPKNnOfTOojtXZEPNfFU+wNPBYRb5c/OjNrZpJWJSVRF0bE0Y1+j1cOzfKhfR+wqqT58w6kEUhaETiP1PF6ZM7hVMPOwC0RMSXvQKyhPQK0kKqzWYPKlkdfDWxC2sz9TgnPPxQ4s/zRmVkzkzQYeBA40E3AO68pkqmIGAfcA+yQdyz1TtLswO3ALyPi1bzjqbSsf5Yb9VrFZUUHzgCOyd531mAkzQr8GZgdWD8ivijhNLsAb0TEi2UNzsyamqQ9Sc2/t4yIO/OOp540RTKVuRHvR5gp2fKja4H7IqJZ9g/9AJgK/DPvQKwp3AnMAayTdyBWXpIWAB4H3iHN6reUcI5uwFE0xz5VM6uCbNnxCcCvgXUi4pmcQ6o7zZRMPQIsmvXysNL8hjSienjegVTRzsCNnSlVbDazsqWkZwFH5x2LlY+k5Uk9pG4GDoiIySWeaitgFPBouWIzs+aVVZC9ktSPdY2IeCvnkOpS0yRT2YfXbaQS19ZFkgYB+wLbR8SkvOOphqxR3Y64ip9V13XA9yWtlHcgNvOypuZ/A46NiDNKHZjJln4eDZR8DjOzaSQNAO4F5gHWi4jPcg6pbjVNMpW5CdjF+xG6JiuReQ2wU0R8knc8VbQu8ElEvJl3INY8shK055GWc1kdk7QDaRBvSERcP5OnWw8YCNw1k+cxsyaXFWT7O/ABsHVEjM05pLrWbMnUM0A/Uidn64Ss+/WfgNMjotlKgw/BhScsH5cDG3tZcn3K9iAcSurFt1FEPFKG0x4NnOkyxWY2MyQtS1p2fAepWXipy44to2ZbLSDpdNLP7T0JHchm8K4E+pJGVpvmzSKpNzAcWDEiPso7Hms+kk4B5oyI/fKOxTovWx78O2AjYHBEfFiGc64C3A0sHhETZ/Z8ZtacJK1Fqsh8VERcnXM4DaPZZqYgLfXb2Y0xO2Vv4IfA3s2USGUGAa86kbIcXQDsKGm+vAOxzslm8m8BVgTWKkcilTkKONeJlJmVStK2pNmooU6kyqsZE4pXgdHA6nkHUssk/RA4hVTCd0ze8eRgCC48YTmKiM9Jy0x/lXcs1jFJcwIPA5OAQeVqaC5pSWAD4I/lOJ+ZNR9JvwJ+D2wSEX/JO55G03TL/AAk/QZYICIOyDuWWiRpbuAF4FcRcVfO4VRdVuHmI2CxiBiRdzzWvLLiLy+Q3ovf5B2PFZb9nh4gLcU7ppz7miT9Afg0Io4v1znNrDlkq7DOBgaTlh2/n3NIDalZk6nFgH/Al1vAnMuRKiSNBd4H/hZB027Gk9QDeBB4LiJ+nXc8lSYxP7A+qVFqAF/CCrPBvzePiC1yDc4MkHQ9acnpmXnHYm1JWhW4h1Sk56Iyn3t+4DVgqYj4spznNrPGJqkPqRLzfKSKfV/nHFLD6pF3ANUm0RdiXXhvAAz8OzCZ9O8wJfuaKHEBcHkEn+YZa05OJiUVx+UdSKVICFiH1Hx4Y2Ai0JP0c0+GF/vD/z0rsVoEL+QYqhnAmcBfJP0+IsbnHYx9S9Jg0s3KLyLizgq8xMHAdU6kzKwrJM1BaqPwCfATf3ZUVlPNTEksDTwGzJJ9FTOOdGM9JIK7qxBaTZD0U+B8YLWI+CLncCpCog+p78v6pDL5RXqOxRTQBNKelX0jmFKtGM1ak/Rn4N6I+EPesVgiaU/gVOCnEfFMBc4/G/B/wCpemmNmnSVpYdKy4weAI9xOofKaJpmSWBb4BymJ6mzhjXHAzyK4pWKB1QhJSwNPAJtFxPN5x1MJEr1IyfRKpHLvndEC/AXYNgJfkCwXWTnba4Cl3RMkX1nLiOOBoaQ9CG9V6HWOAZaJiN0rcX4zazySVgL+DJwTEefnG03zaIpqfhIDSTfRA+jaz9wXuEriB5WIq1ZImoVULvM3jZpIZa4klSzubCIFafZqY1JlQ7NcRMSTpOUa2+YdSzOT1JN0HdkMWKOCiVRfUhXHsypxfjNrPJI2IQ3+HuxEqrqaIpkC9iDNSBVc0rXjjvD66zBmDLzzDqy11gwP9yXtI2pI2SjrFaRZuytyDqdiJBYGtiMlRzNYZhl45BEYORLefhu23rrN0/sDB0vMWuk4zdpxOnB09jdrVZZV+bwXmBdYPyI+q+DL7QE8GxGvVfA1zKxBSNoduJbUzuZPecfTbBo+mcqKDRxBgZtogI02gjPPhJ/9DAYMgHXWgXffnfEUwLoS36l8tLn4FbAkcGCDN+bdnwLJdPfucPfd8Oc/wxxzwD77wPXXw5JLtnn+VGDXKsRpVsz9pGI5mwBI9JdYWGIJibmya51VQFZV7+/AB8BWZeu9J/VDWghpSaS5kZRVVD0COKMsr2FmDStdMnQscAKwXraKwaqs4fdMSaxPKltbsODEU0/BlVfCVVe1e5rxwLkR/Kb8EeZH0jrArcCPI+K9nMOpGImewJfQdmZp+eXhH/9IifQ0Dz0Ezz4Lv/1tm1O9H8EilYvUrH1Sr11h68Ph1i+AdUmVKINUjfJzUj+RayNwT6oykbQMaSP3lcCpMz3olGYW1yclTBsy4+9wxL3wyJ6w+OcRa87U65hZQ8sGXi4BVgM2jYhmrEBdExp+Zor0Jutd6IFu3WC11WDuudPyrg8/hAsvhD592hzah3Tj0jAkLQDcDOzeyIlU5rtA984eLMH3vlfwoe9k1QDNqk7iRzDhTLjq+xAbkm6++5MGinqT3uenA59IHO+ZqpmXFf54DDgxIk4pQyK1CvAeqbnvT2j7O1xgI9jlE1gN6TRSw00zsxlke93vBhYC1nUila9muFDPTvrAamPeeaFXL9huO1h7bVhpJVh5ZTj22ILnma1yIVaXpF6k8uCXRMRDecdTBQOhcGnzN9+Ezz+HI46AHj1g441h3XWhX8FFoUykgd4HVj8kBgF/Ay0Asyil/AX1J+3zPAK4xglV6SRtSyrMs1tEXF2GE25Aqpi6EO3s4e0L3bpDL+CXwK1OqMxsepLmBR4FPgO2iIjROYfU9JrhIj0OCpe0Hjcu/e+FF8Knn8KIEXDuubDppkXP0yjOBkYAp+UdSJWMp8iNy+TJqeDEZpul98Bhh8Gtt8JHHxU8T3ca631gdUBiFeB2iuz7LKI/sA3ed1MSSb8Cfk9qdvmXMpxwBdIocld/h4OyOMzMkLQU8DRwH7BnREzKOSSjOZKpjyhyAzxyZFraN/3CjSKLOIK0NKPuSRpCKuu7WxM1cvuUNNJb0KuvwnrrwVxzwaBBsNhi8NxzBQ8NwCNAVm1/pMhNeAeVSPsDv5S8z6+zJHWT9DvgF8CaEfFSmU59Cen38a3Ro2f8mjwZLrig9fP6A3uS9m2ZWROTtDqpEM5pEXFCgxcNqyvNkEzdSTv7ZYYNg4MOSvumZpsNDjkkVXZrZQxwWeVCrA5J3yeNcm4TESNzDqdqIvgaeJKUDLWxwgrQuzf07ZtmpuafH66+us1hk4Eb3bjXqklieWDZQo91ohIppGv8gRUOsyFI6gPcBPwAWCsi3i/TiRcHVqX17PiAAd9+zTdfWipx222FztCDtOTPzJqUpK1Js9s/j4grcw7HWmn4ZCqCkaT9QQX3zJx8Mjz/PLz1FrzxBrz0Epx6apvDRgF/q2igFSZpNtJSoUMi4pWcw8nD2RBjCz0wdCh88knaO7Xhhmnf1MSJbQ6bCJxf4RjNWjuYIns+TzwRTjopVZ6MgOHD01crvYB9pMJFeCyRNDvwECnh2SQivirj6Q8iJUTFbbttugA98UShR3sCuyP1L/SgmTU2SQcAFwODI+KBvOOxthq+NDqAxErAU3RtvXqmJeCKp+BXg8vWW6TKlDYw3wW8HxEH5RxOLqR+P4J3n4R5epQwhjAFeDmCVSsQmllREl8Ac7X+frduaSLjt7+FvfZKFUjvuisVUhk/vs1pRgGbR1DwTr3ZSVqYVPr8AeCIsi9/lj4CFmz3mEcegccfTxlyYd8AO1CO/VtmVheye7fTga1IidR/cw7Jimj4mSmACP4FnAwUnJlox3jo8Tgc9l/gX5LWLntw1XEMMCdwWN6BVJuk3pJOhXH3wAm/gW5dfQ8E6WZ02wqEZ9aRAYW+2cVKpEH6+7dWJK1EGmi7PCIOq9A+0oHtPrrQQqmE6DXXtHdUN2COcgZlZrVLUm/gemAt0v5NJ1I1rCmSqcyZpGVaLZ08vgV4BnptFjFpN+BQ4BZJ50nqW6EYy07SJsABwPYR0XbxWgOTtDLwPPA9YMWIy84CBpOKSHTmpmkK8BWwXkRjFCCxxtC1SqQjB8LWd0oKf834BbxEmjU6r1KvMbpIw/j/GToUnnwS3nuvrO8RM6tP2baMB0i95zaKiBH5RmQdaZpkKoKI4FhgN+Bt0ixVoRvq0aQb6NOATSLSbFZE3AOsAMxHmqX6cVUCnwmSFgGuBXaOiLa7KRqUpJ6SjiftgTgb2HpaQ7tsqdMPgAeBCaSy6a21ZN+/HVgpgmbcY2a1oWD1yC5UIgVmGwV3rRMR8lf6AvYAPgcq/u8yAD5u9ze8224dzUpB+qwq5z4uM6tBkr5L6kf3KrBDRLgdSx1oij1TrWWNLH9Amm1albSUZhyp/PnvgT9HMLn487UdcBFwNXBCRBS6Ic9VVpnqSeCGiDgv73iqRdL3gGtIN0p7RUTRGxmJBUglkLcjNeOdCnwNXAdcFYFHgyxXEn8k3fi3KWBw4okweHDqkTZpEtxzDzz2WNpH1cpoYO4IJlQ63lonScBvgD2BTSPijSq86PnA/hQqJLL66vDww6ma35h2t+S2APMQhYvomFn9yyou3wecB5zn0uf1oymTqXKQNA+pXPrSpJ5NL+Yc0gwkXUFKEndqhj9IST2AI0gJ8tHAVc3wc1tjy0qjPw+0WVrcowf8/vcwZEgqOnHrrXDkkTBhxpRpInBhBIdXJ+LalV0jLgFWAzaLiE+q9MKLA/8G+rR57LLLoF+/NDtV3CTgCiL2r0yAZpY3SRuSWjMcFBG35B2PdY2TqZmQjXLuTBpF+ANwSi3sS5K0Fymp+GG9ViDsCqWGlleT+oHtGeXqD2NWAyReBFYp8enjgWWbfc+fUlnxW0gzfNtHRHWbb0tPAGvSutdU54wDVibizfIGZWa1QNKuwO9Iy/r+nnc81nVNs2eqEiK5EViZdLPzrKQV84xJ0g9I+71+2uiJlKTukg4lLWe8ltQfxomUNZq96XzhnOmNBS5wIqV5gcdIS3+3qHoilexH16vJkj3nSidSZo1HyTHAqcAGTqTql5OpMsiKO2wBXAD8VdKx2ZKSqpI0F6lB8S+iwT98JS1BukHaGvhRRFwSlSlrbJarCP5JKs3flYRqLHAHaclr05K0FPA0aR/CnhExKZdAIv4NbEnXEqqxpIpev6pITGaWG0ndSY14dwRWj4jXcg7JZoKTqTLJZqmGkWao1gaekbRctV4/+8O8Cbg5Iu6s1utWm6Rukg4E/kGqtrdeRPxfzmGZVVQEDwIbkCrDjSH1jipkLGlZ2FnA7hFFj2t4klYHHgdOi4gTct9DGfEosA7wAe3/DqdVE70A2BEPEpk1FEn9SINdS5IqijZNteVG5T1TFZDtpdqHNHV7FvC7iJhS4dc8Ffgx8JOIKFqJsJ5lpd6vIm3G36PRZ9/MWpPoBmwIHAmsSypOMBXoBXxKagVwbQSjcguyBkjaGvgjqTjQAzmHM6P0+bAeqWDORnz7O+wJjCD9Dq8mYmROEZpZhUiaG7gXeBPYuxb22dvMczJVQZIWJd389ybd/L9VodfZCrgQWC0iPq/Ea+QpS073JiWnZ1OF5NSs1kn0B+YkJVJfA18180zUNJIOAH4NbFlrVVbbSCPUc5I+I0YCI/CHsllDyrYnPADcDPw299lyKxsnUxUmqRtwAHA8cDJwYTn39mR7Ap4ENo+I58p13loh6TvAFcBcwO5eV2xmhWTX2tOBrYDBEfHfnEMyMwNA0o+Au0i9Sf+QczhWZk6mqkTSkqTy3ZOBn0XEux09gbS+fnlgVtJeiPeAB8k2UWflfp8lJWg1+ccp8R1gE2CO7FsjgIciaHeNcDYbtTtpmeQFwJm5bR43s+pKf/+rkZqqDyTtIfoYuI+IcQUO7w0MAxYmzUi54baZVYY0J7ApMDfQnbQy4DEi3il8uLYEriTd+/25anFa1TiZqqKsSMTBwDHAccAf2sxSSQOBPYDDgdlIf6g9SUnYRNLa+ovHwaX90pK3CcDPa2m6WEKkfR1HkBLCyaRlLJDi7Qk8ApwDPNZ6aZKk+YHLgYVIex5erlLoZpantOxtJ+AoYEFSkaSewBTS9U+kpdMXTrtxkTQbaTP3SGCXKJBsmZnNNOmHwGGkypyTSPc13Uj3Nd2BF0kDwPeRbUWQtC9pZdJWjbh6yBInUzmQtCxwDTCKlAh9kD2wCvAw6Q+0fzunGD8Juv0cProevldLNw8SfUjl2dcj/QzFmlQGqWrVQ8CQCCa0aoJ8OXCyN2eaNQlpcVK7g9mAWdo5ciIpuTpUqeT5/cCjwCHeS2lmZZcGwi8CdgP60H4l7DHAqyNh8OypUND2pGXHrjrcwJxM5STrQ3UkcAhw9CR4uUe6IWgvAZnBVBjfDfYg4pYKhtppEr1IP8PKpIp7ndECPA8rDYGXLwSWIRXreL5CYZpZrZEWA14gLenrVMuOKTD+tzD+NDgFOLeWZufNrEGkQd6bgM1pf5D7fwImDIcJ34P/jEz72b+oZIiWPydTOZO0wvxww5uwzIC0nKWrWoC1ifhnuWPrKolhwA5Av649c/IEuHYK7HkRcHxEjK9AeGZWi6Q+wNvAAnSx9+FkmNgDtiLiwYrEZmbNTfo1qTpopxKpaSbC1G7wWI+IDSsTmNUSN+3NWUS8+j7c3qfQbNTo0TN+TZ4MF1zQ+qi+pPW4uZJYgLTXoU0itfDCcN998NVX8MkncOGF0L379Ef06A179IA4x4mUWdPZgbS0r+3n0eyzwx13wJgx8N57sPPOMzzcI5WFP6MKMZpZs5H6kva4z5hIdeLerBd06wE/Rvp+9QK2vDiZypvUvScc1BN6tHlswIBvv+abD8aNg9tua3MGYBOk+aoRbjv2LfbAJZfA55/D/PPDSivBuuvC/vu3PqrbFGCvSgZoZjXpKIrtkbr4Ypg4EeadF3bZBS69FJZbrvVRS/qGxcwqYHso0Luvc/dmkAZ7Dq5siFYLnEzlb1PSH1z7tt02ZSRPPFHsiH3KGVRXSHQHDiRtzGxj0UXh1lthwgT47DN48EFYfvk2h/UFDpH8njRrGtLKwCIFH+vXL133jjsOxo6Fp56Ce+6BoUNbH9mLtPfUzKycjgQGtHtE+/dmPYCdkNo/h9U937jmb1Xar1yV7L47XHttsUf7AGuVMaaumpciiRTA+efDTjtB376wwAIweHBKqAoYSFruY2bNYRUKjfwCLLVUWj7z9tvffu/llwuNxPQAVq9QfGbWjFLhiWU6PK79ezNIJdSXLFNUVqOcTOVvLjqq3rfQQmlt3DXXtHfU7OUMqosGknpJFfT44+n+Z9Qo+PhjeOEFuOuugodOys5lZs1hIMUK78wyS7poTO+bb9LSmrY88mtm5VR0gPh/OndvFvi+puE5mcrf6A6PGDoUnnwybcAu4u+wmqTI4wsWeR3GFLyZkdIs1B13QP/+MOecaU/5mWcW/DG6ATXTM8vMKm4cqWdUW2PGwKyzzvi9WWdNG77bcuEaMyunCXR0j9yJe7OM72sanJOp/H1AKm9e3G67dTTyMXVdGBYRyuML3usPs0woFNgcc6RqfhddlPaRf/UVDBsGm25a8OfoBnzVuX82M2sAH5JmpNt66y3o0QOWWOLb7624Irz2WqGjP6hEcGbWpCKmAl+3e0zH92YAvYGPyxSV1SgnU/m7lfZ+D6uvDgsuWKxSzDTjgD+UOa5Oi6AFuB+Y2vqxESPg3Xdhv/1SOfSBA9MS41deaXOaKcDtEUysfMRmViP+UvSRlpY0pX3SSakYxRprwFZbwXXXtT5yNHBRJYM0s6b0R9IMVVuduzcDeIWID8sdmNUWJ1N5ixgB3EOBRARImce0PivFDQeeK39wXXIORaayt9kGBg2CL76Ad96BSZPgkLa1tyYA51Y4RjOrJRETgUspdsOy//6pcs3nn8NNN6VRmddfb33UFNI11MysnC6hWIGczt2bjQbOqkBcVmMUUfh9YlUkrQb8nQINbzthLHAQEcPKG1TXSAh4g1S1pqtJ+hTg9QjcK8as2UgLAf8htUfoqnHA2UTk3rjczBqQ9CCwAcUK5bTvS2ABIgovZbaG4ZmpWhDxAnAaKTHqinHAg8DV5Q6pqyIIYAug3WGaQk8ljd5sVfagzKz2RXwA/IKO9o62NQF4CTi17DGZmSVDSUlR4UI5xbUAmzqRag5OpmrHacD5dP6GYiwpkRpCjUwvRvA2sB5p02ZnLjyTSQUn1o3gvxUMzcxqWcR1wKF0vupVC/AiMDhbKmhmVn4RXwBrkLZTdOZaE6RB5c2JeL6SoVntcDJVKyKCiGOBXYDXSDcLrROSabM4w4GjgO1q7UYigpeAlUiFNcZTODkcmz12E7BiBG3LUZhZc4n4A7Ap8CwpqSo0ojsaGAGcDqxHxKgCx5iZlU/Ee6T7mqtI9y+FVuCMz77uBX5IxKPVCs/y5z1TtUpaBfgVsAqpIWUL8A5wAfBIrcxGtUdidmAPYAgwBykZHAHcCFwTwcjcgjOz2iUtTbr+rUlqeDmeVP78YuA+Ioo2CTczqxipP7AT8HNgHqA7MBK4G/gDEZ/mF5zlxcmUmZmZmZlZCbzMz8zMzMzMrAROpszMzMzMzErgZMrMzMzMzKwETqbMzMzMzMxK4GTKzMzMzMysBE6mzMzMzMzMSuBkyszMzMzMrAROpszMzMzMzErgZMrMzMzMzKwETqbMzMzMzMxK4GTKzMzMzMysBE6mzMzMzMzMSuBkyszMzMzMrAROpszMzMzMzErgZMrMzMzMzKwETqbMzMzMzMxK4GTKzMzMzMysBE6mzMzMzMzMSuBkyszMzMzMrAROpszMzMzMzErgZMrMzMzMzKwETqbMzMzMzMxK4GTKzMzMzMysBE6mzMzMzMzMSuBkyszMzMzMrAROpszMzMzMzErgZMrMzMzMzKwETqbMzMzMzMxK4GTKzMzMzMysBE6mzMzMzMzMSuBkyszMzMzMrAROpszMzMzMzErgZMrMzMzMzKwETqbMzMzMzMxK4GTKzMzMzMysBP8PRq1rtrcBITYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1080x288 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Gnerate a connected graph with 10 vertices\n",
    "n = 10\n",
    "G = nx.Graph()\n",
    "G.add_nodes_from([0,1,2,3,4,5,6,7,8,9])\n",
    "G.add_edges_from([(0,1),(1,2),(2,3),(3,4),(4,5),(1, 7),(5,6),(6,7),(7,8),(8,9),(9,0)])\n",
    "    \n",
    "k = 9 # Set qubit (vertex) limit\n",
    "S = NaiveLGP(G,k) # Partition G into two subgrahs once\n",
    "sep_node = list(set(S[0].nodes).intersection(set(S[1].nodes))) # Obtain seperation nodes of the partition\n",
    "\n",
    "# Show graph illustration\n",
    "options = {\n",
    "    \"with_labels\": True,\n",
    "    \"font_color\": \"white\"\n",
    "}\n",
    "node_color0 = [\"red\" if i in sep_node else \"blue\" for i in range(n)]\n",
    "node_color1 = [\"red\" if list(S[0].nodes)[i] in sep_node else \"blue\" for i in range(len(S[0].nodes))]\n",
    "node_color2 = [\"red\" if list(S[1].nodes)[i] in sep_node else \"blue\" for i in range(len(S[1].nodes))]\n",
    "\n",
    "fig, ax = plt.subplots(1, 3, figsize=(15, 4))\n",
    "for i, a in enumerate(ax):\n",
    "    a.axis('off')\n",
    "    a.margins(0.20)\n",
    "nx.draw_networkx(G, pos=nx.circular_layout(G), ax=ax[0], **options, node_color=node_color0)\n",
    "nx.draw_networkx(S[0], pos=nx.circular_layout(S[0]), ax=ax[1], **options, node_color=node_color1)\n",
    "nx.draw_networkx(S[1], pos=nx.circular_layout(S[1]), ax=ax[2], **options, node_color=node_color2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that in our case, both subgraphs are under the qubit limit and can be directly solved by QAOA. But if the subgraph still exceeds the qubit limit, we will recursively do LGP through recursively doing DC-QAOA."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Graph Reconstruction (GR)\n",
    "\n",
    "Once QAOA is applied to partitioned subgraphs, we should combine the solutions for the parent graph. Note that since two partitioned subgraphs share some set of nodes, we need to make sure that the labeling sets which these nodes are in are consistent. For example, if in one subgraph, the shared nodes are all in $S_0$, while in the other subgraph, some of the shared nodes are in $S_0$ and some are in $S_1$. In this situation, corresponding solutions cannot be combined. So we provide several possible cuts for each partitioned subgraphs to avoid failure for the combination. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-11T06:51:47.629119Z",
     "start_time": "2021-05-11T06:51:47.585155Z"
    }
   },
   "outputs": [],
   "source": [
    "def GR(str_cnt1, str_cnt2):\n",
    "    com_cnt = []\n",
    "    n = len(str_cnt1[0][0])\n",
    "    com_index = []\n",
    "    for i in range(n):\n",
    "        if str_cnt1[0][0][i] != \"x\" and str_cnt2[0][0][i] != \"x\":\n",
    "            com_index.append(i)\n",
    "    \n",
    "    for (str1, cnt1) in str_cnt1:\n",
    "        for (str2, cnt2)  in str_cnt2:\n",
    "            # Check equality for each bit in common nodes\n",
    "            validity = [str1[i] == str2[i] for i in com_index]\n",
    "            if False not in validity:\n",
    "                com_str = [[0]] * n\n",
    "                for i in range(n):\n",
    "                    if str1[i] != \"x\":\n",
    "                        com_str[i] = str(str1[i])\n",
    "                    else:\n",
    "                        com_str[i] = str(str2[i])\n",
    "                com_cnt.append((\"\".join(com_str), min(cnt1, cnt2)))\n",
    "\n",
    "    # Sort string-count map by counts in reverse order\n",
    "    com_cnt.sort(key=lambda tup:tup[1])\n",
    "    return com_cnt[::-1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We still take the above example for illustration. For two subgraphs partitioned above, we apply QAOA to find their approximate max cuts and then use GR policy to form the max cut for the original graph. The code presented below has achieved this. Note that there is a parameter $t$ controlling how many possible cuts are kept, and if this $t$ is larger than $2^k$, all cuts are provided and the combination can definitely happen. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-11T06:52:48.693469Z",
     "start_time": "2021-05-11T06:51:47.651060Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Max cut for the first partitioned subgraph: \n",
      "{'01010101xx': 0.07887396513978905, '10010101xx': 0.07887396513978903, '01101010xx': 0.078873965139789, '10101010xx': 0.078873965139789, '10101101xx': 0.040906327139884305, '01010010xx': 0.0409063271398843, '01010110xx': 0.04004434869249364, '10100101xx': 0.04004434869249363, '01011010xx': 0.040044348692493625, '10101001xx': 0.040044348692493625}\n",
      "Max cut for the second partitioned subgraph: \n",
      "{'0xxxxxx101': 0.4556003303152275, '1xxxxxx010': 0.45560033031522745, '1xxxxxx100': 0.0254542520553071, '0xxxxxx011': 0.025454252055307092, '0xxxxxx110': 0.010740876742244158, '1xxxxxx001': 0.010740876742244157, '1xxxxxx000': 0.0022930869455346486, '0xxxxxx111': 0.002293086945534646, '0xxxxxx100': 0.002293086945534642, '1xxxxxx011': 0.002293086945534638}\n",
      "Combined max cut for the original graph: \n",
      "{'0101010101': 0.07887396513978905, '1010101010': 0.078873965139789, '1010100100': 0.0254542520553071, '1010010100': 0.0254542520553071, '1010110100': 0.0254542520553071, '1001010100': 0.0254542520553071, '0101101011': 0.025454252055307092, '0101011011': 0.025454252055307092, '0101001011': 0.025454252055307092, '0110101011': 0.025454252055307092}\n"
     ]
    }
   ],
   "source": [
    "# We use direct QAOA to compute approximate max cuts for two subgraphs\n",
    "# In the next section, we will compute these using DC-QAOA because subgraphs might also exceed the qubit limit\n",
    "import paddle\n",
    "from paddle_quantum.QAOA.maxcut import find_cut\n",
    "\n",
    "# Set QAOA parameters\n",
    "p = 3 # Number of layers of QAOA circuit\n",
    "ITR = 100 # Iterations of the training network of QAOA\n",
    "LR = 0.5 # Learning rate in the training network of QAOA\n",
    "# Set graph reconstruction parameter\n",
    "t = 10 # Number of partition strings kept after graph reconstruction \n",
    "paddle.seed(999)  # Fix the seed\n",
    "\n",
    "# Start graph reconstruction procedure\n",
    "S_str_cnt = []\n",
    "for si in S:\n",
    "    siv = list(si.nodes)\n",
    "    # Compute the subgraph's maxcut\n",
    "    tmp, si_str_cnt_relabeled = find_cut(si, p, ITR, LR)\n",
    "    # Make the subgraph's maxcut match the original graph by relabeling \n",
    "    si_str_cnt = []\n",
    "    for str_relabeled in si_str_cnt_relabeled:\n",
    "        strr = \"\"\n",
    "        for i in range(len(G.nodes)):\n",
    "            if i in siv:\n",
    "                strr += str_relabeled[siv.index(i)]\n",
    "            else:\n",
    "                strr += \"x\"\n",
    "        si_str_cnt.append((strr, si_str_cnt_relabeled[str_relabeled]))\n",
    "    si_str_cnt.sort(key=lambda tup:tup[1])\n",
    "    S_str_cnt.append(si_str_cnt[::-1][:t])\n",
    "\n",
    "# Once we have already obatined max cut strings for two paritions, we perform the graph reconstruction (GR)\n",
    "print(\"Max cut for the first partitioned subgraph: \\n\" + str(dict(S_str_cnt[0])))\n",
    "print(\"Max cut for the second partitioned subgraph: \\n\" + str(dict(S_str_cnt[1])))\n",
    "out_cnt = GR(S_str_cnt[0], S_str_cnt[1])\n",
    "print(\"Combined max cut for the original graph: \\n\" + str(dict(out_cnt[:t])))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The top several possible max cuts for the first subgraph include {'01010101xx', '10010101xx', '01101010xx', '10101010xx'} and those for the second subgraph include {'0xxxxxx101,'1xxxxxx010'}, where 'x' indicates those missing nodes in the subgraph. From this example we see that the possible maximal cuts of separated vertices 0 and 7 in the first subgraph, '01010101xx', belong to $S_0$ and $S_1$, respectively. And they also belong to $S_0$ and $S_1$, respectively, in the maximal cut of the second subgraph, '0xxxxxx101'. so we can integrate these two maximal cuts and get the maximal cut of the original graph, ' 0101010101', as shown in the third diagram below.\n",
    "\n",
    "Graph illustration is shown below. The left and middle subgraphs are subgraphs with approximate max cuts, where red and blue nodes represent $S_0$ and $S_1$ and dashed lines represent cuts."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-11T06:52:49.021611Z",
     "start_time": "2021-05-11T06:52:48.704511Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1MAAADnCAYAAAD7CwxiAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAABKA0lEQVR4nO3dd5xcZfXH8c832XRC6L03KSIBUaQH6VUkSA8gKiX0jkovAoIUAVF+0kEQkCqKAtJVuoKA9CZVasqmkvP747kxm92Z3dnJzNwp3/frtS/xzp07ZzOzd+65z/Oco4jAzMzMzMzMeqdP3gGYmZmZmZk1IidTZmZmZmZmZXAyZWZmZmZmVgYnU2ZmZmZmZmVwMmVmZmZmZlYGJ1NmZmZmZmZlcDJlZmZmZmZWBidTZmZmZmZmZXAyZWZmZmZmVgYnU2ZmZmZmZmVwMmVmZmZmZlYGJ1NmZmZmZmZlcDJlZmZmZmZWBidTZmZmZmZmZXAyZWZmZmZmVgYnU2ZmZmZmZmVwMmVmZmZmZlYGJ1NmZmZmZmZlcDJlZmZmZmZWBidTZmZmZmZmZXAyZWZmZmZmVgYnU2ZmZmZmZmVwMmVmZmZmZlYGJ1NmZmZmZmZlcDJlZmZmZmZWBidTZmZmZmZmZXAyZWZmZmZmVgYnU2ZmZmZmZmVwMmVmZmZmZlYGJ1NmZmZmZmZlcDJlZmZmZmZWBidTZmZmZmZmZXAyZWZmZmZmVgYnU2ZmZmZmZmVoyzsAy4/ESsAewFLAYOBj4CHgNxGMyzM2M2tdEkOBXYB1gLmBduA14PIIXsgzNjOrXxKLA3sBKwBDgU+BJ4ArI/g4z9hKJg0Etgc2AuYFJgFvAdcQ8USeoVlhioi8Y7AakhCwHfBDYEWgHzMn1eNII5ZXAz+N4LWaB2lmLUliaeBoYDdgGjCkw8NTgKnAv4DTgVsj8BeYmSGxIfAjYC3SNUz/Dg+3Z9tuA34SwTO1j7AE0kLA4cAPsi1DOzw6DZgIvAmcSUqsvqhtgFaMk6kWItEGXEZKpob0sPsU0t2QbSK4r9qxmVlrk9gIuAUYSM+zJsYDNwB7RzC12rGZWX3KbhCfDBxGmmHTnS9I1zXfj+C6asfWK9JqwD2ka7P+Pew9njSLaDsiJlQ7NOuZk6kWkZ1wrgO2pucTTkftwKYRPFyVwMys5UmsB/yR3p+bbgFGeYTKrDVJnAYcTM83iDtqB/aM4MbqRNVL0krA34DZAJX4rAnAY8BGRPiGUs6cTLUIidHAT+ndCWe6McASEXxa2ajMrNVJzAW8wcxTWko1Hjg8gl9VNCgzq3sSWwA30rubMNO1A8MjeLmyUfWS1J80dW9+Sk+kpmsHLiDimIrHZb3ian4tQKIPcCwFEqmrr4Z334XPP4cXX4Tvfa/gIdpIhSrMzCrte0DfQg8svzzcey989hm8/DJsu22XXYYAx2Uj72bWWk6gQCJV4nVNP+CQ6oZXkm1J57GZz2H9+8Ovfw1vvAFjxsDTT8Nmm3V+7mBgf6RBtQjUinMy1Ro2Ig0fd3H66bDEEjBsGGyzDZx6Kqy2WpfdBgNH+oLFzCopu9FTcK1D375w223w+9/DXHPB3nvDNdfAsst2OcwwYIPqR2tm9UJieWDlQo+VeF3TD9hTKmu2TiUdQ6FR+bY2ePttWH/99IsceyzccAMsvnihY+xQ7SCte06mWsMhFEmmnn8eJk9O/x2RfpZeuuAxhgLrVic8M2tR61Nk6vHyy8NCC8G558K0aXDfffDIIzBqVJddh1Afd5jNrHb2pUihml5c10wDRlYpvp5JywHLF3ysvR1OOgnefDP9AnfeCa+/Dl/9auc9ZyPdkLIcOZlqDcvTzVzciy6C8ePTcPh778Ef/lBwtz7AMlWKz8xa07IUmeJXiARf/nLXzaSeMmbWOr5CGl0qqMTrmtnI97pmGWBySXvONx8stxw891yhR5eoYExWBidTraHbxZn77w9Dh8I668DNN8OkSQV3awNmr0ZwZtayhlLk7vKLL8KHH8KRR6YZLxtvnGa8DC58Nst7qo6Z1Va3BWtKvK4Bnl1NUnT42Rug07Y7sm13dNyebdu7075bS1qo07ZLsn2f7LDtXWC2iakVRPfa2uDaa+HKK9OJsauej2FV5WSqNYzvaYdp09IUmkUWgf32K7jLVODzSgdmZi1tDBTuEzV1aio4seWW8P77cPjhacnAf/5T8DjjqhijmdWfMT3tUMJ1DbDykxGhDj+XAHTatnW2beuO27Ntl3Ta946IeLfTtr2zfb/aYdtCwNiBqRFvcVKqqDF5MhxwQLG93GsqZ06mWkI8nybd9qytrejc4gBeqmRUZtbyXiKtWyjo2WdhxAiYZ55UyGqppeCxx7rsFsDz1QvRzOrQPyhxilw31zXjyPe65iV6atB76aUw//wwcmS6w1TYa5UOzHrHyVQTk7SopLnhhtsKDU7NOy/suCMMGQJ9+sAmm8DOO6dSxF19PgX6HC7p25IGVD14M2sFD9LNiPfKK8OAATBoUBqZWnBBuOKKLruNB86pXohmVocupsCNmN5d1wBwczWD7FbEq8CzRR+/+GJYYQXYemuYWHQAayzwsypEZ73gZKrJSBoiaU9J95Lu3HwNXrgUhrzbed+INPT9n//Ap5/C2WfDIYfAHXd0Oex4mPoTiDtJncbfljRU0iBJLpduZmWJIICzSc0nuxg1Ki0e//BD2HDDtG5qctd70R8DD1U3UjOrJxG8AvFUge2lXtdMBi6NyH2K3E9JCdHMFlsM9t0Xhg9P85zHjk0/u+wy027ZlKObqh+mdUdR2uwvq2OS+pJ6Sb0KTAIuAq4Cfh8RE9M+/AA4l/IWan8GLB6R5ihLWiAi3pd0HDAKuBq4JiJen9Xfxcxai8QcEG+AhvX+2V9MhL6jI7i80nGZWf1KN3L/uiOseg0MLLkiaAftwMoROU+Rk9qAN4AF6eUAxySYej6MPxq2j4h7qhGelcYjUw1M0jySzgbeBk4BFoyItyNim4i4aXoilfk1cCNF7gB3ox3YdHoiBRAR72f/eSopmZofeCAbFVtUKueiyMxaUQSfwT2jYHLRBQGFTZ0E1wBtn1QjLjOrP5L6SDoEuDRizeth4ImUUGSrkwnALrknUgARU4ENSeu3ejO6MWEA3HsS7Az8UtKCVYnPSuKRqQYjaQFgF+BR0lzbo0mjQi/0/Fz6kkatdqPnEarJpCozW0TwSAlx9YmIadlJ7kTgLtLo2J8i4ouenm9mrUfSsqRz2N4QawN3AgPoaVF2uni6AgZfBRNuBbaNiK6lKcysaUhaGricNBCwZ0S8IiHSOeQ4YBDd9NQkVQ6dDOwWwS3VjrdXpC8D95F6X/VU6nw88GdgZyImSWqLiKmSjgceiIgHqhytdeKRqQYhaUFJfwReAL4MjImIMRHx41ISKYAIvgD2A3YgrTGYSNdqOGOzn4uAL5eSSKVjx7Tsf88DlgLuBw4A+kj6uqRVvb7KzKaT9A1SAYq/R8S0CB4indt+QbpL27nc+fQbPA8A34nggIj2x4BVgcclfUmSv9PMmkyHa4eNgNuA9SPiFUjrLiM4A9iEdDNmIl3LjY8njUZdBXy17hIpgIh/ASuS1lB9Std1VFNJv8NTwF7A9kRMSk+N6aP6/wCuk3SepG77i1pleWSqTmUXBesBu5NGeW4BRgK3R0Rvp+oVeQ2WJo1SLUkaqfov8Dfgxogeeh/06nW0O3AS6YR2FXBuREyp1PHNrLFIWgJ4jHR3+Q9dH2cQ8B3gG8C8pHPHa8DVEXRZm5ldbN0FfATsFdlFhpk1NkmLAZcBZ0bE3T3vz4Kk5QdfAoaRCtQ8Dfym43KFuib1A7YBRpCWUUwE3gKuI+K57p+quYELgKcj4qwqR2oZJ1N1Jkui5iNN4/uMlHxc22GdUkPKfq+1gc2BHwPbAYOBWyLCDTfNWoSkpSPiVUmLREThFrzlHXcwcC3pAmq7iPisUsc2s9rKbpDsBZxBKv19docRGOtBds21FikpO77TGnqrME+JqANZIYkDJD0GHAR8AGwREatExM8aPZGCNA0wIh6KiB9FyuDHAjsC/5F0pVJFGzNrUtnC8TOAW7I5/hVLpACyEfvtSaPrc1by2GZWO1ki0BdYE9ggIs5wItU72dKLl4ClgSclrZ5zSE3NI1OFpD/kTYGjgK+RFjVOJc1jvQz4BbN4IZA1vh0EtAEvA38gjULd20onDUnzkU6Wv5V0FDA3cFX0MJRdwoH7kCrkHEk6IQ8GviC9h1eQ3sM3Z+k1zKwk2fnuMtKU4m0i4qMqv56AXwIXRkTxpphmVjeyv9vdgEOANVrpWqhasn/TnYDRwAgXBKsOJ1OdSbuShpQHA0ML7DF9qPQB4LtEvFf6oSXShf3upDuox0XExZIGRUTejeNyJ2kFYA9gV1K593WAiN5+SKXvAOeR3r/Z6FrdZxKpBOkjpPfw7VkK3My6JWle4HjgqFqd6yTtDJwP7BQRf6nFa5pZeSTND/yKVMBqj4h4OueQmkp2/dkPuJK0/uwf+UbUXJxMdSSdAhxGSqR6Mn2kal0iXuz+sFqKNKXyI1LpyxtI66DemrWAm1PWhHjFiHhW0q+ARUijdrf3eCEm/Yi0JqvU9/BzYERWScfMKkjSIsAJwP4R0blyaC1efwPgt6Q7ss/X+vXNrGfZNP+VSTeZT3YBmerIEqrdgbNIRSrOcDGwynAyNV3qj3QapV2ETxekCnir0Gldk6SBpA/t7sBywGERcU1lgm0dkoYA3yb9Oy5J+recC/h0ejn2DjvvDZxL79/DT4Dhszp108xmkLQyqVTxz4Gf9XqEuXJxLAS8BywDvJJXHGY2M0nzkFohPBcRJ+UdT6vIbnL9H3BaRDycdzzNwMkUgLQoaaFeT43SCpkC3EbEd5TKWW5GGu34O+nDehNwl7P/WSdptogYJ+mXpH/nq4GrI+Il0hSBNyjvPZwK3E3EFpWL1qx1ZVN2ngEOjojr6yCePqRpvU8BB3ndgFm+JG0LXEyqwHmclzrUliRFREg6nLR2/2yfF8vnan7JaLrvmg3LLAMTJsDVV3d+pF/AVktKlwDvkDpxD4mIyRGxR0Tc4USqMjqUUN8P2JbUG+tSSXoBjit4W6B/f/j1r+GNN2DMGHj6adhss857tQEbkO5gm9kskLRERHwAfK0eEin4X2WrzUi9Z37nhpZm+ZDUP/vPpYDvRMQRTqRqr8MI/e9IBdcelvSlHENqaB6ZSn/YH5J6kxT3pz/BoEHw5pswatRMDwVMuAn+skO6C/tq9YK1gqS+42HskFQdcWaDB8ORR8IVV8Bbb8EWW8B118HKK6f3coaJwNlEHFejqM2aSjYf/0hgb+DL9djXJLuQO4M07fCdvOMxazjS0qS/8ZVIRZ4+J434/h89/E1J2oI0GrVeuJpu3chG7vcDvhQRB/W8P0NJhcI2JlVgngi8CVwOPBpByyUWTqakDYGbgdmL7rPjjrDddvD882mEqlMylXmFiGWrFKV1R1oLuIvC1Re7+uc/4aST4OabOz/yHyIWrXB0Zk0vKxpzPrAuqUdeXScqWby/AH7qG2BmJZA2AY4DVifNaurf4dHpN07uB04m4m8zP1VDSRV2NwD2ioj7qxytlSkrmHYxMLrzuVFiSeBHpERqGml20HTTgAmk9alnAFdE0DLTBj3ND+ajuyl+Q4fCySfDYYf1dJy5KxmU9cr8UOKdkPnmg+WWg+cKtrFyo0+z8gwlXVytV++JFEC2NuAfwEOSvpZzOGb1SxLSqcAtpHYlA5k5kSLbNpA0XewepNEdnj6I9P38PrCKE6m69ybwJ+DvkkZno1ZIrAf8E9iTNAtoSKfn9cm2LUO6sfZ7qVfFwBqak6lUd794MnXKKXDppfBOj9cHbZUMynqlX0l7tbXBtdfClVfCiwWr2fs9NOsFSfNkBWEmRcTeEfF53jGVKiIuBvYF/iBp8bzjMatTp5Ca6JZyYaxsv7M+kQ7Kzg03RsS4iPhxRIytYpxWARHxRUScQzbLAJhTYg3gj6SbZqVcJw0B1iclVC1xXeVkKvWKmlbwkVVWgY02gnPPLeU4Pknk5xN6GpmSUvGQyZPhgAOK7TWu2ANmNjOltRN/Jf391d36qFJExO3AahHxpqQl847HrK6kqX2H0nUUoieDB8F5G8MCpClh1mAi4t8RsRXERBj3EL1rOQNp9GoNUp/BpudkCh4FBhR8ZMQIWGKJVLjgvffgiCNg5Eh48snOe04F7qlqlNadJyn2Hk536aUw//zp/Zs6tdAe00gNlc2sB5KGAQ8C50bEjxq5d1NEvJ01Db1D0ilZIQ0zS2ukil9EF69yzACIP8O4RhqttoJ2hkFdKlLvvz88/jhMnAiXX170uYOBg6Qers+agAtQAEi3ANvQObkcNAhm71CX4ogjUnK1337w0Ucd92wH1iLin1WP1QqTrgF2Avp2eezii2H48DTKOH58sSOMBzYk4tGqxWjWBCQtmiUgS0XEa3nHUymS5gPuAF4AfuCWFtbS0sjzv+iud2M3VY4zE4GFiPi0OkFaNUmI1IN1mc6PffvbMG0abLpp+gh897tFDzMW2CeC66oXaf48MpWcTapCMrMJE+CDD2b8jBuX0vCZEylIlfycSOXrHGBSl62LLQb77puSqfffh7Fj088uu3Te813gseqHada4JO1DWpg8ezMlUgAR8SHwTdLNsTnyjcYsd3vT3TXijjvCZ5/Bvfd2d4xpwM6VDctqaDiwYKEHbrkFbrsNPv64x2MMJU0VbWotsTCsBH8F/saMSjWFnXRSoa0TgIOrEpWVLuIppHtIfQ9m9Jt66620Xqp7E4AD8TCtWUHZ1LdTSKO/IyJiTM4hVUVEjAdGS+ov6VfASRHxbt5xmeVgJbpW7UumVzn+5jfh+9/v7hiDgeWqEJvVxlJQkfLmTb8e1SNTQHYRvS3wIr1bSD0B2BeX+qwXO5KmJfSmm3o7cCgRf6pOSGZNYTCwMLBWRLycdzA1MAV4HfirpBXzDsYsB8X7NpZe5RhgropFZLU2G5XJEwb1vEtjczI1XbojuRZw3xcwsWCJghnGZz87EXFV9YOzkkRMBNYjNfBtJxUGKWZ8ts+eRPyqBtGZNRxJs0u6EOgbEd/NpsI1vUjOAI4F/pKtpzJrJYULR/SuyjFAzxPBrF6NpVi1695pr8Ax6pqTqY4i2onY4nbY9t3UtGwi6YQy/Wc88AZwJLAgqayu1ZOIiURsB3wDuIo0StXxPRwHvH0p3LtaSoZvzC9Ys/olaSFSxb6+tMCXYSERcQ3wjYj4UNJiecdjVkNPUWimTulVjiF93z5b1Sitmv5NZZYDvVCBY9Q1V/PrJOs10i8iXiKV//0KaTHyJFIH72e9tqaBSLOT3sM5gcnAB8A/BaOA3SNiozzDM6tHkgaRpsz+H3BmI5c+rwRJA4HngQtJ5eBb+t/DWoC0MPAKndeRl17lGNJNmPmymT/WgCT+AazSeXvfvtDWBiecAIssAj/4Qeo680XXFVZjgd0iaOrBBydTnUi6CnghIk7POxarHkn9SF8U20fE43nHY1YvJC0UEe9KWrZF1keVRNKiwB9JPQUPj4hKLMw2q1/SH4FNgeJVnE44IfWb6loafQpwKRH7VS9AqzaJnYFf0WkN3QknwIknzrzviScWrNP2EbBAREUKWdQtJ1MdSFqcNLS9dER8lnM4VmWSfgBMCq97MwNA0g7Az4GVI+K/ecdTbyTNQWqlcYS/I6zpSWuSbh4Ub9xbXDswHN+QaWgS/UnFeBag90uDxgPHRVDyArtG5WSqA0nnApMj4ui8Y7HakdQnIiqxyNKsIWWlzw/NfraMiGdyDqmuSRpMSqqOiwgvsLfmJY0GzqJ3CdUEYBcibq1KTFZTEl8CHoNpQ6FPj71mMu3A7cAuETR9ouECFDM7GTgz7yCsdiTtCPwy7zjMctaf1KBxLSdSJZlIWgvwSLbO1qw5RfyCdJOlnZ57Dk3J9nMi1UQieBF2ORs+mwZRSvug8cD1wKhWSKTAydT/SNoJmC8iPsk7Fqupe4DtJS2SdyBmtSZpkKRzgCERsXtEvJ13TI0gIqZlMxguBB5SKlZk1pwiLiG1jrl+EsS0rr0cx5GSqEtJU/turXGEVkWSvgrXHQT7rgM6jbQOauzMe8XkdI+p/W/Ad4DvR3TbnqapeJofIGko8BqwZkS8knc8VlvZxWRExOF5x2JWK5LmBm4D3gb2jIhJOYfUkCQtHRGvSlo4IkrqYmrWqB6V5l0jXSwvR6qS+zGp8udvXbWv+WTfE08AR0XWSkaiL7AFsDZpLVU78BZc9eeI3Z/KLdgcOZkCJB0GrBERO+Ydi9VeNiq1fUScl3csZrWQVbN8ilSd7hivGZw1Si0YngeOj4jL8o7HrBqU1k89EBHP5R2LVZ+kvsCdwLMRcWQJ+w8Cvh8RF1Q9uDrT8slUtvD6MWCfiGjJjNoSSQtGxHt5x2FWTZLmj4gPJK0QEU3fTLFWJC1HSk6vAk52LyprJpL6A+8Cq0XEW3nHY9Un6WRgXWDjiOhxyp6kPsCbwOYR8a9qx1dPWn7NVPaFt5YTqdYmaQngH1mVLrOmJGlz4FlJizqRqqyIeIm0rmRpyislbVbPNiX14HQi1QIkbQ18F9iplEQK0lpSUuGJnasZWz1q6WRKUl9J1wKz5R2L5Ssi3gAeAb6XcyhmVSHpe8DlwLYuNFEdEfFBROwO9JF0viR/t1iz2AT4Td5BWPVJWoZUTGSHiPigl0//DbBDNuurZbR0MgVsCywDfJZvGFYnzgSOyNaTmDWNbO77CGD9iPhrzuG0ggnAIOABSQvkHYxZBRwE/DrvIKy6stk5NwMnRsTfyjjEP4C1W22ac8smU1nWfAxwRqu96VZYRDwKHJh3HGaVIqmfpNOBeSJiVES8mHdMrSCbFrMPcCvwcLYw26whSdoY2CYipuQdi1VPdl18CSkhuricY2TX04MlfaeCodW9lk2mgHlJiylvyzsQqyt3ABtlCynNGlY2xex24CukJopWQ5GcAmwZERM8QmUN7CBg9ryDsKrbH1gZ2HcWBxn6ARdIaqtMWPWvlS8YP4qIb7kksBVwKrBV3kGYlSu7GfAnUg+pb0XEuJxDalkR8aKkeYF/ShqZdzxmvZH1GVqPNMpqTUrSWsDxwHYR0T4rx4qIl0nfPRtUIrZG0JLJlKTVgT/nHYfVn+xuzJnAD1ttAaU1B0nzZjeJ9ia1fGiZLvT1KiL+C2wGnC/p4LzjMeuFDYC7ImJs3oFYdWSj5jcA342IVyt02N8AO1ToWHWvJftMSboReMRNWq2QbLH+i8CoMhdgmuVC0jrATcCIiPh33vHYzCQtDvwE2CsiJuUdj1kpJPWPiMl5x2GVlxXcuge4PyJOqOBxZyPdn26JKeYtNzKVNVYcgavSWBER8QXpbtyjecdiVqpsCtnNwO5OpOpTRLwZEbsCQySdJWlg3jGZFSNpYUk/dCLV1M4A2oGTK3nQbGr5apK+Ucnj1quWS6ZIzRSP8BoC607Wh2cHSV/OOxaznmRTUrcBNo0IT2Guf+3AEsCfJM2ZcyxmxexEah9jTUjSDsC3gV2zm8iVtjxweBWOW3daapqfpGHA5IiYkHcsVv8kHQmsGhG75B2LWSFZoYkfA5dHxH/yjsdKl713ZwMbAqu77LTVG0lPAMdExD15x2KVJWlF4AFgk4h4ukqvMRfwOrBoRIypxmvUi1YbmToWqNicUGt6vwI2kbRU3oGYdSZpAGmR78akkQ5rIBExLSIOA3aOiCmS5ss7JrPpJC1IaiFzX96xWGVJmh24BTiyWokUQER8QkrYtq7Wa9SLlhmZyqZSvAoMj4i38o7HGoOk04BPIuJnecdiNl02re9OYBxpjdTEnEOyWSBpMeAJYDdP07R64cITzSf77vgd8GFE7FuD15sH+LRK0wjrRislUz8Glo2IPfOOxRqHpDaXlrZ6kk2d+JTUXPFf7pXXHCStS6rEeFREXJl3PNa6sgvu04FTvb68uUg6ChgJrFeriqKSdieV1/+wFq+Xh1aa5ncXcEreQVhjiYipkjbLTkBmuZI0HHgGWC0innEi1Twi4iFSpdlNsvYMZnlZDfgO0BJlrVuFpG8ChwLb17g1w6akz1PTaolkStLXgTcq2IzMWsuLwFFZAROzXEjamNRs/NCIeDLveKzyIuKFrHT6PJJOkdSWd0zWknYGrotWmbrUAiQtClxLmkr8do1f/jdAUxfyavpkKmtIdgOwbN6xWGOKiNeBPwH75B2LtbSdgJERcWPegVjVtQNfB26VNCTvYKzlbEa6ALYmkBUrugk4LyLuzSGEu4EvSVokh9euiaZfMyVpV+AHETEi71iscUn6CrBPROyfdyzWOrK1C4cCt0bEa3nHY7WT3Qi8hNTnZ31P6bRaceGJ5iLpYmB+0s24XC76JS0QEe/n8dq10NQjU9mFyBGkDs9mZcvWp+zvtQxWK9kUr18CuwHujddisr5TewH7R8Q0SXPnHZM1P0n74ka9TUPSnsA3gT1znrb5uaTROb5+VbXCyNQSwJue+2uzKvss3UZa/N/UZT4tf5J+A8wFfCcixuYdj+VH0gqkfj/fjoi/5R2PNSdJ/YF3Sd9xbiHT4CStSlpnu35EPJ9zLH2At4BNI+K5PGOphmYfmTqc1CPIiZRVwpukvj4j8w7EmpekObJR9TOBrZ1IWUS8QBqlul3StjmHY81rE+DfTqQaX9ZC43ekke1cEylITcqB60jFTZpO0yZTktYBRpMW8prNsiwpPwM4JrvYNasoScsBT5LuJP4zm+plRkT8AdgcGOnzj1WJC080gWwU6Frgloi4Ie94OrgO2LEZz19NMc1PYgAwGBgTwRdpm+4Afh8Rv8o1OGsq2UnqOOCnETEhbWMQMID0+fMicZshfV5mByaRfV6K76o1gVuAH0fEpbUIzxqTpIWBPYAzui1MkdZ4zg60U9u+MlbnJAQMAvoBYyOYll3ktvkmTmOTdBKpZ91G9fReZp+v+SPifdJ/DyEN6oylwZORhh2ZklhS4hyJz0mjT+8DUyRel945BuZYEXAXeauoiJgWESfBTWtLXCExARgLfED6/D0rsVuW4FsrkgYijUL6FzCF9NkYizQB6TKklYs8c3dgLydSVoKJwJbA1VnZ4xmkYUgHIr1B+vy9D7QjfY70M9LaT2tREqtKXEMqajMG+BCYIn3+Clz5U4imGzVoJZK2Ar4H7FhPiRRAAPfAiLekJ4HJwKfAR8AUpPuQNs9uQDachhuZkpiPNHy5DikZ7F9gt/EQfUC/AI6ePlplNqskloP4LUxcBQZMgz6FqvtNX+NyHPDzCBrrj8zKM6OM+YnZlqEF9ppKusB9AdiRiFck7Q08EBEv1iROawqSBpG+CwcAWwX0Bc4i9cObRrrr29kk0jXNQ8AuRHxUo3AtZxIrAtcDS5M+MwW+uyZPhP6TgWMiuLimAdosk7Q08Ddg24j4a97xzEQaAVw5DeYBBhfJmMYB44HvE/H7WoVWCQ2VTEksAfyV9Gb0K+Ep7aQvja0jqKsM3RqPxNeAe4DZKG1Ut500Orq/E6omlxKp6WXMB5fwjGkBY7eDW2+FNYHNsubQZiXLWjV8+WV4YTH4c3/4GqV9/iaT7giviYsNND2JtUiN54cApYw8tZPOZ0f4u6sxSBpMuj7+dURcmHc8M5F2AK4gTSstxQTgMCJ+WbWYKqxhkimJuYCngYUpeEelqHbgdmAXnxSsXBJLkwoDDOvlU8cDP43g5MpHZXVDOoU0KlVoNKCosfDFH+FrO0Q8XZ3ArOlJeh/uHgbfHFTahfJ0XwBvA6sS8Vl1grO8SSwPPEbhkfLujAdOiuCsykdllZStRbqSdJN3VF1VsJY2AO6k9ERqugnArkTcUvmgKq+R5iYeT+rgPFMidd99MGECjB2bfv797y7PGwxsRWpaZlauX1Lgy2j//eHxx2HiRLj88oLPGwL8MBtVtWaUplYcTqFEas454eabYdw4eOMN2HnmqrCzATvA6bUI05rWxgvANwomUosvDnfeCZ98Au+9BxdcAH3/9xXaF1iQNB3ZmtevKXBuKuHaaQhwssTCNYjRZs1+wCrA3nWWSE2vKtg1kerhuzF7zhV0XhNapxoimcqqpX0PCi/qP+AAGDo0/Sy/fMFDDAGOqF6E1swkFgPWpcDfy7vvwqmnwmWXdX8I4IDqRGd14ECKjZZfdBFMngzzzw+77goXXwwrrvi/h5WeN4JUnc2sHEdSbGrfL34BH34ICy4Iw4fD+uvD6NEd9xgA/KBRLlisdySWAb5KkWu9Eq6dAPatUnhWAVkV2BOBkRFRb62ANibdM+yqh+/GjGiQvp4NkUwBO8AsTdETMMJ3WKxMo4s9cMstcNtt8PHH3T5/ALC3K/w1oVQE4HsUKoQzeDCMHAnHHQfjx8Mjj8Dtt8OoUZ33DNKdRbPekRYjFWMqPL1vySXhhhtg0iT44AO46y5YaaVCe25fxSgtP8Vv9JRmILC/VNIadasxSfMDNwDfi4hX8o6ngKMolEyV/t04FDi6BnHOskZJpkbRzXzf00+H//4XHn443XgrYhqwdRVis+a3M0VGRXshSBc91lzWhSK9xZZbDqZOhZdfnrHtn/8sdDE7ENi1SvFZc9ua7m40nnce7LQTDBoECy0Em2+eEqqZDSWV5bfmsyPdFOsq8dqpL/D1KsRms0BSG6k64xURcUfe8XQh9QfWp9CNntK/GwG+hDRflaKsmEZJpuYv9sDRR8NSS8HCC8Mll8Add6T/X8BAUhVAs96as0LH8eev+cxLsVGB2WaDMWNm3vb552lOTVdzVDguaw3zkr7bCnvwwXSBMmYMvPMOPPFEqh/ZVd1frFhZihZM6sW1U+Dvrnp0OqnVwYk5x1HMnFCkinbvvhsn0wCfv0ZJporG+dhjaf3a5Mlw1VVptHCLLQruqu6OY9aNSnxu/PlrTn0olkyNGwezzz7zttlnT6u9Cx/HrLeKf26kNAp1880wZAjMPXda9H3mmb07jjUyXzs1IUnbk6bm7hoR9dpHtQ/FRs17990YNMDnr+4DzHS/IqWDiPQdUsAkUrdls94q+BfeS9Pw568ZfQJFmoK/9BK0tcEyy8zYtsoq8NxzhfYeU2ijWQ8+Id257WquuVI1vwsvTFfMn3ySSo4WvmL+pJpBWm7GlbpjN9dOgb+76oakFYCLge0jouRr4xx8SrHlEb37buxPA3z+GiWZuoXU82Amw4bBJpvAgAGp2usuu8B66xWaEg6kE8LdVY7TmtOdwNRCD/TtO+Pz1/G/C+hPaqhnzeURin1htLenUYGTT04LbtdaC771Lbj66s57TgEaqtu71Y27KZbMf/wxvPYa7LdfOikNGwZ77AHPPNN5z/HAzVWO0/JxFwU+H728duoHPF7lOK0EkoaS/laPjogn846nWxETgS4nG6A3340AHwDvVjHSimiUZOoKClSk6dcvlaX+73/ho4/gwANh221nXtPWwbMRdO2kYNaz8ygy9/fYY1OPqR/+MBWimTgxbetkKnBDhEcfmk5qdvo7il3Qjh6dFv9/+CFcd126sH3++c57TQXOr2qc1pwingNeLPr4dtvBZpulL8lXXoEpU+DQQzvv1Qe4qopRWn7OIc3KmUkvrp2mAFdGdL2ZbbWVNea9HHgwIrpvxlI/zqTYzJ7SvhvHA2dTT72zilADxAiAxDXATpRX5nMcsGcEv6tsVNYqJJ4CVi3z6e3AWhH8s4IhWb2QVgMeolivn549RsQaFYzIWom0I6kxa+F+Lt2bClxLxJ4VjcnqhsTzwAplPn0CsGpENwm71YSkI0htgtaNiC4Jcl1KFf3+C8ze065FTAAWIKLub0Q3ysgUwLH0Yv5vB5OA54DbKhuOtZjRMLngVL8etAO3OJFqYhFPAXeQ3uveascNnW3W3DwW3irz6moccEJFo7F6sy/porS32oFrnUjlT9IGwOGkdVKNkUgBREwm9Tor57txPHBsIyRS0EDJVARvAJuQhgxLHU6bCLwObBZReM2LWU8kDQM9CY8fBNGbL6V24O/Ad6sUmtWP3UnrCnrzpdEO7EyE1yNY2QTDXk/9XN4kfeeVIkjfpRsT8WbVgrPcRfAg6Tuot99dD+Jm4rmTtAhwLbBbRLyVdzy9FnEVqYx7b78bLwPOrUpMVdAwyRRABI8BawJv021SNTVg6hTgAeDrEXxWoxCtyUhajFRkYPuItS8GbQnxKYwr3Kg1mUS6qLke2DSiSK8Fax7pDtzGpG70EymwTqGDsaTqRJsTcXsNorMmJKmPpDOAm1eBjwfAV4CHSRfNxW4eTk+i3gLWIOKJ2kRreYrgt8C2wOd0W532iymk89eVwFa+CZ0vSQOAm4CfR8S9ecdTtohTgYOA9mndJ1XtpM/ficDBjbBWarqGWTPVkUQf0oXLUcDapAuX6bXo2+CDO+DeiyN2uT+/KK3RSVqFVGXtnIg4d8b2ebaBbX4Ol34CWoFUmjhI6/mmkcqWXhTB23nEbTmTFgdGk6bX9CEVpxCpouNzwBnAbUQ4ybayZBdZlwFLAVtHxEcdHlwJOBTYhZRUTSN9/gaQRhvOAu4lorsbQtaEJAYA2wHHAMsx03fXZMFTv4dvHB7BO3nGaYmki4CFgO2iES/WO5OGPg1HLgkHzJHOR9O/A9tI7UHOBi4nou5LoXfWkMlURxLzkj5sg0l3Xd6MYLyk5YA5I+LRXAO0hiXpeODfEXFDp+0PAhdHxHUSCwHzk04MnwGvRRTp+2KtJS2+XQqYg3TD5wMi6r7Eq9U/ScsDPwb2iYjCd3qlIcDiwDDSHd936Jh0WUuTWASYjxl9fF4DTQX61HEj2JYhaXfS3/jXI+LzvOOpFEl9I93gWQKYm3TD8RPgtUa+wdPwyVQxkrYBjge+1hQZvdWMpFHAGxHxUIHH1gauBpaLCE+BMLOaydZP7AWc4u81qzRJl5BKb1+TdyytTNJwUg+5EZHaHzQFScuS+mR9pdnOXw21ZqqXfg8MAjbMOxBrDEp+BJwCFOssfjRwlhMpM6slSSuTGn+XU5nNrBQPADvnHUQrkzQnqXfhgc2USGV2Au5vtkQKmnhkCv43TLp9RGyTdyxW/ySdCHwL2DIKTMeS9GXgHmDJiF5V9TMzK5ukL5F6mR0cEdflHY81J0mzAe8AS4enhNacpD6kNhsvRUSX7tqNLGs6/DywV0T8Le94Kq3Zk6l+wOBmmm9qlae0tuALYGHgv1Gkr4Gkq4AXIuL0WsZnZq1L0lyk9Zhfjohncg7Hmpyko4FbI8L9pWpM0gmk2VQbRpMVKJI0mFR86WCPTDUgSYsCu0bEGXnHYvVH0nykO0GXRcSvutlvceAp0h27z2oUnpm1qOxO7lGkqTFfjQZenG2NRVJfF6GoLUlbAJcAq0fE+3nHU2nN/plq5jVT030CHJZVPzL7n2wx5F+Bu0gnse4cDvzaiZSZVZukvsCFpPLmWzmRslrJkvjnsh6LVgOSlgIuB3Zq0kSqD/B8NrjRlJo+mYqI8aQvpSPzjsXqzrbAmRFxQnfDzpLmBXYDzqtRXGbW2hYjla1eLyLc88dqJvsufBDYMe9YWoGkQaSCE6dFxMN5x1MlawOTI6Jpe282/TQ/AElzA38h1euflHc8lq+sbP6EiLi7xP1PBuaPiH2qG5mZtTJJ85AaPp/q0SjLi6QNSM3qV807lmaWjQJeDvQDdmvGtUQAki4G3mrm9eZNPzIFEBEfA8OdSJmk/YBfkpoUlrL/UGA/4KxqxmVmrU3S0qRpxwOBpryosobxIPBXSQPyDqTJ7QN8Fdi7WROpzHvA9XkHUU0tMTIF/5uD/gdg54j4JO94rPYkHUK667tZRLxW4nMOA9aICE95MLOqyNanPAqcHBEX5x2PGTR/0YA8SVqDVPxq7Yh4Oe94qqVVPkMtMTIFkL2Z7wD75x2L1Zak/llZztuAtXqRSA0ADiOV8zQzqzhJcwBvA1s7kbJ6kRVF+Ec2Fc0qKKsifCPw/WZOpDKXSfpO3kFUW8skU5mzgAOzvkLWAiQNA/4IHBARr/eyEeGuwL8i4unqRGdmrUzSPsDfgb4R8UTe8Zh18DppyunqeQfSTCS1kaa8XRURt+cdTzVlN7G/RZo22tRaKpmKiBeAK4CmLc9oM0haBHgIeAH4WS+f2xc4Go9KmVmFKTmVVGV264iYmndMZh1la3h+A+ycdyxN5jRgKnBC3oHUwFbAYxHxQd6BVFtb3gHUWkQcJalNUpu/wJret4BrgLPKWNy5LalIxQOVDsrMWt4CwHDStOMPc47FrJhrgZF5B9EsJI0klZxfvRXWEZGSxgvyDqIWWqYARUeSrgb+HBFX5x2LVV5W1nVwRNxZ5vMFPEbq+3BrJWMzs9YlaXbgIOD0FrmYsibQKkUEqknS8qTpblu0wpTebHbPtCavUvg/LTXNr4OrgaOzrszWRCTtAvwWaJ+Fw3wTmA1o6vnMZlY7khYiXUwtDHhRvzUESbsCv8g7jkaWtVi5GfhhKyRSmT2B8/MOolZaNZm4G5gMbJl3IFY5kr5LWuO0YUTcNwuHOgY4000zzawSJM1L6iF1PTDaU8ytgTwIbO+eU+XJZrpcCjwSEZfmHU8N7UwLLZNoyWl+AJJGABMj4u85h2KzKBtOHgjMRVo3+59ZONbqwC3A0hExuUIhmlmLyqb2jSX1q/P3jTUcSQ8CZzd79blqyHpV7gysGxET846nFiQtCDwPLBQRE/KOpxZadWSKiLgfeFaSK/s1sKz05u+AH0XE27OSSGWOBn7mRMrMZpWkHYB/AAOdSFkDOy/vABqRpPWBo4DtWyWRysxNWnPeEokUtPDIFICk7wHbRYSn+zUgSfOQOoi/Cuw1qwmQpOWAR4AlI2JcBUI0sxaUTe05NPvZMiKeyTkks1mSfab7uBBFaSQtDDwO7BERd+cdTy21YsGSlh2ZylwDDJf0lbwDsbJsBdwPjKrQSNKRwEVOpMxsFs0BbEoqfe5EyprBL4Gd8g6iEUjqD9wIXNiCidQywFNZ8t0yWnpkCkDSUcDwiNgl71isNJK+BiwSEbdU8JgLA88Cy0XER5U6rpm1DkmDgENI60um5ByOWcVkVf12joit8o6l3km6AFiUNPOppQpZSToWWCAiDsg7llpq9ZEpSHdbzso7CCuNpK2APwCVHkI+BLjKiZSZlUPS3KRKsV/B363WfG4D1s2m11sRknYjjUrv0YKJlIBdgd/kHUuttfwJPyLGAG9I2ibvWKx7Wffw/wO2rmRVIUlzAt8DzqnUMc2sdWQV+x7JfnaNiEk5h2RWUdn09zMAJ1NFSFoFOBcYGRGf5x1PDgYDfwH+lncgtdby0/wAJM0PvACsEBEf5B2PzSy72zEYGAoMjYiXK3z8HwPLRsSelTyumTU/SbNFxDhJ60TEw3nHY1ZNJ0l9T4D1gKVJze3HkK6f/k4LX1BmN2UfB46LiOvyjqfaJJYDvkFaHzoZeA+euTviK+25BpYTJ1MZSRcBn0fEj/KOxWaQ1I80GjU2Ig6swvEHA68DG0TE85U+vpk1L0mbAZcAK7fonWhrFdLcAd/7CE6bGyZl05ragKlAAB8BPwWuIWJsbnHmQFIf0jTIVyPikJzDqRqJNmBrUrn3VUjLLdqAaRBToX02mHA5zHN6BK/mGWutOZnKSFqSNDS5uKdo1Ids6sxNwCRgp4gYX4XX2B/YOCK2rfSxzax5SdoL+Alpkflf847HrGqkbwK3An1Js0SKGU/6vt6YiKdqEFldkHQcaZ3UBs1aeEZiPuBeYAnSiGQxU0hJ1jERnF+D0OqCk6kOJM0VEZ/kHYclkrYHNgIOiIipVTh+P+BlUqLmhppmVhJJQ0g97vaLiBfzjsesatLo683AoBKfEaSk6ptEPF61uOpENjp9KbB6RLyXdzzVIDEv8DQwH9CvxKe1A2dGcHLVAqsjTqY6kNSXNEz94xbrVl1XJK0IrBQRN1b5dXYFfhARI6r5OmbWHLIbMAcDF3gGgzU9aXngCWBIGc/+DFiBiPcrGlMdyWY0/R3YPiIeyjueapDoQ0qklgf69/Lp7cAeEdxU8cDqTMtX8+so69i8PLB73rG0KknrAfcBA6v8OgKOAU6v5uuYWXOQNBtwO7ABaZ2AWbP7Md19F++4Izz/PIwbB6+8Auus0/HRgUDT9hrKesrdBPykWROpzCbAkhRIpJZfHu69Fz77DF5+GbbdtstzBwNnSTR9A18nU12dARwlqS+SkOZFWhZpMVKxAqsSSZuQTk67RsTVlTkmkphHYhmJxaX/3WHbgrRw9s+VeB0za16SBgL3A28D36rG+k2zuiLNAWxPWifV1UYbwZlnwne/C0OHwnrrwWuvddxjILA/aTS3qWQ3Yy8CXgJ+nnM41XYUqZLyTPr2hdtug9//HuaaC/beG665BpZdtsvz5wHW6bK1yXiaXyeStADccxf8YxXYGZiLtKBOpLmi95Ka/N7fymVAKyk7MQ0i3cVYKCKemfVjMgewB3AEMC8z3sP+wH0wciG447SIydfP6muZWfOSNDgi2iWNAB4In/etFUgHkwqsFL6J/MgjcOmlcNll3R1lLLAXEU01zUvS3sBBwDey/ltNSWIx4EUKjE6utBL8/e8pj57uT3+CRx+F44+fadcA7ojgW9WNNl8emepI6hNwxruw1iqwL7AgMIBUuWQI6UJ8M9JUj9eRhucWa5PI1qmdC/wyIj6a1UQqG4k6CXgPOA1YhJnfw34QG8OVK8GksyW+Nou/gpk1KUnrAC9Kmi8i7nciZS1kfYolUn36wOqrw7zzpvldb78NF1wAA7tccw8F1qhynDUl6evAqaQqnk2bSGVWJfWQKokEX/5y18002WegECdT06WL+puA/ZWy8GJT+kS6MF8MeBhp/RpF2HSyOce/BYaTFnXP4vHoA1wLHE56D4ssmpVgNoEWBu6X2GRWX9vMmoukkaQqZt+PiA/zjsesxuYq+sj880P//rD99rDuujB8OKy6Khx7bKG9561WgLUmaV7gRlLhqpfyjqcGhlEkT3jxRfjwQzjySGhrg403hvXXh8GFr5zLKWDSUJxMzXAhaaFdqW+6sn1/j7RS1aJqbuuT7npsGhGfVuB4ZwHb0Ls/3MHAzRLDK/D6ZtYEJPUnLZ7fNCL+lHc8ZjmYUPyR7KELLoD334ePP4ZzzoEttii0d1OsL5TUBlwHXBsRt+UdT41MJE3T62Lq1FRwYsst00fg8MPhhhvgP/8peJySR7calZMpIEuG9qDYRXj3FWuGkBIxK5GkJSTtHBF3kYpNzHKJYYmlgdF0eg/Hjp35Z+pU+HnX5aKDgV/Nagxm1tgk9ZF0EKla3zcj4um8YzLLyWuk5qtdffZZmtrXcdZr4Rmwk4DXKx9aLk4hJRbH5R1IDb3T3YPPPgsjRsA888Bmm8FSS8FjjxXctWnL40/nZCo5hGKNyHquWCPgG0hLVD3KJiDpq8AjwNwAFVyDcCAFqg4NHTrjZ4EF0g21G7t2rxKwssTyFYrFzBpMVrHvOlIFswFeH2Ut7lJSMlTY5ZfDgQemdVNzzAGHHppKu80sSH9TDU3St0kFyXbOWui0ir+RekUVtPLKMGAADBqURqYWXBCuuKLLbuNpgQEHJ1PSUGBXivUNOekkOPnkVKIkAt59N/3MrA9N3E+hUiStBfwROCAiKvbHJTEI+B49dOYeOTLN8X2ocEeINiqwbsvMGk82hecu0o2VTSo07discUU8RXejSqecAo8/Di+9BC+8AE8/DaedNuPpKZG6n4huRzfqnaTlSDNXvhMRH+UdTy1FMA04B6YVnPI5ahS89166rtpww7RuanLXCX19gGuqHGruXBpd2hy4Hpi9y2N9+qShjOOPh+9/P1WqufXWtOJu4sTOe79FxOLVD7gxZXd9BwLLRUThgeCyj80GwK0Ueg87uPdeePDBlB8X8WEE81cyNjOrb5IGRsRESRsC90XEtLxjMqsL0q6kRKLXBQTaIT6FLRZO0/kbUtao+1HgvIj4v7zjyYP0tVHw4JUwqJzGuxOBayP4fqXjqjcemUrTzQp/SHpXsWZYFWNsWEqOB66LiM8qnUhl5u5ph8UWS5Vmrryy2926NKYzs+al1N7i35IWj4h7nUiZzeQ3wB/orhhFYe194eeLwJ8knZiN7jSUrP/lr4G/tW4ipcPgiePhzmPp/WdgKvAWcGjlI6s/Tqa604uKNZ/BMEknAkh6V1JkP09m2y7psC0kLSRp607b9s727bjtjmzbHR23Z9v27rTv1tlxO267JNv3yQ7b3s22ndhp369mPx23zdLvBDxIqrC3bbV+J9j+RhjTbSI0ahQ8/DC88UZ5HwUzay6SNgb+DBwZEW/mHY9Z3UlTl3YD/kTpVfnGA1cNgMOy//8R8FdJh0hqpGvOg4FlaMElHJK2kLQgaXre8IjtfwLsR+kJ1UTgVWD9CMZWKcy64ml+3U3zA3jrLfjxj+Hqq9P///a34bjjYLXVuuzpaX4zk7QucCSwSzWb25Uyze/FF+GMM9Ka2W54mp9ZC8iahf8ZODEiCq+iNLMkJUEHAT8EBtF1FkeQkqiPgeOJuGrmp2sZ4ArgFxHxm6rHO4uya5cbgW9ExBs5h1MzkoYB5wIbACMjrZvr8DjrAmcDK5PWmXdepz6ONNPrMuBHETR7U+P/cTKVClB8QDpBdHXSSbD55qmY/pQpcPvtcP/9aR3VDJOBC4g4ourxNgBJCwAbRsS1klTtqlhZAYoPSc2Uu1hzTbj77lTNb1zxP+0pwKUR7FedKM0sb9nUnR+QKoyNc8U+s15ISdWmpArIy5DWUo0D/gWcAzxEkb+p7AYGwGbA4sAv63FaraSFgMeB70bEn/OOp1ayIjzPkGYTHRkRRUeUssrHBwPfJN3Enkwqf34x8NuIXk8JbHhOpgCk/wP2pFBFv7Y2OP982GWXVHTihhvgqKNg0kwVQycCK9BCdzCKkfQlUsW+yyPilNq9LueR+kx1qej3y1+mrty7797tISYAq0bwYlUCNLNcZRcLFwFfBzaPiKbvfWJWbyQtD1wJjAX2ioi3cg7pf5Sadd8H/CEiTutp/2aQFdnYISIuk7RQRHQpV209czIF05v2Pk6x0anuBfAAERtUNqjGI2kV0tzqH0ZE9xPqKv7aLE26OzawjKcH8HgEa1Q2KjOrB9lajVuBAcD23d11NbPqym5sHAmsFxGb5x3PdJLOB5YEtq3HUbNKk7Q+cDnwALB3REzJOaSG5WRqOuliYBS9LwE6DliDiOcrH1TjkDSANCr01Yh4IJ8Y+BmwD71/D8cD60Twj4oHZWa5kjQgIiZJ2oRU+twXDGZ1IJv6N4S0Tuf4yLEnlaRdgJOB1SPis7ziqBVJG5FGCPeJiC7dlq13GqmySrUdQFqQXGrFmukLLrdyIqWDgN9HxLi8EqnMkcDtdNOxu4B24NtOpMyaT1aS+V+SVoyIPzuRMqsfEfEFaYr9W8DTkkZl6xprStJXgPOB7Zo9kZK0lqQNSNMZV3YiVRlOpqZLf9TbAxeS1kAVuyAP0mjUW8A65Js85EpSH0lnk0pm/iDveLJu3buSqs1MoHhiPP09/A8wIoK7axOhmdWKpDVJi6nPjBa/4WVWryJiSkScRCpMcSiwbC1fX9IcwO+AQyLimVq+di1JGijpp6TfdXBEfBERn+QdV7PwNL9CpDmB7wKHkxrCTiElnv2Au4GzSOukWvofT9JqpMRl+3r7o5QYBuxBGq2aj1RtZvp7+BfSe/iXCFr6PTRrRtnd7T8CP4+IP+Qdj5n1bHr1X0knAP+OiN9W+fWmr6V8IyIOquZr5U3S9UBfYHRE/DfveJqNk6nupC/kuYE5gEnAx0T0ZgpZU1JKNreKiKtrUfp8VkgImAuYk5RQfRxR8lROM2swknYD7gQ+q+dzk5kVJunrpPU8z5Iu/j+q0uv8GNgC2CAiJlfjNfKUrWU/GPgFqVr15z4nVoen+XUnIoj4iIhXiHjbiRRIWgx4GFi13hMpgAgigo8jeCWCt5xImTWnbNrxGcBxwOz1fm4ys8Ii4jFgNdJyigOr8RpZQZrRwHeaNJEaTqpSvTYwMCJ8c6mKPDLVyqTBpLsyC5NKin8OPEHEE4V317KkRYs/i4hzaxanmbWUbER5deCrwDDSOtZ3gDsLNYTMpvVdQWokuk1EfFy7aM2sWrK/7dVIjYIP7mlJgUQ/UmPhJUmVAscAzwEPTp/WL2kJ4O+k/koPVi34WZFa9qxNmhk1BfgAuJOIz3t+qhYAngaOBq52ElV9XZvUWvNLSdHBpEbF04D+pLm0U4BAehs4E7ieiAnpKepPKtiwVyt1BTez2pEYDOxEughYmBnrHL8gTdOVxGXABRG8kp6j/hExWdK1wEORnbPMrPFla6j+DXwKPCtp74i4s/N+EguTimHtTzpv9Cdd404hnT8+kzgbbrgOuIlUmKa+EimpH/Bt0vlvBVKxrP6k67RJQBvSDcA5FCiWIenLwPoRcZGkZSLCM3FqxCNTrUban1R8oS/pj7SYccAnwAjBuqT+Tev4DoeZVUPWePt+0p3Y2brZdTLp4ugw0J2k9VHfz6YGmVmTkjQCOB7YPCImzdjOtsC1pCRqYDeHGA/j2uDbD8A9m9XV9Yw0H3AvsATdn/+mkhLEs4ETSMlmG3AEqWjaMRFxaZWjtU6cTLUS6SjgBGBwKbsHfDEBJn8FPn4VNomIF6oboJm1IomlgCdIU/pKXMv7xUQ4fiL85FTgnLq6MDKzqsmq8N0M/AJiDuBySryuSYM8Gg8aEUHBJQ01J81LmpY3H2kkvhTjgcuJOFDSEaQlG3tFxBvVCdK642SqVUhbAjcCg3rztC8gAt5pg2XocCfIzKwSJAYCLwML0euiSFMnQ9u3IrirCqGZWZ1KBSS+fiU8OA8MKGfJyqfA8hF8WOnYeiWtCXsSWInuZwt1MQUmPwY/XQdOBaZExLRqhGg9czW/1nEmhRKpOeeEm2+GcePgjTdg551nergvqC1NuxlZiyDNrOXsQDrHzPR9NHbszD9Tp8LPf975qW39gTNqE6aZ1Yu0dvvBp6Ff386PlXbuYCBp+ULeNiA1Ku6aSC2+ONx5J3zyCbz3HlxwAfSd8ev2g/5rwH4BU51I5cvJVCtIzXWXLPjYRRfB5Mkw//yw665w8cWw4oqd95qNtCDSzKzSjqbAGoGhQ2f8LLAATJgAN95Y8PnLSnyl2kGaWf2QWAAGbAB91PmxEs8dg4CDJbokYzV2JKnqYFe/+AV8+CEsuCAMHw7rrw+jR8+0S1tKwjatepTWLSdTreEQYECXrYMHw8iRcNxxMH48PPII3H47jBpV6BjLIPmCxcwqRmJV0oLrbo0cma4pHnqo4MP9gUMrG5mZ1bkflLJTCeeOLSsZVK9ICwEjgC4JIQBLLgk33ACTJsEHH8Bdd8FKK3XeayhwVFXjtB45mWoNa0KBuy/LLZfGv19+eca2f/6z0B8rpFWbq1YpPjNrTauRyv92a4894Kqrij7cRjrHmVnrWJfuK/cBPZ47ZiPf65qVSSXPCzvvPNhpJxg0CBZaCDbfPCVUXa1SpfisRE6mWsPQgltnmw3GjJl52+efp7HxrtpI6xrMzCplGD1Ur1pssTS75coruz1O4XOcmTWrOXvaoYRzh4B5KhlULw2j2KgUwIMPppvbY8bAO+/AE0/ArbcW2rPESoZWLU6mWkPhOx/jxsHss8+8bfbZ04rNrqYB7ZUOzMxa2gRSz6iiRo2Chx9O9XG6MbGCMZlZ/evxeqS0c8etwwEkvSspsp8ns22XdNgWkhaStHWnbXtn+3bcdke27Y6O27Nte0///1vDb6cUmjWUdkyjUDffDEOGwNxzp4JhZ55ZaO/JPf1bWHU5mWoNbxXc+tJL0NYGyywzY9sqq8BzzxXaeyrwn2oEZ2Yt621SA8qidt+9x1EpKHaOM7Nm9SrpJm9RJZw72mHb3wBExEIRoeznq9m2vTtsU0S8GxF3dNp2SbZvx21bZ9u27rg923bJ9P9/B6xadFh+rrlSNb8LL0xFwj75BC6/HLbYotDeH3T7W1rVOZlqDRcCXYeb2tvTXY+TT07FKNZaC771Lbj66kLHmAbcU+U4zay1/Lm7B9dcExZeuGgVv+nGks5xZtY6fkUa2S6oxHNHH+CGCsfVG/+EIn2uPv4YXnsN9tsvlUMfNiwtAHvmmc57tgMXVTlO64GTqdZwC8Xu4IwenRY3fvghXHdd+sN9/vnOe00CLiKi2zvIZma9EcFk4GKKTEXeY48ZbfC68QVwe+WjM7M69hjwbrEHSzh3TANuj+DjKsRWmogg9QAdX/Dx7baDzTaD//4XXnkFpkyBQ7sULu0DXFHVOK1HSu+lNT3pVOAwCjXu7dkEYFki3qlsUGbW6iQWA/5N+eemsyI4obJRmVm9k9iTNCpduE9T99qBERE8XtGgekuaDXiPAr32SjAJuIGI3SsblPWWR6Zax8nAM3RXhrOwdmAvJ1JmVg0RvAXsQ+8L3EwCngZOq3hQZtYIrgT+SDfT/YoYD5yWeyIFEDEO2Jbe/w5TSWtFD6h0SNZ7TqZaRcRkUpfspyn9j3YCcDAR11ctLjNreRFcTRo5L/Xc1A48CWyeTRU0sxYTQQC7AndRbKpcV+3AecDpVQqr9yLuBXYmxVbKdLGJwCvAekSM6Wlnqz4nU60k4nNgfeAM4FMKFaVIlbUmAH8DNifi17UL0MxaVQS/ArYAHiWdgwqt0RwLfEy6EBoRgS8kzFpYdjNle+Bo0hqqsXRNSL4gJSrPAbtGcGyWiNWPiNtIjYjvJSVLhWYRjc1+LgK+RsT7tQvQuuM1U61KagO2BkYDi5I6iX8OPAT8nIiXcozOzFqYxJeAg4G1SY0tJ5KmtFwE3BnB1BzDM7M6JCHgm8BBwLKkZrZjgaeA8yN4KsfwSictCuxPmk00B+nG0nvAJcBNRPR2uYZVmZMpMzMzMzOzMnian5mZmZmZWRmcTJmZmZmZmZXByZSZmZmZmVkZnEyZmZmZmZmVwcmUmZmZmZlZGZxMmZmZmZmZlcHJlJmZmZmZWRmcTJmZmZmZmZXByZSZmZmZmVkZnEyZmZmZmZmVwcmUmZmZmZlZGZxMmZmZmZmZlcHJlJmZmZmZWRmcTJmZmZmZmZXByZSZmZmZmVkZnEyZmZmZmZmVwcmUmZmZmZlZGZxMmZmZmZmZlcHJlJmZmZmZWRmcTJmZmZmZmZXByZSZmZmZmVkZnEyZmZmZmZmVwcmUmZmZmZlZGZxMmZmZmZmZlcHJlJmZmZmZWRmcTJmZmZmZmZXByZSZmZmZmVkZnEyZmZmZmZmVwcmUmZmZmZlZGZxMmZmZmZmZlcHJlJmZmZmZWRmcTJmZmZmZmZXByZSZmZmZmVkZnEyZmZmZmZmVwcmUmZmZmZlZGf4fpfjM8bFtWsUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1080x288 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Computed max cut for two subgraphs\n",
    "strr1 =  '01010101xx'\n",
    "strr2 = '0xxxxxx101' \n",
    "strr = '0101010101'\n",
    "\n",
    "# Show graph illustration\n",
    "options0 = {\n",
    "    \"node_color\": [\"red\" if strr1[i] == '0' else \"blue\" for i in S[0].nodes],\n",
    "    \"style\": [\"solid\" if strr1[u] == strr1[v] else \"dashed\" for (u, v) in list(S[0].edges)]\n",
    "}\n",
    "options1 = {\n",
    "    \"node_color\": [\"red\" if strr2[i] == '0' else \"blue\" for i in S[1].nodes],\n",
    "    \"style\": [\"solid\" if strr2[u] == strr2[v] else \"dashed\" for (u, v) in list(S[1].edges)]\n",
    "}\n",
    "options2 = {\n",
    "    \"node_color\": [\"red\" if strr[i] == '0' else \"blue\" for i in range(n)],\n",
    "    \"style\": [\"solid\" if strr[u] == strr[v] else \"dashed\" for (u, v) in list(G.edges)]\n",
    "}\n",
    "\n",
    "fig, ax = plt.subplots(1, 3, figsize=(15,4))\n",
    "for i, a in enumerate(ax):\n",
    "    a.axis('off')\n",
    "    a.margins(0.20)\n",
    "nx.draw_networkx(S[0], pos=nx.circular_layout(S[0]), ax=ax[0], **options, **options0)\n",
    "nx.draw_networkx(S[1], pos=nx.circular_layout(S[1]), ax=ax[1], **options, **options1)\n",
    "nx.draw_networkx(G, pos=nx.circular_layout(G), ax=ax[2], **options, **options2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### DC-QAOA\n",
    "\n",
    "To find the max cut for a large graph with the number of vertices exceeding limits, DC-QAOA recursively divides it through LGP policy and conquers sub-solutions with GR policy described above. It adopts the divide-and-conquer paradigm to deal with the max-cut problem for large-scale graphs. \n",
    "\n",
    "The input graph is separated into exactly two subgraphs through LGP policy if the vertex size is larger than the qubit limit. Otherwise, its max cut would be directly calculated by QAOA. Each subgraph will recursively call DC-QAOA until its max-cut solution are returned, i.e. at some step, the vertex size of the subgraph is under the limit and can be directly processed by QAOA and returned values would be combined and then be given to the upper layer. \n",
    "\n",
    "The code below provides the way to run DC-QAOA and LGP and GR policies are both applied in the `DC_QAOA` function. Note that QAOA would return a series of possible cuts, sorted by frequency count.   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-11T06:54:04.788730Z",
     "start_time": "2021-05-11T06:52:49.025846Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "First t possible approximate maxcut for graph G: {'0101011011': 408, '0101010011': 393, '1010101100': 386, '1010100100': 350, '1001010100': 257, '0101010101': 250, '1011010100': 248, '0101001011': 244, '0110101011': 235, '0100101011': 225}\n",
      "Max cut found by DC-QAOA algorithms: 0101011011\n"
     ]
    }
   ],
   "source": [
    "def DC_QAOA(g, p, t, s, k, ITR, LR):\n",
    "    if len(g.nodes) > k:\n",
    "        # get exactly two subgraphs with LGP policy\n",
    "        S = NaiveLGP(g, k)\n",
    "\n",
    "        S_str_cnt = []\n",
    "        for si in S:\n",
    "            siv = list(si.nodes)\n",
    "            # Compute the subgraph's maxcut recursively\n",
    "            _, si_str_cnt_relabeled = DC_QAOA(si, p, t, s, k, ITR, LR)\n",
    "            # refilling str_cnt1 and str_cnt2\n",
    "            si_str_cnt = []\n",
    "            for str_relabeled in si_str_cnt_relabeled:\n",
    "                strr = \"\"\n",
    "                for v in g.nodes:\n",
    "                    if v in siv:\n",
    "                        strr += str_relabeled[siv.index(v)]\n",
    "                    else:\n",
    "                        strr += \"x\"\n",
    "                si_str_cnt.append((strr, si_str_cnt_relabeled[str_relabeled]))\n",
    "            si_str_cnt.sort(key = lambda tup:tup[1])\n",
    "            S_str_cnt.append(si_str_cnt[::-1][:t])\n",
    "        # Reconstruct string-count map with GR policy\n",
    "        out_cnt = GR(S_str_cnt[0], S_str_cnt[1])\n",
    "    else:\n",
    "        if len(g.nodes) == 1:\n",
    "            return [(\"0\", 99999), (\"1\", 99999)]\n",
    "        _, out_cnt = find_cut(g, p, ITR, LR, shots=3000)\n",
    "        # Transform {str:cnt} dictionary into [(st, cnt)] tuple list\n",
    "        out_cnt = [(k, v) for k, v in out_cnt.items()]\n",
    "        # Sort string-count map by counts in reverse order\n",
    "        out_cnt.sort(key=lambda tup:tup[1])\n",
    "        out_cnt = out_cnt[::-1]\n",
    "\n",
    "    # retain only top t (str,cnt) pairs by sorted order\n",
    "    out_cnt = out_cnt[:t]\n",
    "    # resacle total numer of counts to s or around\n",
    "    cnt_sum = sum(cnt for (str, cnt) in out_cnt)\n",
    "    out_cnt = [(k, int(s * v / cnt_sum)) for (k, v) in out_cnt]\n",
    "\n",
    "    return out_cnt[0][0], dict(out_cnt)\n",
    "    \n",
    "# Set QAOA parameters\n",
    "p = 2     # Number of layers of QAOA circuit\n",
    "ITR = 100 # Iterations of the training network of QAOA\n",
    "LR = 0.5  # Learning rate in the training network of QAOA\n",
    "\n",
    "#Set DC-QAOA parameters\n",
    "s = 3000  # Multiplier to make frequency bigger\n",
    "t = 10    # Number of partition strings kept after graph reconstruction \n",
    "k = 5     # Maximum qubits/vertices limit\n",
    "\n",
    "# Using DC-QAOA\n",
    "max_cut, out_cnt = DC_QAOA(G, p, t, s, k, ITR, LR)\n",
    "print(\"First t possible approximate maxcut for graph G: \" + str(out_cnt))\n",
    "print(\"Max cut found by DC-QAOA algorithms: \" + str(max_cut))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Applicability"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**DC-QAOA described above can approximate max cut for a graph if and only if its and one family of its recursive children's pseudo-connectivities are all smaller than $k$, where $k$ is the qubit limit (vertex limit).**\n",
    "\n",
    "The pseudo-connectivity here is defined as the minimum number of vertices that need to be removed to separate the remaining vertices into exactly two isolated subgraphs. A graph's recursive children mean its two partitioned subgraphs and then their partitioned sub-subgraphs until some are smaller or equal to the qubit limit. \n",
    "\n",
    "Cycles are applicable examples if the qubit limit is larger than $1$. Its pseudo-connectivity is $2$ as removing two vertices can partition cycles into two paths. For the example of $C_6$ (cycle of six) with qubit limit 4, partitions described above are $P_5$ (path of five) and $P_2$. $P_2$ is under the qubit limit, and $P_5$ can be partitioned into $P'_4$ and $P'_2$ and thus have pseudo-connectivity $1$. Then $P'_4$ and $P'_2$ are under the qubit limit. So this family of the $C_6$ has all its recursive children's ($C_6$, $P_5$) pseudo-connectivities smaller than $4$. \n",
    "\n",
    "Non applicable examples include complete graphs with $n$ larger than $k$, because no matter how you remove vertices, the rest graph is still connected.  Another example for which different qubit limits influence its applicability is shown below. \n",
    "\n",
    "The left graph is the original graph. If the qubit limit is $2$, then the number of separation nodes can only be $1$ and no vertice could partition the graph into $2$. If we remove vertex $4$ (as shown in the middle), the graph will be partitioned into three while others will leave a connected graph. Besides, the graph's pseudo-connectivity is $2$, not smaller than $k=2$, and so DC-QAOA fails in this case. However, if the qubit limit is $3$, then its pseudo-connectivity is under $k$ and we can remove vertex $0$ and $3$ as shown in the right. Then the rest  are two components, with one under the qubit limit and one can be partitioned into two further (by removing vertex $4$ and with pseudo-connectivity being $1 < k$). So DC-QAOA succeeds here. \n",
    "\n",
    "![limitations-connectivity](./figures/dcqaoa-fig-applicability_example.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "People can try the example shown above by adjusting the parameter $k$ in the code below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-11T06:54:45.016920Z",
     "start_time": "2021-05-11T06:54:04.825816Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "First t possible approximate maxcut for graph G: {'11010': 535, '01001': 515, '00001': 495, '11110': 495, '10110': 488, '00101': 472}\n"
     ]
    }
   ],
   "source": [
    "G = nx.Graph()\n",
    "G.add_nodes_from([0, 1, 2, 3, 4])\n",
    "G.add_edges_from([(0, 4), (1, 2), (1, 4), (2, 4), (3, 4)])\n",
    "k = 3\n",
    "_, out_cnt = DC_QAOA(G, p, t, s, k, ITR, LR)\n",
    "print(\"First t possible approximate maxcut for graph G: \" + str(dict(out_cnt)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Performance\n",
    "\n",
    "We have compared DC-QAOA with classical max-cut approximation algorithms to test its performance. We provide the performance below of finding max-cut by DC-QAOA and the upper bound for max cut provided by SDP. We take five random graphs of 10 vertices as an example to show the performance (qubit limit for DC-QAOA is set to be 5). \n",
    "\n",
    "In order to run the following test, users need to install cvxpy: `pip install cvxpy`. **Windows users may encounter an error at runtime if they install cvxpy using pip. Instead, we recommend Windows users to create a new conda environment and install cvxpy with conda.** For more details please see [https://www.cvxpy.org/install/](https://www.cvxpy.org/install/)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-11T06:54:47.679276Z",
     "start_time": "2021-05-11T06:54:47.341106Z"
    }
   },
   "outputs": [],
   "source": [
    "import cvxpy as cvx\n",
    "import networkx as nx\n",
    "\n",
    "def sdp_solver(G):\n",
    "    \"\"\"\n",
    "    Use SDP to find the upper bound for the max-cut problem.\n",
    "    \"\"\"\n",
    "    n = len(G)\n",
    "    adj_mat = nx.adjacency_matrix(G).toarray()\n",
    "    Y = cvx.Variable((n, n), PSD=True)\n",
    "    cut_size = 0.25 * cvx.sum(cvx.multiply(adj_mat, 1 - Y))\n",
    "    problem = cvx.Problem(cvx.Maximize(cut_size), [cvx.diag(Y) == 1])\n",
    "    opt_val = problem.solve(cvx.SCS)\n",
    "\n",
    "    return opt_val"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-11T07:02:35.277145Z",
     "start_time": "2021-05-11T06:54:47.682769Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of node = 10\n",
      "Node = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]\n",
      "\n",
      "Random graph 1\n",
      "Edges = [(0, 1), (0, 5), (0, 6), (0, 7), (1, 3), (1, 4), (1, 6), (1, 8), (1, 9), (2, 3), (2, 4), (2, 6), (2, 8), (3, 5), (3, 6), (3, 8), (4, 5), (4, 7), (4, 8), (5, 8), (5, 9), (6, 7), (6, 8), (6, 9), (7, 8), (7, 9)]\n",
      "SDP upper bound: 21.03787674994891\n",
      "DC-QAOA node partition: 1001100011, max cut = 21.0\n",
      "\n",
      "Random graph 2\n",
      "Edges = [(0, 1), (0, 2), (0, 4), (0, 8), (0, 9), (1, 2), (1, 3), (1, 5), (1, 6), (1, 7), (1, 8), (2, 3), (2, 6), (2, 8), (3, 8), (4, 8), (4, 9), (5, 9), (6, 8), (6, 9)]\n",
      "SDP upper bound: 16.028382190698274\n",
      "DC-QAOA node partition: 1011011100, max cut = 14.0\n",
      "\n",
      "Random graph 3\n",
      "Edges = [(0, 2), (0, 5), (0, 6), (0, 7), (1, 2), (1, 3), (1, 4), (1, 5), (1, 6), (1, 7), (1, 9), (2, 3), (2, 5), (2, 6), (2, 9), (3, 4), (3, 6), (3, 7), (4, 6), (4, 7), (4, 8), (5, 6), (5, 7), (6, 7)]\n",
      "SDP upper bound: 17.0971057189657\n",
      "DC-QAOA node partition: 0010111100, max cut = 15.0\n",
      "\n",
      "Random graph 4\n",
      "Edges = [(0, 1), (0, 2), (0, 3), (0, 4), (0, 6), (0, 7), (0, 8), (0, 9), (1, 7), (1, 8), (1, 9), (2, 3), (2, 6), (2, 7), (2, 9), (3, 4), (3, 5), (3, 7), (3, 8), (3, 9), (4, 7), (4, 8), (5, 6), (5, 8), (5, 9), (6, 7), (6, 9), (7, 8), (7, 9)]\n",
      "SDP upper bound: 20.153931079390457\n",
      "DC-QAOA node partition: 1111101000, max cut = 17.0\n",
      "\n",
      "Random graph 5\n",
      "Edges = [(0, 1), (0, 2), (0, 6), (1, 2), (1, 4), (1, 5), (1, 7), (1, 8), (2, 3), (2, 4), (2, 5), (2, 6), (2, 8), (2, 9), (3, 6), (3, 7), (3, 8), (4, 7), (4, 9), (5, 7), (5, 8), (6, 7), (6, 9), (7, 9), (8, 9)]\n",
      "SDP upper bound: 18.82486962324589\n",
      "DC-QAOA node partition: 1010000110, max cut = 18.0\n"
     ]
    }
   ],
   "source": [
    "n = 10\n",
    "iter = 5\n",
    "print(\"Number of node = \" + str(n))\n",
    "print(\"Node = \" + str([i for i in range(n)]))\n",
    "\n",
    "value_dc_qaoa_ls = []\n",
    "ubound_sdp_ls = []\n",
    "for i in range(iter):\n",
    "        print(\"\\nRandom graph \" + str(i+1))\n",
    "        # Generate random graph\n",
    "        G = nx.erdos_renyi_graph(n, 0.1, 100 * i, directed=False)\n",
    "        while nx.is_connected(G) == False:\n",
    "            G = nx.erdos_renyi_graph(n, 0.5, directed=False)\n",
    "        print(\"Edges = \" + str(list(G.edges)))\n",
    "        # SDP upper bound calculation \n",
    "        ubound_sdp = sdp_solver(G)\n",
    "        ubound_sdp_ls.append(ubound_sdp)\n",
    "        print(\"SDP upper bound: \" + str(ubound_sdp))\n",
    "\n",
    "        # QAOA parameters\n",
    "        p = 2     # Number of layers of QAOA circuit\n",
    "        ITR = 100 # Iterations of the training network of QAOA\n",
    "        LR = 0.5  # Learning rate in the training network of QAOA\n",
    "        # DC-QAOA parameters\n",
    "        s = 3000  # Multiplier to make frequency bigger\n",
    "        t = 20    # Number of partition strings kept after graph reconstruction \n",
    "        k = 5     # Maximum qubits/vertices limit\n",
    "\n",
    "        try:\n",
    "            cut_dc_qaoa, out_cnt = DC_QAOA(G, p, t, s, k, ITR, LR)\n",
    "            cut_dc_qaoa1 = [\"solid\" if cut_dc_qaoa[u] == cut_dc_qaoa[v] else \"dashed\" for (u, v) in list(G.edges)]\n",
    "            value_dc_qaoa = cut_dc_qaoa1.count(\"dashed\")\n",
    "            value_dc_qaoa_ls.append(value_dc_qaoa)\n",
    "            print(\"DC-QAOA node partition: \" + str(cut_dc_qaoa) + \", max cut = \" + str(float(value_dc_qaoa))) \n",
    "        except Exception as e:\n",
    "            value_dc_qaoa = 0\n",
    "            value_dc_qaoa_ls.append(value_dc_qaoa)\n",
    "            print(\"DC-QAOA fails with error message '\" + str(e) + \"'\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-05-11T07:02:35.547590Z",
     "start_time": "2021-05-11T07:02:35.280626Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEWCAYAAABhffzLAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAABR60lEQVR4nO3dd3hUZfbA8e9JIwmBUEJNAqG30EJoCkgRFQGRFn5YALEXxHXV1V3X3ta+9gaCiohUBUFXmgiikITQezMJJbRQ0sv5/XEnEEPKpEwmybyf55knM7eemSRn7n3ve88rqophGIbhOtycHYBhGIZRvkziNwzDcDEm8RuGYbgYk/gNwzBcjEn8hmEYLsYkfsMwDBdjEr9hOICItBGRGBE5LyIPOjsew8jNJH6j2ETkkIiki0hAnumbRERFJMRB++0hIktFJFFETovIBhG5zc51Z4jIC0UsoyKSJCIXRCReRN4UEfcShvsYsEpVa6jqOyXchmE4hEn8RkkdBMbnvBCRjoCvo3YmIr2BlcAvQEugLnAvMKSMd9VZVf2AQcBNwJ3FjNPD9rQpsL0kAeTahmE4hEn8Rkl9CUzI9Xoi8EXuBURkqO0s4JyIxIrIM7nmjRORgyJS0/Z6iIgcE5F6BezvNWCmqv5HVU+qJUpVI2zrTxKRtXn2ryLSUkTuAm4GHrMdzS8u6s2p6i7gVyDUtq1htqabRBH5TUQ65drPIRH5h4hsAZJEZCUwAHjPtr/WIuIvIl+IyAkROSwiT4qIW67Y14nIWyJyCnjGdobygYgss21jnYg0FJG3ReSMiOwSka65YnhcRPbbmpZ2iMjIXPMmichaEXndtu5BERmSa34dEflcRI7Y5i/KNW+E7X2fs23/Ott0fxGZJiJHbWdHL5Ti7Mgob6pqHuZRrAdwCLga2A20A9yBOKyjXAVCbMv1BzpiHWB0Ao4DN+bazixgBtbR+xFgWAH78wWygAGFxDQJWJtnmgItbc9nAC8U8b5yL98eOAbcDnQFEoCetvc60fYZVMv1ecQAwYCPbdpq4I5c2/4C+A6oAYQAe4Dbc8WeCUwBPAAfW7wngW6AN9bZzkGsL1t34AWspqSc7Y8FGts+63FAEtAo1/YzsM5e3LHOlI4AYpv/AzAHqA14AlfZpvcAzgKDbdsNBNra5i0EPgaqA/WBDcDdzv7bNA87/4edHYB5VL4HlxL/k8DLwHXAz7akdTHx57Pe28BbuV7XAv4EtgIfF7K/QNt22xayTFkl/nPAGWC/Lbm6AR8Cz+dZdneuBHkImJxn/sXEb0u26UD7XPPvBlbniv3PPOvPAD7N9XoKsDPX645AYiHvJQYYkWv7+3LN87W914ZAIyAbqJ3PNj7O/fvKNb0BkIbtS842bTy5vojMo2I/TFuiURpfAmuAZuRp5gEQkZ7AK1jNJV5ANWBuznxVTRSRucDDwOhc6/0T+Kft5Ve2+dlYSWqXI95ILmGqui/3BBFpCkwUkSm5JnthHWHniC1kmwFYR9KHc007jPWFVtj6x3M9T8nntV+uGCdgfU4htkl+tv3mOJbzRFWTRSRnmTrAaVU9k8/+g4Gl+Uxvans/R23bAesLsrDPwKhATBu/UWKqehir+eF6YEE+i3wNfA8Eq6o/8BFwMVOISBdgMjAbuNjzRVVfUlU/2+MeVU0G1pPryyEfSeS6uCwiDfOGW4y3llcs8KKq1sr18FXV2XZu/yRWU0vTXNOaAPFlEZ/ti+lT4AGgrqrWAraR67MuRCxQR0RqFTCvRQHT04CAXJ9HTVXtUJL4jfJnEr9RWrcDA1U1KZ95NbCOJlNFpAdWLxkARMQb62j+n8BtQKCI3FfIfh4DJonIoyJS17aNziLyjW3+ZqCDiHSxbfuZPOsfB5oX/+0BVlK9R0R6iqW67cJ1DXtWVtUs4FvgRRGpYUvUD2O9/7JQHeuL4wSArYtrqJ2xHQWWAR+ISG0R8RSRfrbZ04DbRGSQiLiJSKCItLWt8z/gDRGpaZvXQkSuKqP3YziYSfxGqajqflWNLGD2fcBzInIeeAor+eV4GYhV1Q9VNQ24BXhBRFoVsJ/fgIG2xwEROQ18gq0pQlX3AM8By4G9wNo8m5gGtLf1yllUzPcYiXVh9D2s9v99WO3mxTEF66zkgC22r4HpxdxGQfHtAN7AOis6jtX+v64Ym7gV64xkF9ZF7Ids292A9aX8FtZF3l+4dNYyAau5awfWZzIPqynOqARyruobhmEYLsIc8RuGYbgYk/gNwzBcjEn8hmEYLsYkfsMwDBdTKW7gCggI0JCQEGeHYRiGUalERUWdVNXL6l9VisQfEhJCZGRBPQYNwzCM/IjI4fymm6YewzAMF2MSv2EYhosxid8wDMPFmMRvGIbhYkziNwzDcDEm8RuGYbgYhyV+EQkWkVW28T+3i8hU2/SxttfZIhLuqP0bhmEY+XNkP/5M4O+qGm2rWx4lIj9jDRAxCmtYN8fLygT3SnG7gmEYRrlw2BG/qh5V1Wjb8/PATiBQVXeq6m5H7Te3XRuXk/Raezi+ozx2ZxiGUSmUSxu/iIQAXYE/irHOXSISKSKRJ06cKNF+v4/zJSvlPOn/e7pE6xuGYVRFDk/8IuIHzAceUtVz9q6nqp+oariqhterd1mpCbvc0DuUDzNvwGv//+BQcQYkMgzDqLocmvhFxBMr6c9S1fwG43aotg1rEtkwghNSF/35KTCjjRmGYTi0V49gjXO6U1XfdNR+inJjj5a8lj4KiY+Eg784KwzDMIwKw5FH/FdiDeI8UERibI/rRWSkiMQBvYEfROQnB8bA8M6N+cGtP581ewuaXeXIXRmGYVQKDuvnqKprASlg9kJH7Tevmt6eXNsxiLe3e3JTRha+7grunuW1e8MwjArHJe7cHRcezIW0THb88AG8Fw7pSc4OyTAMw2lcIvH3aFaHZgHVWRRbHc4cgvUfODskwzAMp3GJxC8ijA0P4qsjDUlqdh2s+y8knXR2WIZhGE7hEokfYExYEO5uwqwat0FGEqx5zdkhGYZhOIXLJP76Nb0Z0KYen+70ILvLrRD5OSSfdnZYhmEY5c5lEj9ARHgwJ86n8WvQXXDnSvCt4+yQDMMwyp1LJf4BbesT4FeNL7elQsNQa2JmunODMgzDKGculfg93d0Y0y2IVbsTSDiXCksfhdnjnB2WYRhGuXKpxA8QER5EVrYyPzoeajeD/Suth2EYhotwucTfvJ4fPULqMDcyFg2fDLWawM9PQ3a2s0MzDMMoFy6X+AEiugdz4GQSG+OSYeC/4dgW2Dbf2WEZhmGUC5dM/Nd3bIhfNQ++2fgnhI6Bhh1h7VumbLNhGC7BJRO/r5cHwzs3ZunWo5xLz4KRn8CERSAF1ZQzDMOoOlwy8QOM6x5MakY2izcfgQbtwa++dcRvuncahlHFuWzi7xzkT5sGNfh2Y6w1ITMNpl8Ha151bmCGYRgO5rKJX0SI6B7M5riz7Dp2DjyqgX8QrH8fzh9zdniGYRgO47KJH2Bk10A83YU5OUf9A5+ErAxY/bJzAzMMw3Agl078dap7cU37hizcFE9aZhbUaQbhkyH6Szixx9nhGUbVZHrPOZ1LJ36w+vQnJmfw847j1oSrHgNPX1j/nnMDM4yq5MIJOPCL9XzrXJg+BPavMl8CTuKwxC8iwSKySkR2iMh2EZlqm15HRH4Wkb22n7UdFYM9+rQMoLG/96XmnuoBcOsCGGIu8hpGmUhPgq8jYM4tkHIGxA3OHIQvb4Rp18Den80XQDlz5BF/JvB3VW0P9ALuF5H2wOPAClVtBaywvXYadzdhTHgwa/edJO5MsjUxuAd4eltdO80fpGGUXFYmzJ0ER2Ng5MfgUxs6joEHY2DoG3D+KMwaAwvudHKgrsVhiV9Vj6pqtO35eWAnEAiMAGbaFpsJ3OioGOw1tlsQAPOi4i5NPLUf3u0Ge35yUlSGUcmpwg9/g73/g+tfh7bXX5rn6Q3d74Ap0TD8HQgdbU1POw87F5vaWQ5WLm38IhICdAX+ABqo6lHbrGNAgwLWuUtEIkUk8sSJEw6NL7iOL31aBjA3Mo7sbNsRfq0m4OEFy5+B7CyH7t8wqqTdyyD6C+j7d+h+e/7LeHhBt4nQZoj1evM3VpPQR1da9bPM/55DODzxi4gfMB94SFXP5Z6nqgrk25aiqp+oariqhterV8/RYRIRHkx8Ygrr9tsGYXf3hEFPwYmdsHm2w/dvGFVOmyEQ8aVVCNFe4ZNh9DQr4c+bDB/0gs1zTJNrGSsy8YvIWHumFbCuJ1bSn6WqC2yTj4tII9v8RkCC/eE6zjUdGlDL15Nvci7yArS7AQLDYdVLkJHivOAMozI5sBpO7rNqX7W/oXg1sNzcrWsA9/0OY2eAmyfEfHVpG+YLoEzYc8T/hJ3T/kJEBJgG7FTVN3PN+h6YaHs+EfjOjhgcrpqHOzd2CeTn7cc5k2Sr1yMCg5+Dc/GwdZ5zAzSMyuDIJph9Eyx9pHTbcXODDiPhnrUw1nZJ8Gw8vN8DomaYmlqlVGDiF5EhIvIuECgi7+R6zMDqsVOUK4FbgYEiEmN7XA+8AgwWkb3A1bbXFcK47sGkZ2WzcFP8pYkhV8KkpdD1FucFZhiVwZlDMCsCfOvCyI/KZptubuBbx3qecgaq1YDFU+GdrrDhU8hILZv9uBjRAk6dRKQz0AV4Dngq16zzwCpVPePw6GzCw8M1MjKyXPZ1w3trSc/MZtnUvkjeU9TMNKumj2EYf5V0CqZfA0kn4fb/Qb02jtmPqjVU6i+vQuzv4B8M928AL1/H7K+SE5EoVQ3PO92joBVUdTOwWURmqao9R/hVQkR4ME8u2saWuLN0Dq51acbuH+G7++GuVVaPH8MwLln9MiTGwsTvHZf0wWp+bTkIWgyEQ79aTUs5SX/7Imh5NVTzc9z+qwh72vj3isiBvA+HR+YkN3RpjLenG3MiY/86o2Go1cd41UvOCcwwKrLBz8GE76BJr/LZnwg06wdXTrVen9oPcyfC2x3h1zcg9Vzh67s4exJ/ONDd9ugLvAN85cignKmmtyfXhzZiccwRUtJz9SH2D4Ked1v9jI9tc16AhlFRqFrt7KlnraPupr2dF0vdFnD7zxDYDVY8Z30BrP6PdbBmXKbIxK+qp3I94lX1bWCo40NznojuwZxPy2Tp1qN/ndH3YfCuCSuedU5ghlGRrHvb6r0T/aWzI7EE94Bb5sGdq6DpFfDbu1aZdTDdQPOwpx9/WK5HuIjcQyHXBqqCns3qEFLX9/LmHp/a1l2Ie/8HCTudE5xhVARbvrXuag8dDb3uc3Y0fxUYBuNnw9QYq0eQKswYZsWbdNLZ0VUI9iTwN3I9zwQOAREOiaaCEBHGhgfz2k+7OXgyiWYB1S/N7HE3hPSB+u2cF6BhONOB1bDoPgjpCzd+aHW5rIiqB1g/0y9AjQaw9m3442Pr7uArHrSmuagCu3NWJOXZnTPH8XOp9H55BXdf1YJ/XNc2/4UyUq1iU4bhKrKz4IPe1h22ty0Dn1rOjsh+J/ZYF363fgvuXjDhe2jS09lROVRB3TkLu4HrYRG5rLKSiNwuIg+VcXwVToOa3gxsW5/5UXFkZuVTKXDN6/BRn0ttiIbhCtzc4daFcPO8ypX0Aeq1hlEfwwOREH47NO5iTT/4q9UV1YUUdo52M/BFPtO/BCY7JpyKJSI8mITzaazenU910IYd4dRe6/Zxw6jqUs5YR8vZWeAfaD0qq7ot4LqXrJsxs7Ph+wesO4G/nwKnDzo7unJRWOL3UNXLDmdVNR0oRtWlymtA2/oE+FW7/CIvQKtroOmV8Mt/IO1C+QdnGOUlI9Wqv7P6larXqcHNDSYugW6TrK7a73aDhffC6Sp7qxJQeOJ3E5HLrn7kN62q8nR3Y3S3QFbuSiDhXJ6aIDkF3JJOmPF5jaorOxsW3g1//mZdyG0Y6uyIyl6tYBj6OkzdYt2rs30BnNxrzasE10BLorDE/xrwg4hcJSI1bI/+wBLg9fIIriKICA8mK1uZHx1/+cygcKt084ZPTLEoo2r635OwYxFc84JVLrkqq9kIrnsZ/rbdOqMHWPk8fDuxyt20WVitni9E5ARWkbZQrAFTtgNPqeqycorP6VrU86N7SG3mRsZyz1XNLy/cdq2thIPp3WNUNacPQOQ06Hkv9H7A2dGUn5xuoACevrBvhfXl13YYXPUYNOrstNDKiunOaYe5kbE8Om8L397dmx7N6uS/kCpkpoKnT/kGZxiOlLATAtpU3L765SH5NPzxEfz+EaSdhUFPW3fxVwLF7s5pXDK0UyP8qnkwZ2MBXb5U4ZubYdG95RuYYTjCobWwaZb1vH471076YN39O+Cf8LetMOBJqwIoWD2A/vzDubGVkIv/Ru3j6+XB8M6NWLr1KOdT8+m3LwINOsD2hRAfVf4BGkZZSdgJ39wE6/5rjT9hXOLtD1c9Co06Wa9/e8cag2DmDdaXZSViEr+dIsKDScnIYvHmo/kvcMUUa+Shn5+usj0BjCru3BH4agx4+FjFzsygQ4W75gW45kXry3LGUPj8eji4xtlR2cWeIm1fioh/rtdNRWSFY8OqeLoE16J1A7/8+/SDVbXzqn9Yg0Psc7mPx6jsUs/CrLHWz5vnmsGG7OFVHa54AB7aAkNetZp+dv1gzVOt0AeA9hzxrwX+EJHrReRO4GfgbYdGVQGJCBHhwWyOTWT3sQJqfHe7DWo1hfXvlm9whlFau3+EE7th3BeXmjIM+3j6WP3/p8ZY1wLAavr5dADsXlYhvwDsqcf/MXAH8B1W185+qrq4qPVEZLqIJIjItlzTOovIehHZKiKLRaRmaYIvb6PCgvB0l4Iv8np4wf99DeOq7Dg1RlXVeRw8sMEa0tAoGY9q1nUAsHr4JZ+G2f8HH/eFHd9bN8NVEPY09dwKTAcmADOApbaB2IsyA7guz7TPgMdVtSOwEHi0OME6W53qXlzTviELN8WRlpmV/0INQ6FaDcjKhMz08g3QMIprzWsQu8F6Xqe5c2OpSloNhilR1t3O6cnw7a3w5Y3Ojuoie5p6RgN9VHW2qj4B3APMLGolVV0DnM4zuTWQc/XjZ9u2K5WI7sGcSc5g+Y6EghdKPg0f9IKNn5ZfYIZRXBs+hZUvwLb5zo6kanL3hC43wf0bYNRn0Hm8NT0r0xoYPivTaaHZ09Rzo6om5Hq9AehRwv1tB0bYno8FggtaUETuEpFIEYk8cSKf6phO0qdlAI39vQu+yAtWv99awdbRVEpiucVmGHbbuRiWPgptrr9097nhGO4e0GksdLEl/j3LrIHh3+8Om75ySml3e5p6vEXkfhH5wNZuPx34qIT7mwzcJyJRQA2gwLYQVf1EVcNVNbxevXol3F3Zc3cTxoQH8+veE8QnphS84NXPWKVs1/233GIzDLv8+QfMv8OqNTV6mlVj3yg/bYbCuFng5Qff3Q/vhkHk5+V6BmBPU8+XQEPgWuAXIAgo0dD1qrpLVa9R1W7AbGB/SbbjbGO7BQFWKYcCNeoMHSPg9w+t/tGGUVFs+gJqBsL4OeDl6+xoXI+bG7QbBnevgZu+her1rDyRUwesHHoB2ZP4W6rqv4EkVZ0JDAVKNF6ZiNS3/XQDnqTkZw5OFVzHlytbBDA3Mo7s7EJ+SQP/BZoFkdPLLzjDKMrwd6xhE6vXdXYkrk0EWl8Ld6yAST9YZ15p5+HDK60vgvRkh+3ansSf0wCVKCKhgD9Qv6iVRGQ2sB5oIyJxtmEcx4vIHmAXcAT4vGRhO19E92DiE1NYt/9kwQvVDoHJP0L/J8otLsPIV9p5WHAXnI23EowLDzRe4YiAn605O+kk+NSGHx+H/3Zy2EFjgWWZc/lERGoD/wa+B/yAp4paSVXHFzCrSjR6X9O+Af4+nszZGEvfVoVcgwjsZv1MTzan1YZzZGVYNeUPrIZOEZV72MSqrk4zuO0HOLQO1rwKSaccspsiE7+qfmZ7+gtgOvraeHu6M7JrIF//8SdnktKpXd2r4IXjImHWGBj/DTTpVX5BGoYqLJ4K+1fADe9eqixpVGwhV0LIdw676avAxC8ihRacVtU3yz6cyiUiPJgZvx1iUUw8t13ZrOAF67cD92rw81Mw+adLF3EMw9FWvQQxs+CqxyFsgrOjMYrLQSWxC9vq68AtQF2s5p0aeR4ur33jmnQM9GfOxlgKHdDGqzr0fxxi/7hUxMkwHC092fp763qr9fdnGDaFJf6uwP+wevE0BdYBz6nqs6r6bHkEVxlEdA9m17HzbI0/W/iCXW+Fuq1gxbNOvWPPcCFevjB5GQx7y5xlGn9RYOJX1c2q+riqdgGmYd1xu0NEbiiv4CqDGzo3ppqHW8GF23K4e8DVT8PJPbB/ZfkEZ7imuCirB09GilU0zN3T2REZFYw9d+7Wwzr67wjEAYUUqXE9/j6eXN+xEd/HHCElvYDCbTnaDoM7V0Lra8onOMP1nNoPX4+FP3+HtAvOjsaooApM/CIyWUR+BOYCAkSo6mBV/b3coqskxnUP5nxaJsu2FTA6Vw6RXN07kxwfmOFaLpyAr0ZbPXluWXCpb7hh5FHYEf9nQGOs8gzXAp+JyPc5j3KJrpLo2awOIXV9i27uyRH9JbzdyariaRhlIT0Jvo6A88esMgABLZ0dkVGBFdaPf0C5RVHJiQhjw4N57afdHDyZRLOA6oWvEBQOKadhzetwnamMaJSBM4fgbCyMmQ7B3Z0djVHBFXZx95fCHuUZZGUwplsQbgLfFla4LUf9dlad7o2fwpnDjg/OqPoadIAHY6Dt9c6OxKgEHHN3gAtqUNObAW3qMz8qjswsO+626/9PEDdY9aLjgzOqrjWvwepXrHb9an7OjsaoJEziL0MR3YNJOJ/G6t12DBzjHwg974HtC612WcMork2zrBG0Th90diRGJWMSfxka2LY+AX7VCh+dK7c+f4P7focaDR0bmFH17FsOix+E5v2tGjzmBi2jGAqr1bMYKLAOgaqaG7ny8HR3Y3RYIJ+tPUjC+VTq1/AufAWfWtYDrD7X5lTdsMeRGJgzAeq1g4gvwaOQAoGGkY/CevW8Xm5RVCFjw4P5eM0BFkTHc89VLexb6Ye/Q3wU3LHSYUWZjCrk5F7wqw83zwXvms6OxqiECkz8pudOybSs70d409p8uzGWu/s1R+w5BQ8Mh42fwY6FEDra8UEalZOq1aTTaSy0vwE8qjk7IqOSsqdkQysRmSciO0TkQM6jPIKrrCK6B3PgZBKRh8/Yt0KnCGgQCiueh8wCx583XFlGCnwxAnYutl6bpG+Ugj3tCp8DHwKZWDd1fQF85cigKruhHRtR3cvd/jt53dzh6mfgzEGImuHI0IzKKDsL5t8BB9dYzw2jlOxJ/D6qugIQVT2sqs9glWo2ClC9mgfDOzfmhy1HOZ+aUfQKYI2MFNIXNnzisFF3jEpIFZb9A3Ytgetehg43OjsiowqwJ/GniYgbsFdEHhCRkVgDsxRKRKaLSIKIbMs1rYuI/C4iMSISKSI9ShF7hTauezApGVks2VJE4bYcIla3vDt+Nhd4jUvWvW3d4d37Aeh1r7OjMaoIezLMVMAXeBDoBtwKTLRjvRnAdXmmvQo8a6vx/5TtdZXUJbgWrRv42d/cA9ZAyz61rSN+U73TUIXzx60L/oOfd3Y0RhVSZOJX1Y2qekFV41T1NlUdZU9pZlVdA+QtP6lATv8zf+BIsSOuJESEiPBgYmIT2XP8vP0rZmXCZ4Pgf/92XHBGxZedbZ0FXvcyjPzEnAW6IFUl+s8zpGeWfdOvPb16wkVkoYhEi8iWnEcJ9/cQ8JqIxGLdJ/BECbdTKYwKC8LTXYp31O/uYdXsj5oBJ/c5LDajAju2DT7sDQk7reTvXtjtNkZVE3cmmXdX7GXgG78w6oPf+GWPHSVgismev6hZwKPAVqC0Xz33An9T1fkiEoE1pOPV+S0oIncBdwE0adKklLt1jjrVvRjcvgELouN47Lo2VPNwt2/Fqx6DmK9h5fMQMdOxQRoVS2IszBoDCFSr4exojHJyIS2TZVuPMj86jt8PWA0lvZrX4d7+Lejdom6Z78+exH9CVctq4JWJWNcMwBrZ67OCFlTVT4BPAMLDwwssHVHRRYQHs3TrMZbvSGBop0b2reRXH66YAr+8Yo2fGtTNsUEaFUPKGSvppyfB5B/BP8jZERkOlJWt/Lb/JAui4/lx2zFSMrJoFlCdvw9uzY1dAwmu4+uwfduT+J8Wkc+AFUBazkRVXVCC/R0BrgJWAwOBvSXYRqXSt1U9Gvt7Mycy1v7ED3DFAxA5zerRYRJ/1ZeRCt/cbI2Ze+sCq76+USXtSzjPvKh4Fm2K59i5VGp6ezAyLJDRYUGENall393+pWRP4r8NaAt4cqmpR4FCE7+IzAb6AwEiEgc8DdwJ/FdEPIBUbE05VZm7mzCmWxDvrtpHfGIKgbV87FuxWg2Y8D0EtHJsgEbFkJUObh4w8iNo1s/Z0Rhl7HRSOos3H2F+dBxb4s7i7ib0b12Pfw9rz6B29fH2tLMZuIyIauGtKCKyW1XblFM8+QoPD9fIyEhnhlAqsaeT6fvqKv52dWumXl2CRJ6eBB7e1h2+RtWTlQHunlZPHtN7p8pIz8xm5a4E5kfHsWpXApnZSvtGNRndLYgbOjemXg3Hl90QkShVDc873Z4j/t9EpL2q7nBAXC4huI4vV7asy9yoWKYMbImbWzFO5c4cgmnXWCUdutzkqBANZ1n/gTUYzy3zTaXNKkBV2RJ3lvnRcXy/+QiJyRnUq1GN264MYVRYEO0aVYzfsT2JvxcQIyIHsdr4BVBV7eTQyKqYiPBgpn4Tw2/7T9GnVYD9K9ZqCjUDYeWL0GEUeBZR49+oPLYvhJ/+Ce2GgVd1Z0djlMLRsyks3BTP/Kg49p9IopqHG9d0aMiosED6tgzAw71incnZk/jz3n1rlMC1HRri7+PJnMjY4iV+ERj8LMwcbtXxufJBxwVplJ9D62DBXRDcE0Z9aprxKqHk9Ex+3HaMBdHxrNt/ElXoEVKHO/s25/pOjajp7ensEAtU2AhcNVX1HFCM206Ngnh7unNjl8bM3hhLYnI6tXyLMWpSs37QcjD8+gaE3WqVdTAqr4Rd8M14qB0C42eDp50X/A2ny85Wfj94ivlR8SzbdpTk9Cya1PFl6qBWjOwaSNO6lePMrbAj/q+BYUAUVi+e3A3TCjR3YFxVUkT3YGauP8yiTfFMurJZ8Va++mn4qC9smw/d73BMgEb5cPOwhk0c9Qn41nF2NIYd9p+4wMLoeBZuiic+MYUa1Ty4oXNjRncLIrxp7XLpglmWiuzVUxFU9l49uQ1/dy2Z2crSB/sU/4/l2FZrwJZK9kdm2GSkWL2zRC6NpmVUWInJ6SzecpT5UXHExCbiJtCvdT1GhQVxTfsG5d4FsyRK3KtHRFao6qCiphn2iegezL8XbWNb/Dk6BvkXb+WGHa2fqedMD5DKJjMdvh4HdZrD8LdN0q+gMrKyWb37BAui41ixM4H0rGzaNqzBv65vx4gujalfs2p0riisjd8bqxxzgIjU5lJTT00gsBxiq5Ju6NyYF5bsYE7kn3QM6lj8Dez5H8ydBHcshwbtyzw+wwFU4fsH4OAv0Hm8s6Mx8lBVth85x7yoOBZvPsKppHTqVvfill5NGd0tkPaNala6ppyiFHbEfzdWNc3GQHSu6eeA9xwYU5Xm7+PJ9R0b8d2mI/zr+vb4eBXzdDEo3GojXvEs3DTHMUEaZWvFs7BlDgz8N3Qxib+iOH4ulUWb4pkfHcee4xfwcndjcPsGjAoLpF/renhWsC6YZanAxK+q/8UqrzBFVd8tx5iqvIjwYBZusnoFjAorZiEu3zrQ5yErmRxaByFXOiRGo4xs/AzWvgXdboO+f3d2NC4vJT2L/+04xvzoeNbuPUG2QliTWrw4MpRhHRvj71txu2CWJXv68X8mIg8DfbB68/wKfKSqqQ6NrArr1bwOTev6MmdjbPETP1hD8G34FJY/Dbf/bNqLK7LaIRA6Bq5/3fyenCQ7W9l46DTzo+NYuvUYF9IyCazlw/0DWjIqLIhmAZWjC2ZZsifxz8Tqy59z1H8T8CUw1lFBVXU5o3O99tNuDp1MIqS4f3iePjDgCfj+QTiyCQLDHBOoUXKpZ8HbH1pebT2McnfoZBILNsWzcFMcsadTqO7lzvUdGzEqLIiezeoUr3RKFWNP4g9V1dxXEVeJiKnbU0qjw4J443+7+TYylseua1v8DXS+CYJ6QP0SrGs41sm9MP06uOYF06Zfzs6mZPDDlqMsiI4j8vAZRKBPywD+PrgN13RogK+XGc0M7Ev80SLSK2ecXRHpCVSNTvVO1NDfm/5t6jMvKo6HB7cufi0Pd49LST/n6NJwvvPH4atR1vMmPZ0bi4vIzMrm170nmRcdx887jpOemU3L+n7847q2jOwaSEP/qtEFsyzZk/i7YVXo/NP2ugmwW0S2Yoq1lUpEeDArdyXwy54TDGrXoGQbWfMabPgMHow2hb6cLe0CfB0BSSdh0hKrz77hMDuOnGNBdByLYo5w8kIatX09ualHE0aFBdIx0L/KdcEsS6ZImxMNalefAD8v5myMLXniD+kHK1+wyvte9WjZBmjYLzsL5k607q4ePxsCzahpjnDifBrfxcQzPzqenUfP4ekuDGxbn9FhQfRvUx8vj6rbBbMs2ZP4mwM548BtV9VVDozHpXi6uzEqLIjpaw9y4nxayQZmaNIT2g6Ddf+F8NugejEqfxplR9ygeX9oNxxaX+vsaKqU1Iwslu88zvyoONbsPUlWttI5uBbPjejA8E6NqV29GAUPDaDwO3cDsYZXTMUq1AYwVkT+A4xU1fhyiK/KiwgP5pM1B1gQHcfdV7Uo2UYGPQW7e1nNPkP+U7YBGkVLOml94V4xxdmRVBmqStThM8yPjmfJliOcT82kkb83d/drzqiwQFrWr+HsECu1wo743wM+VNUZuSeKyATgA2CEA+NyGS3r+9GtaW3mRMZyV7/mJWuXrNcGut566e7Qan5lH6iRv6gZ8L9/w+SfTAmNMhB7OpkF0fEs2BTH4VPJ+Hi6MyS0IaO7BdGreV3cXbgLZlkqLPG3V9WReSeq6hci8q+iNiwi07HKOieoaqht2hwgZ/zeWkCiqnYpbtBVzbjuwTw2bwtRh88QHlLCMr2DnoKBT5qkX15O7rWS/u8fQosBEFCCsZQNAM6nZrBs6zHmRcex4eBpRKB387pMGdiKIaENqV7NdMEsa4V9ovleJRERN8CeAjMzsM4avsiZoKrjcm3nDeCsXVFWcUM7NuLZ77czZ2NsyRN/Ttu+qtW906dWmcVn5PHtRNixyKqZ1OFGGP6ONVi6YbesbGXtvpMsiI7jp+3HSM3IpnlAdR69tg03dg0ksJYZnMaRCkv8S0TkU+AhVU0CEJHqwFvA0qI2rKprRCQkv3litWdEAAOLHXEVVL2aB8M7N+a7mCM8Nbw9NUozZNvs8ZCZChMWlVl8Lu/YNmsAnIFPWkMkNukFjbtAl5vBr76zo6tU9hw/z/yoOBbFxHP8XBr+Pp6M7RbMqLBAugTXMl0wy0lhif8x4GXgsIgctk1rglXC4Z+l3G9f4Liq7i3ldqqMiO7BfLMxliVbjjK+R5OSb6hZX2sA7/0roYX5Xi2xtAtWso+eCfFR4O4FoaOsMRF63evs6Cqd1IwsXvhhB1/9/icebkL/NvV5ZnggA9vVp5pHxR/QpKopcgQuEfEBWtpe7lfVZLs3bh3xL8lp4881/UNgn6q+Uci6dwF3ATRp0qTb4cOHC1q0SlBVrnlrDdWrebDo/lJU3MxMg/fCwbsW3PULuJl+zcV2fAdMGwzpF6BeWwibCJ3/zwyTWEL7Ei7wwNfR7Dp2njv6NOOe/i0I8CtB12Wj2Eo8ApeqpgBbyzAQD2AU1h3Bhe33E+ATsIZeLKv9V1Qiwrjuwbzww072HD9P6wYl7K7mUc3q2bPgTuuItZOppVeklETYOteqntn9jku9pDqMhOAepqpmKcyLiuPfi7bh4+XO57d1Z0Ab0zRWETjjcPBqYJeqxjlh3xXayK6BeLoLczbGlm5DoWOsJomoGWUSV5WkCofXw8J74I22sPQR2POTNc/NHYa8Yt0cZ5J+iVxIy+ThOTE8MncznYP9WTa1r0n6FYjD+kmJyGygP9bQjXHA06o6Dfg/YLaj9luZ1fWrxtXtGrBwUzz/uK5tyW8/d3ODiC+hRqOyDbAq+elf8Pv74FXDqqAZNtG6YGuU2vYjZ3ng600cPpXE365uzQMDW5r+9xVMYXfuFlrkXVWji5ifbz1aVZ1kV2QuKqJ7MMu2HWP5zuNc37EUibtOM+tnZppVR8bLt2wCrIyys+HQGoiaCVc9BvXbQcfR1s/QUaa4XRlRVb78/TAvLNlJ7eqefH1nL3o1r+vssIx8FHbEn3Ph1RsIBzZjDbjeCassc2/Hhuaa+rWqRyN/b+ZsjC1d4gdIOw8fXgGdxlldEV3N+WOw6SvY9CWcOWRd8A4dZSX8wG6mkFoZOpucwWPzN/PT9uMMaFOP18d2pq65gFthFTbm7gAAEVkAhKnqVtvrUOCZconOBbm7CWO6BfHeqn0cSUyhcWluZKlWAwLDYf371kXLGg3LLtCKLiMV3usOaeegaR8Y8KRVQM3T1GYva1GHz/Dg7E0knE/lyaHtmHxlM5ce3aoysKcRuU1O0gdQ1W1AO8eFZIztFoyq1SOi1Ab9G7LSYfUrpd9WRZYYC6tetm5gAyvBD/8vPBAFt/1g9W4ySb9MZWcrH6zeR8TH63Fzg7n3XMEdfZubpF8J2HNxd4uIfAZ8ZXt9M7DFcSEZTer6ckWLunwbGcsDA1qW7h+pTnMInwwbp0Hv+6tWTZmsDNi9DKK/gH3LrWktBlo3X1Xzs5p1DIc4cT6Nh7+N4de9JxnasREvj+5IzdLccW6UK3uO+G8DtgNTbY8dtmmGA43rHkzcmRTWHzhV+o31e8waoH3TV0UvWxnk3HS4bT58eysc3w79HoWpm+HWBaZQnYOt23eS69/5lQ0HT/PSyI68d1NXk/QrGXtu4EoVkY+Apaq6uxxiMoBrOzTE38eTORtjubJlKQdX8asHd66q3Ef7mWmwc7F1b0LboVbZhHbDrbGGWw62xiA2HCozK5u3l+/l/dX7aFHPjy9v70HbhjWdHZZRAkX+t4jIDcBrgBfQTES6AM+p6g0Ojs2leXu6c2OXxszeGEticjq1fEs5ylC91tbP5NPgU7vy3JiUsMuql7N5NqScgVpNwMt2RO9VHdoMcW58LuJIYgpTv9nExkNniAgP4pkbOuDrZb5sKyt7mnqeBnoAiQCqGgM0c1xIRo6I7sGkZ2azaFMZDXZ2ZBO83fHSHaoVVVbGpedLH4ENn0Kzq+DWhfDgZgi71XmxuaDlO45z/Tu/suPIOd4e14VXx3Q2Sb+Ss+e3l6GqZ/OUS63ytXMqgg6N/QkNrMmcyDgmXhFS+pK1DULBrwEsfwZaDbZKE1QkRzdbF2q3L4R710ONBjD0DfCpYzVXGeUqLTOL/yzbzfR1B+nQuCbv3RRGswBzs1tVYM8R/3YRuQlwF5FWIvIu8JuD4zJsxoUHs/PoObbFnyv9xtw9rZG6Tuy0mk4qgrQLEDkdPr4KPu4H0V9Cy6utMQXAKphmkn65O3QyiTEfrmf6uoNMuiKEBfddYZJ+FWJP4p8CdADSgK+xRs2a6sigjEtu6BJINQ835kT+WTYbbD/CumN11UuQkVI22ywuVUi1fZGlnIElD1vNO0NehUd2w6hPoHZT58Rm8F1MPMPeXcufp5P5+NZuPHNDB1Mzv4qxJ/EPVdV/qWp32+NJwFzYLSf+Pp4MCW3IdzFHSM3IKv0GRWDwc3AuHvatKP32iiP5NPz+kVVGYp6tR3CtYLh/A9y7DnrebV14NpwiJT2Lx+dvYeo3MbRpWIOlU/tybQcXutvbhdiT+J+wc5rhIBHdgzmfmsmybUfLZoMhfeCBSGg3rGy2V5S4KJh/p1X++Md/WGMGtMt17FCvdeXpZVRF7Tl+nhHvr2VOZCz3D2jBN3f1MuPeVmGFVeccAlwPBIrIO7lm1QQyHR2YcUmvZnVpUseXORtjGdk1qGw2mtOnP/m0Y0aWunACvGtaSf7gL7DnR6s3TthEaNSp7PdnlIiq8s3GWJ5dvB2/ap58MbkHfVuZaypVXWFH/EewqnCmAlG5Ht8D1zo+NCOHm5sQER7E7wdOc/hUUtlteNNX8FYHq85NWcjOtpqPvp0Ab7aDHd9b03vcCX/fbfXQMUm/wjifmsGU2Zt4YsFWwpvWYenUPibpu4jCqnNuBjaLyNeqmlHQckb5GNMtmDd/3sO3kbE8em3bstlos6usWv2rXoKRH5Z8O5np8Nt/ra6YiX9a7fQ97oRA25AO1Uo4jKThMFviEnng603EJ6bw6LVtuPeqFqa4mguxp40/RETmicgOETmQ83B4ZMZfNPT35qrW9ZgXFUdmVnbZbLRWsHVBdfNsq95NcWRlXlrH3RO2LYDaITB6mnV0f93LULdF2cRplBlVZdrag4z+8Dcys7KZc1cv7i9tIUCj0rEn8X8OfIjVrj8A+IJLlTqNcjSuezDHz6WxZu+Jstton79ZbfHLn7Fv+TOHYeUL8HYoTLsW0pOsC7N3rICJi6HjGKtd36hwziSlc8fMSJ5fsoP+beqzdGpfwkMccH3HqPDsuXPXR1VXiIio6mHgGRGJAp5ycGxGHgPbNiDAz4s5G2MZ2LZB2WzUtw70eRhWv2wl9YL6zx+JgRXPwf6V1uuWV0O3ieBuS/KuPLRjJbDh4GkenL2J00npPDO8fdncCW5UWvYk/jQRcQP2isgDQDxg6t46gZeHG6PCgpi+9iAnzqdRr0YZHVn3vNs6UvfP02Po5D5r4PY6za2j+hO74Kp/QNdbrGYio8LLylY+WLWPt5bvoUkdXxbcdwWhgf7ODstwMnuaeqYCvsCDQDfgVmBiUSuJyHQRSRCRbXmmTxGRXSKyXUReLUnQriwiPJjMbGXhpjIYnSuHp4+V9FXhbDxs+RY+HwrvdYM1r1vLNOoMD22FAU+YpF9JJJxL5dZpf/DGz3sY3rkxSx7sa5K+AdhXj3+j7ekFijcAywzgPaxrAgCIyABgBNBZVdNEpH4xtmcALev70a1pbb7ZGMudfZuX7en6D3+HyGnW89ohVl2fLjdfml/RiroZBfplzwkenhNDcnoWr47pxNhuQaZpx7iosBu4FlNIFc6i6vGr6hoRCckz+V7gFVVNsy2TYH+oRo5x4cE8Nn8LUYfPlO3FuS43QXYGhI6GkH5WM49RqWRkZfPG//bw0S/7adOgBu/d1JVWDUx3WuOvCjvif90B+2sN9BWRF7FuDHsk1xnFX4jIXcBdAE2aNHFAKJXX0E6NeHbxduZsjC3bxB8Ubj2MSinuTDIPzt5E9J+JjO/RhKeHt8fb05ylGZcr7AauXxy0vzpAL6A78K2INFfVy84sVPUT4BOA8PBwU/8/l+rVPBjWqTGLtxzh6Rs64FfNDIrh6n7cdpTH5m1BFd67qSvDOjV2dkhGBVbkubyIHMx941Ypb+CKAxaoZQOQDZRyQFnXFNE9mOT0LJZsPuLsUAwnSs3I4qnvtnHPV9E0C6jODw/2NUnfKJI9h4q5z/29gbFYR+0lsQjrJrBVItIaaxzfkyXclksLa1KLlvX9mBMZy//1ME1hrujAiQvc//Umdh49x519m/HotW3x8jDXZYyiFflXoqqncj3iVfVtYGhR64nIbGA90EZE4kTkdmA60NzWxfMbYGJ+zTxG0USEceHBbPozkb3Hzzs7HKOcLYiOY9i7azl2NoXpk8L519D2JukbdivyiF9EwnK9dMM6A7CnG+j4AmbdYl9oRlFGhgXynx93MWdjLE8Oa+/scIxykJSWyVPfbWd+dBw9mtXhv//XhUb+pm6+UTz2NPW8ket5JnAQiHBMOEZxBPhV4+p2DViwKZ7HrjOn+VXdjiPneGB2NAdPJjF1UCseHNQKd1NczSgBe47cB5RHIEbJjOsezI/bj7Fi53GGdGzk7HAMB1BVvvrjT55fsoNaPp7MuqMnV7QwfSKMkrOnV89LIlIr1+vaIvKCQ6My7NavdT0a1vRmTmQZDaZiVChnUzK4b1Y0/160jd7N67J0al+T9I1Ss6dtYIiqJua8UNUzWEMyGhWAu5swNjyINXtOcPRsirPDMcrQpj/PMPSdX/l5x3H+eX1bPp/UnQA/U/LaKD17Er+7iFz8axMRH8D89VUgY7sFk60wL7IMC7cZTpOdrXz8y37GfrQegLn39OaufmaELKPs2HNxdxawQkQ+t72+DZjpuJCM4mpS15crWtRlTmSsGU2pkjt1IY2Hv93ML3tOMCS0Ia+M7oS/j6ezwzKqGHsu7v5HRLYAg2yTnlfVnxwbllFc47oHM/WbGNYfOMWVLU0bcGX02/6TPPRNDIkpGTx/Yyi39GxiKmoaDmFXkRdVXQYsc3AsRilc26EhNb09mLMx1iT+SiYzK5t3Vu7j3ZV7aRZQnRm39aB945rODsuowgory3ye/MsyC6Cqav4yKxBvT3du7BrINxtjOZucgb+vaR6oDI6eTWHqNzFsOHia0WFBPDeiA9VN0T3DwQq8uKuqNVS1Zj6PGibpV0wR4cGkZ2azKCbe2aEYdli56zjX//dXtsWf5c2IzrwR0dkkfaNc2H2rp4jUF5EmOQ9HBmWUTGigPx0a12TORtOnvyJLz8zmhSU7mDwjkob+Piye0odRYUFFr2gYZcSeG7huEJG9WKUafgEOYdr7K6xx3YPZcfQc2+LPOjsUIx9/nkpm7Ee/8dnag0zs3ZSF911Bi3p+zg7LcDH2HPE/jzVwyh5VbYbVu+d3h0ZllNiIzoF4ebiZo/4KaPHmIwx951cOnkzio1vCeHZEqBkhy3AKexJ/hqqeAtxExE1VV/HXGv1GBeLv68mQ0IYsioknNSPL2eEYWIOlPLFgK1Nmb6JlAz9+eLAv14WaukqG89iT+BNFxA9YA8wSkf8CSY4NyyiNceHBnE/N5Mdtx5wdisvbe/w8I95bx+wNf3LPVS349u7eBNfxdXZYhouzJ/GPAJKBvwE/AvuB4Y4MyiidXs3rElzHxzT3OJGq8u3GWIa/t5aTF9KYObkHjw9pi6e7KZ1tOF+Bf4Ui0lJErlTVJFXNVtVMVZ0JRAO1yi1Co9jc3KzRudYfOMXhU+bkrLydT83goTkxPDZ/C2FNarNsal+ual3P2WEZxkWFHX68DZzLZ/pZ2zyjAhvTLRg3gbmmcFu52hZ/luHvrmXx5iP8fXBrvry9J/Vrejs7LMP4i8ISfwNV3Zp3om1aiMMiMspEQ39vrmpdj3lRcWRmZTs7nCpPVfl83UFGffAbaZnZfHNXb6aYEbKMCqqwxF+rkHlFDvIpItNFJME2sHrOtGdEJF5EYmwPU9ffgcZ1D+bYuVTW7D3h7FCqtMTkdO76MopnF++gX+sAlj7Ylx7N6jg7LMMoUGGJP1JE7sw7UUTuAKLs2PYM4Lp8pr+lql1sj6X2hWmUxMC2Dahb3ctc5HWgjYdOc/1/f2X17gT+Paw9n04Ip3Z1L2eHZRiFKqwwyEPAQhG5mUuJPhzwAkYWtWFVXSMiIaUN0Cg5Lw83RoUF8vm6Q5w4n0a9Gmb8nLKQmZXNr3tPMj86jmXbjhFU24f5915Bp6Bazg7NMOxSYOJX1ePAFSIyAAi1Tf5BVVeWcp8PiMgEIBL4u20ox8uIyF3AXQBNmpjSQCU1rnswn/56kIWb4rirXwtnh1Op7Tx6jvlRcSyKOcLJC2nU9vVk0hUhPHR1K2p4V8xqqBkZGcTFxZGamursUAwH8vb2JigoCE9P+/4ORTW/ystlw3bEv0RVQ22vGwAnsco9Pw80UtXJRW0nPDxcIyMjHRZnVTfqg3WcTclg+cNXmYE9iunE+TS+i4lnfnQ8O4+ew9NdGNS2AaPCAunfpj5eHhW7X/7BgwepUaMGdevWNb/7KkpVOXXqFOfPn6dZs2Z/mSciUap6WaWFcq0BazuLyAnoU2BJee7fVY3rHsw/5m8l+s8zdGtqLjoWJTUji+U7j7MgOp5f9pwgK1vpHFyL50Z0YHinxpWqDT81NZWQkBCT9KswEaFu3bqcOGF/J45yTfwi0khVj9pejgS2Fba8UTaGdmrMs4t3MGdjrEn8BVBVog6fYX50PEu2HOF8aiaN/L25u19zRoUF0rJ+DWeHWGIm6Vd9xf0dOyzxi8hsoD8QICJxwNNAfxHpgtXUcwi421H7Ny7xq+bBsE6NWLLlKE8N74CfGezjotjTySyIjmfBpjgOn0rGx9OdIaENGd0tiF7N65p++EaV5LAGSlUdr6qNVNVTVYNUdZqq3qqqHVW1k6rekOvo33Cwcd2DSU7P4octR5wditOdT83g242xRHy8nr6vruLtFXsIrOXD62M7E/nk1bw5rgtXtgwwSb+MuLu706VLFzp06EDnzp154403yM6+dFPhhg0b6NevH23atKFr167ccccdJCcn57utRYsW0alTJ9q2bUtoaCjz5s37y/zMzEzq1avH448//pfpZ8+eZcKECbRs2ZIWLVowYcIEzp7965gVDz30EIGBgX+JrcpS1Qr/6Natmxqlk52drQNfX6Uj31/r7FCcIjMrW1fvTtAHZ0drmyeXatN/LNEBr63S91bu1bgzyc4Oz2F27Njh7BC0evXqF58fP35cBw0apE899ZSqqh47dkybNGmiv/3228Vl5s6dq8eOHbtsOzExMdqiRQs9cOCAqqoeOHBAmzdvrpGRkReXWbp0qV5xxRXavHlzzc7Ovjh99OjR+vTTT198/dRTT+mYMWMuvs7KytImTZpoz549deXKlaV/006Q3+8aiNR8cqpDe/WUFdOrp2x8smY/Ly3dxfKH+1XqNuvi2HP8vK0LZjzHz6Xh7+PJDZ0bMyoskC7Btap8+/fOnTtp164dAM8u3s6OI/mV3yq59o1r8vTwDoUu4+fnx4ULFy6+PnDgAN27d+fkyZM8/fTTADz33HNF7uvWW29lwIABTJ58qSPgtGnTWLFiBV9//TUAEyZMYNiwYXz44Ye8+OKLXHHFFezbt4/Bgwezb98+3N2tgW+ysrJo2bIly5cvp0WLFqxcuZLXX3+dcePGsW7dOj755JNifxbOlvt3naOgXj0Vuy+aUaZGhQXh4SZV/k7eUxfS+HzdQYa9+yvXvLWGaWsP0jGwFh/eHMaGfw3i+RtD6dqkdpVP+hVV8+bNycrKIiEhgW3bttGtWze71tu+fftly4aHh7Njxw7A6sG0fPlyhg8fzvjx45k9ezYAO3bsoEuXLheTPlxqftq+fTsAs2fPZvz48YwcOZIffviBjIyMsnirFZa5yudCAvyqcXW7BiyIjufRa9tW+D7oxZGWmcXKnQnMj45n9e4EMrOVjoH+PD28PcM7NybAz9y1XNSReWW3ZMkSBgwYgI+PD6NHj+b555/n7bffLnK99PR0li5dyptvvkmNGjXo2bMnP/30E8OGDXN80E5iEr+LGdc9mB+3H2PFzuMM6Vi5h/9TVWJiE5kfHcfizUc5m5JB/RrVuL1PM0aFBdGmoWs0Z1U2Bw4cwN3dnfr169OhQweioqIYMWLEZctde+21HD9+nPDwcD777DPat29PVFQUnTt3vrhMVFQU4eFWS8bs2bNZu3YtISEhAJw6dYqVK1fSvn17YmJiyM7Oxs3NOtjJzs4mJiaG9u3b89NPP5GYmEjHjh0BSE5OxsfHp0onfqdfuLXnYS7ulp3MrGzt+eJynTj9D2eHUmJxZ5L1vZV7dcDrq7TpP5ZomyeX6oOzo3X17gTNzMouegMupKJd3E1ISNDBgwdfdnH3999/v7jM/Pnz8724u2nTJm3ZsqUePHhQVVUPHjyoHTt21F27dunZs2e1Xr16mpqaenH56dOn62233aaqqiNHjtRnn3324rxnn31WR40apaqq48eP16+//vrivAsXLmi9evU0KSmpDN59+SnOxV1zxO9i3N2EMd2C+GD1Po6eTaGRf5EVtiuEpLRMlm07xoLoONYfOIUq9GhWh3v6tWBIx4YVtlaOASkpKXTp0oWMjAw8PDy49dZbefjhhwFo0KAB33zzDY888ggJCQm4ubnRr18/rrvu8sK+Xbp04T//+Q/Dhw8nLS2NQ4cOsWrVKtq0acPMmTMZOHAg1apdatIbMWIEjz32GGlpaUybNo0pU6bQooVVr6p3795MmzaN5ORkfvzxRz766KOL61WvXp0+ffqwePFixo0b5+BPxzlMrx4X9OepZPq9toq/D27NlEGtnB1OgbKyld8PnGJ+lFUFMyUji6Z1fRnVNYhRYYFm0HI75NfTo6p4/PHH+eOPP/jpp5/w8qo8ZTQcpTi9eswRvwtqUteX3s3r8m1ULPcPaIlbBbtRaV/CBeZHx7FoUzxHz6ZSw9uDG7sGMjoskG5NTW8cw/LKK684O4RKyyR+FzWuezAPzYnh9wOnuKJlgLPD4UxSOou3HGF+dDybYxNxdxP6tQrgX0PbcXW7Bnh7uhe9EcMw7GISv4u6LrQhNb7zYE5krNMSf3pmNqt3JzA/Oo6VuxLIyFLaNarJk0PbcUOXxtSvYQYpNwxHMInfRXl7unNjl0DmRMbyXHIG/r7lc3FUVdkaf5YF0fF8v/kIp5PSCfCrxsTeIYwKC6J945rlEodhuDKT+F3YuO7BfPn7Yb7bHM+E3iEO3dfRsyks2nSEBdFx7E24gJeHG4PbN2BMWBB9WwXg4V51biYzjIrOJH4XFhroT/tGNZmzMdYhiT85PZOfth9jQXQ8a/edRBXCm9bmpZEdGdqpEf4+pgumYTiDOcxyceO6B7P9yDm2xZ8temE7ZGcr6/ef4pG5m+n+wnL+NmczB08mMWVgK1Y/0p95917BTT2bmKTvQl588UU6dOhAp06d6NKlC3/88QcA/fv3p02bNhfLLD/wwAMkJiZeXC+nnk5oaChjx44tsFSzM/Xv35/y7Gp+6NAhQkNDi16wCOaI38Xd2CWQF5fu5NvIWEID/Uu8nYMnk1gQHceC6HjiE1Pwq+bB0E6NGB0WRPeQOhWuy6hRPtavX8+SJUuIjo6mWrVqnDx5kvT09IvzZ82aRXh4OOnp6TzxxBOMGDGCX375BQAfHx9iYmIAuPnmm/noo48u3vjlLFlZWX8p9lZZmcTv4vx9PRkS2pBFm+L55/XtitVt8mxyBou3WO320X8m4ibQp1U9HruuDde0b4iPV+X/B6lyPh96+bQON0KPOyE9GWaNvXx+l5ug682QdAq+nfDXebf9UOjujh49SkBAwMU7agMC8u9B5uXlxauvvkrLli3ZvHnzX+rxAPTt25ctW7Zctl7uks/z5s1jyZIlzJgxg0mTJuHt7U1kZCTnzp3jzTffZNiwYcyYMYOFCxdy9uxZ4uPjueWWWy6Whv7qq6945513SE9Pp2fPnnzwwQe4u7vj5+fH3XffzfLly3n//ffp06fPX2L48ssvueOOO8jMzGT69On06NGD06dPM3nyZA4cOICvry+ffPIJnTp14plnnsHPz49HHnkEgNDQUJYssYYeHzJkCH369OG3334jMDCQ7777Dh8fH6Kioi6Wor7mmmsK/bztZZp6DMaFB3Mu1WqPL0pGVjYrdh7nvllRdH9xOU8u2saFtEyeGNKW9U8M4ovJPRjRJdAkfQOwElVsbCytW7fmvvvuu3g0nx93d3c6d+7Mrl27/jI9MzOTZcuWXSyiZq9Dhw6xYcMGfvjhB+655x5SU1MBa8Sv+fPns2XLFubOnUtkZCQ7d+5kzpw5rFu3jpiYGNzd3Zk1axYASUlJ9OzZk82bN1+W9MEq6hYTE8MHH3xwMUE//fTTdO3alS1btvDSSy8xYcKEy9bLa+/evdx///1s376dWrVqMX/+fABuu+023n33XTZv3lys918Yc8Rv0Kt5XYLr+PDNhlhGdAm8bL6qsv3IOVsXzHhOXkinTnUvburZhDHdgujQuKa5m7ayKOwI3cu38PnV6xZ5hJ+Xn58fUVFR/Prrr6xatYpx48bxyiuvMGnSpHyXz11CJqfGD1hH/Lfffnux9h0REYGbmxutWrWiefPmF79QBg8eTN26dQEYNWoUa9euxcPDg6ioKLp3735x3/Xr1wesL6TRo0cXuJ/x48cD0K9fP86dO0diYiJr1669mLgHDhzIqVOnOHeu8EFwmjVrdvH9duvWjUOHDpGYmEhiYiL9+vUDrMFoli1bVqzPIT+OHGx9OjAMSFDV0Dzz/g68DtRT1ZOOisGwj5ubENEtmDd+3sPhU0k0rVsdgIRzqSyKiWdBdDy7jp3Hy92NQe3qMyosiP5t6uFpumAadnB3d6d///7079+fjh07MnPmzHwTf1ZWFlu3br1YbyZ3G39Bch9w5BzR5zcv9+v8pqsqEydO5OWXX75sH97e3oW26xe0n/x4eHj8ZUzf3DHnLjDn7u5OSkpKgdspLUf+584ALiuxJyLBwDXAnw7ct1FMY8KDcBP46vfDfL/5CBOnb6DXyyt4aekuvD3def7GUDb8axAf3tKNwe0bmKRv2GX37t3s3bv34uuYmBiaNm162XIZGRk88cQTBAcH06lTJ7u336BBA3bu3El2djYLFy78y7y5c+eSnZ3N/v37OXDgAG3atAHg559/5vTp06SkpLBo0SKuvPJKBg0axLx580hISADg9OnTHD582K4Y5syZA8DatWvx9/fH39+fvn37XmwqWr16NQEBAdSsWZOQkBCio6MBiI6O5uDBg4Vuu1atWtSqVYu1a9cCXNxmaTnsiF9V14hISD6z3gIeA75z1L6N4mvk70O/1vX49FfrD7Gxvzf39m/BqLAgWtTzc3J0RmV14cIFpkyZQmJiIh4eHrRs2fIv49nefPPNVKtWjbS0NK6++mq++654aeGVV15h2LBh1KtXj/Dw8L+M7dukSRN69OjBuXPn+Oijj/D2tkqA9OjRg9GjRxMXF8ctt9xycSCXF154gWuuuYbs7Gw8PT15//338/2Sysvb25uuXbuSkZHB9OnTAXjmmWeYPHkynTp1wtfXl5kzZwIwevRovvjiCzp06EDPnj1p3bp1kdv//PPPmTx5MiJSZhd3HVqW2Zb4l+Q09YjICGCgqk4VkUNAeEFNPSJyF3AXQJMmTbrZ++1rlNy2+LPMjYzl2g4N6dW8rumCWQVU5bLMhZk0aRLDhg1jzJgxf5k+Y8YMIiMjee+995wUmeNUyLLMIuIL/BOrmadIqvoJ8AlY9fgdGJphExroX6q+/IZhVA7l2aunBdAM2Gy7+BEERItID1Utuh+hYRiGnWbMmJHv9EmTJhXYo8iVlFviV9WtQP2c10U19RiGUTZU1XS3reKK22TvsK4ZIjIbWA+0EZE4ESleJ1zDMErN29ubU6dOFTsxGJWHqnLq1KmLF6/t4chePeOLmB/iqH0bhmEJCgoiLi6OEydOODsUw4G8vb0JCgqye3lz565hVGGenp40a9bM2WEYFYy5C8cwDMPFmMRvGIbhYkziNwzDcDEOvXO3rIjICaCkt+4GABWxy6iJq3hMXMVj4iqeihoXlC62pqpaL+/ESpH4S0NEIvO7ZdnZTFzFY+IqHhNX8VTUuMAxsZmmHsMwDBdjEr9hGIaLcYXE/0nRiziFiat4TFzFY+IqnooaFzggtirfxm8YhmH8lSsc8RuGYRi5mMRvGIbhYqpM4heR60Rkt4jsE5HH85lfTUTm2Ob/UcCwkM6Ia5KInBCRGNvjjnKIabqIJIjItgLmi4i8Y4t5i4iEOTomO+PqLyJnc31WT5VTXMEiskpEdojIdhGZms8y5f6Z2RlXuX9mIuItIhtEZLMtrmfzWabc/x/tjKvc/x9z7dtdRDaJyJJ85pXt56Wqlf4BuAP7geaAF7AZaJ9nmfuAj2zP/w+YU0HimgS8V86fVz8gDNhWwPzrgWWAAL2APypIXP2xhvIs77+vRkCY7XkNYE8+v8dy/8zsjKvcPzPbZ+Bne+4J/AH0yrOMM/4f7Ymr3P8fc+37YeDr/H5fZf15VZUj/h7APlU9oKrpwDfAiDzLjABm2p7PAwaJ40ensCeucqeqa4DThSwyAvhCLb8DtUSkUQWIyylU9aiqRtuenwd2AoF5Fiv3z8zOuMqd7TPIGfXc0/bI24uk3P8f7YzLKUQkCBgKfFbAImX6eVWVxB8IxOZ6Hcfl/wAXl1HVTOAsULcCxAUw2tY8ME9Egh0ckz3sjdsZettO1ZeJSIfy3rntFLsr1tFibk79zAqJC5zwmdmaLWKABOBnVS3w8yrH/0d74gLn/D++DTwGZBcwv0w/r6qS+CuzxUCIqnYCfubSt7pxuWis2iOdgXeBReW5cxHxA+YDD6nqufLcd2GKiMspn5mqZqlqF6yxtXuISGh57LcodsRV7v+PIjIMSFDVKEfvK0dVSfzxQO5v5iDbtHyXEREPwB845ey4VPWUqqbZXn4GdHNwTPaw5/Msd6p6LudUXVWXAp4iElAe+xYRT6zkOktVF+SziFM+s6LicuZnZttnIrAKuC7PLGf8PxYZl5P+H68EbhBrHPJvgIEi8lWeZcr086oqiX8j0EpEmomIF9bFj+/zLPM9MNH2fAywUm1XSpwZV5524Buw2mmd7Xtggq2nSi/grKoedXZQItIwp11TRHpg/f06PFnY9jkN2KmqbxawWLl/ZvbE5YzPTETqiUgt23MfYDCwK89i5f7/aE9czvh/VNUnVDVIreFo/w/rs7glz2Jl+nlViaEXVTVTRB4AfsLqSTNdVbeLyHNApKp+j/UP8qWI7MO6gPh/FSSuB0XkBiDTFtckR8clIrOxensEiEgc8DTWhS5U9SNgKVYvlX1AMnCbo2OyM64xwL0ikgmkAP9XDl/eYB2R3QpstbUPA/wTaJIrNmd8ZvbE5YzPrBEwU0Tcsb5ovlXVJc7+f7QzrnL/fyyIIz8vU7LBMAzDxVSVph7DMAzDTibxG4ZhuBiT+A3DMFyMSfyGYRguxiR+wzAMF2MSv1GpiEiWrWriNhFZnNMvuwy2O0lE3iuLbeXZroeIvCQie3NVfPxXGW5/hoiMKavtGa7BJH6jsklR1S6qGorVn/l+ZwdUhBeAxkBHW6mAvtjuTcjNduOX+X80yoX5QzMqs/XYCqGJSA8RWS9WPfPfRKSNbfokEVkgIj/ajrpfzVlZRG4TkT0isgHrZqic6SEistJWqGuFiDSxTZ8hIh+KyO8ickCsWvfTRWSniMzIG5yI+AJ3AlNUNRWsKpqq+kyu/ewWkS+AbUCwbfuRkqdevIgcEpFXRWSrWDXlW+baVT/bez5gjv4Ne5jEb1RKtrsvB3GpBMYuoK+qdgWeAl7KtXgXYBzQERgn1gAmjYBnsRJ+H6B9ruXfBWbaCnXNAt7JNa820Bv4m23fbwEdgI4i0iVPmC2BP20lkwvSCvhAVTuo6mHgX6oaDnQCrhKRTrmWPauqHYH3sKo55mhkew/DgFcK2ZdhACbxG5WPj608wTGgAVYFRbCKVs0Va/SunGScY4WqnrUdde8AmgI9gdWqesI2VsKcXMv3xhoQA+BLrKSaY7Gt5MFW4LiqblXVbGA7EFJY4LYzjBgRiZVL5X4P2+r354gQkWhgk+095P5Cmp3rZ+9c0xeparaq7rB9JoZRKJP4jcomxdZW3hRrRKWcNv7ngVW2tv/hgHeuddJyPc+idDWqcraVnWe72flsdx/QRERqAKjq57bYz2LVbgJIyllYRJoBjwCDbGcbP+R5H1rA89xxOHpwIaMKMInfqJRUNRl4EPi7XCpTm1MGeZIdm/gDqymlrliljcfmmvcbl4pg3Qz8WooYpwHviYg3XGyi8ipglZpYXwRnRaQBMCTP/HG5fq4vSUyGAVWkOqfhmlR1k4hsAcYDr2JVXnwS60i5qHWPisgzWAk0EYjJNXsK8LmIPAqcoHSVNv+FdTayTUTOY1XInAkcwertkzumzSKyCet6RSywLs+2atvebxrWezaMEjHVOQ2jEhBrkI5wVT3p7FiMys809RiGYbgYc8RvGIbhYswRv2EYhosxid8wDMPFmMRvGIbhYkziNwzDcDEm8RuGYbiY/wcxKxz5n6RYrgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "plt.plot(value_dc_qaoa_ls, label=\"DC-QAOA\")\n",
    "plt.plot(ubound_sdp_ls, label=\"SDP upper bound\", linestyle=\"--\")\n",
    "plt.title('Max-Cut Performancce')\n",
    "plt.xlabel('Random Graph')\n",
    "plt.ylabel('Calculated Optimal Max Cut')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the line graph above, we have verified that max-cut generated by DC-QAOA is under and close to this bound.\n",
    "\n",
    "## Applications\n",
    "\n",
    "The Max-Cut problem belongs to quadratic unconstrained binary optimization (QUBO), which has a wide range of applications, from finding ground states for the spin glass problem to modeling NP-hard problems [5]. The Max-Cut problem itself inherits widespread usefulness.\n",
    "\n",
    "It is extensively related to various fields, including VLSI circuit design, statistical physics. Both the problem of minimizing the number of vias subject to pin preassignments and layer preferences and the problem of finding ground states of spin glasses with exterior magnetic field in the Ising model can be reduced to the Max-Cut problem [6]. \n",
    "\n",
    "Also importantly, the Max-Cut problem provides a prototypical testbed for algorithmic techniques that can be applied to many interesting problems. For example, SDP relaxation of the Max-Cut problem is adopted in designing data clustering algorithms [7] and addressing the phase retrieval problem [8, 9].\n",
    "\n",
    "More detailed investigations on the Max-Cut problem can be found in [10-12]."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "## References\n",
    "\n",
    "[1] Akshay, V., et al. \"Reachability deficits in quantum approximate optimization.\" [Physical Review Letters 124.9 (2020): 090504.](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.124.090504)\n",
    "\n",
    "[2] Li, Junde, Mahabubul Alam, and Swaroop Ghosh. \"Large-scale Quantum Approximate Optimization via Divide-and-Conquer.\" [arXiv preprint arXiv:2102.13288 (2021).](https://arxiv.org/abs/2101.03717)\n",
    "\n",
    "[3] Goemans, Michel X., and David P. Williamson. \"Improved approximation algorithms for maximum cut and satisfiability problems using semidefinite programming.\" [Journal of the ACM (JACM) 42.6 (1995): 1115-1145.](http://www-math.mit.edu/~goemans/PAPERS/maxcut-jacm.pdf)\n",
    "\n",
    "[4] Burer, Samuel, and Renato DC Monteiro. \"Local minima and convergence in low-rank semidefinite programming.\" [Mathematical Programming 103.3 (2005): 427-444.](https://link.springer.com/article/10.1007/s10107-004-0564-1)\n",
    "\n",
    "[5] Kochenberger, Gary, et al. \"The unconstrained binary quadratic programming problem: a survey.\" [Journal of Combinatorial Optimization 28.1 (2014): 58-81.](https://link.springer.com/article/10.1007/s10878-014-9734-0)\n",
    "\n",
    "[6] Barahona, Francisco, et al. \"An application of combinatorial optimization to statistical physics and circuit layout design.\" [Operations Research 36.3 (1988): 493-513.](https://www.jstor.org/stable/170992?seq=1)\n",
    "\n",
    "[7] Poland, Jan, and Thomas Zeugmann. \"Clustering pairwise distances with missing data: Maximum cuts versus normalized cuts.\" [International Conference on Discovery Science. Springer, Berlin, Heidelberg, 2006.](https://link.springer.com/chapter/10.1007/11893318_21)\n",
    "\n",
    "[8] Candes, Emmanuel J., et al. \"Phase retrieval via matrix completion.\" [SIAM review 57.2 (2015): 225-251.](https://epubs.siam.org/doi/10.1137/110848074)\n",
    "\n",
    "[9] Waldspurger, Irene, Alexandre d’Aspremont, and Stéphane Mallat. \"Phase recovery, maxcut and complex semidefinite programming.\" [Mathematical Programming 149.1 (2015): 47-81.](https://link.springer.com/article/10.1007/s10107-013-0738-9)\n",
    "\n",
    "[10] Deza, Michel, and Monique Laurent. \"Applications of cut polyhedra—I.\" [Journal of Computational and Applied Mathematics 55.2 (1994): 191-216.](https://www.sciencedirect.com/science/article/pii/0377042794900205)\n",
    "\n",
    "[11] Deza, Michel, and Monique Laurent. \"Applications of cut polyhedra—II.\" [Journal of Computational and Applied Mathematics 55.2 (1994): 217-247.](https://www.sciencedirect.com/science/article/pii/0377042794900213)\n",
    "\n",
    "[12] Poljak, Svatopluk, and Zsolt Tuza. \"Maximum cuts and largest bipartite subgraphs.\" [DIMACS Series 20 (1995): 181-244.](https://arxiv.org/pdf/1810.12144.pdf)"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Raw Cell Format",
  "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.7.11"
  },
  "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
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}