QAOA.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Quantum Approximate Optimization Algorithm (QAOA)\n",
    "________________________\n",
    "\n",
    "This demo code reproduces the QAOA for MaxCut problem example from [Qiskit](https://qiskit.org/textbook/ch-applications/qaoa.html) and [Pennylane](https://pennylane.ai/qml/demos/tutorial_qaoa_maxcut.html) tutorials. \n",
    "\n",
    "Here I only present the explicit Tequila code. Check the above tutorials for details about the problem and theory."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import everything we need\n",
    "import tequila as tq\n",
    "import random\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Some useful functions\n",
    "\n",
    "# Exact ground state from a general Hamiltonian\n",
    "def exact(ham):\n",
    "    ham_matrix = ham.to_matrix()\n",
    "    energ = np.linalg.eigvals(ham_matrix)\n",
    "    return(min(energ))\n",
    "\n",
    "# Exact eigenstates and eigenvalues of a general Hamiltonian\n",
    "def exact_state(ham):\n",
    "    ham_matrix = ham.to_matrix()\n",
    "    eigenValues, eigenVectors = np.linalg.eig(ham_matrix)\n",
    "    idx = eigenValues.argsort()[::-1]   \n",
    "    eigenValues = eigenValues[idx]\n",
    "    eigenVectors = eigenVectors[:,idx]\n",
    "    return(eigenValues,eigenVectors)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Problem Hamiltonian\n",
    "def HP(param,graph):\n",
    "    \"\"\"\n",
    "    param (list): Problem Hamiltonian parameters\n",
    "    graph (2D list): Hamiltonian connectivity\n",
    "    \"\"\"\n",
    "    ham = tq.QubitHamiltonian() \n",
    "    for k in range(len(graph)):\n",
    "        ham -= param[k]*0.5*(1-tq.paulis.Z(graph[k][0])*tq.paulis.Z(graph[k][1]))\n",
    "    return(ham)\n",
    "\n",
    "# commuting part (Pauli string) from the problem Hamiltonian\n",
    "def HP_str(q1,q2):\n",
    "    ham_prob = \"Z({})Z({})\".format(q1,q2)\n",
    "    return ham_prob\n",
    "\n",
    "# commuting part (Pauli string) from the mixing Hamiltonian\n",
    "def HM_str(q1):\n",
    "    ham_mix = \"X({})\".format(q1)\n",
    "    return ham_mix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def QAOA_circuit(qubits,layers,graph):\n",
    "    \n",
    "    circuit = tq.QCircuit()\n",
    "    \n",
    "    #init state (eigenstate from the mixing Hamiltonian, |+>^n state for this example)\n",
    "    circuit += tq.gates.H(range(qubits))\n",
    "    \n",
    "    for l in range(layers):\n",
    "        beta = tq.Variable(name=\"beta{}\".format(l))\n",
    "        gamma = tq.Variable(name=\"gamma{}\".format(l)) \n",
    "        for k in range(qubits):\n",
    "            circuit += tq.gates.ExpPauli(angle=beta, paulistring=HM_str(k))\n",
    "        for k in range(len(graph)):\n",
    "            circuit += tq.gates.ExpPauli(angle=gamma, paulistring=HP_str(graph[k][0],graph[k][1]))\n",
    "    \n",
    "    return circuit\n",
    "\n",
    "# VQE hardware-efficient: layer= single-qubit Ry gates + alternating CNOTs\n",
    "# for comparison with QAOA\n",
    "def VQE_HE(qubits,layers):\n",
    "    \n",
    "    circuit = tq.QCircuit()\n",
    "    \n",
    "    for l in range(layers):\n",
    "        for i in range(qubits-1):\n",
    "            th = tq.Variable(name=\"th{}{}\".format(i,l))\n",
    "            circuit += tq.gates.Ry(target=i, angle=th)\n",
    "        # Layer alternating CNOTs\n",
    "        ent0 = np.mod(l, 2)\n",
    "        for ent in range(ent0, qubits - 1, 2):\n",
    "            circuit += tq.gates.CNOT(control=ent, target=ent + 1)\n",
    "    \n",
    "    return circuit\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "parameters =  [1, 1, 1, 1]\n",
      "graph:  [[0, 1], [0, 3], [1, 2], [2, 3]]\n",
      "exact energy =  (-4+0j)\n",
      "QAOA) 1 layers (2 parameters): -2.0\n",
      "QAOA) 2 layers (4 parameters): -2.999999635412692\n",
      "QAOA) 3 layers (6 parameters): -3.9999990931931464\n",
      "QAOA) 4 layers (8 parameters): -3.9999994309843405\n",
      "QAOA) 5 layers (10 parameters): -3.9999378551577838\n",
      "QAOA) 6 layers (12 parameters): -3.9999883750793974\n",
      "+51.0000|1010> +49.0000|0101> \n",
      "VQE) 1 layers (3 parameters): -1.999999992820456\n",
      "VQE) 2 layers (6 parameters): -3.999999646750707\n",
      "VQE) 3 layers (9 parameters): -3.9990771404706233\n",
      "+100.0000|0101> \n"
     ]
    }
   ],
   "source": [
    "qubits = 4 # number of qubits\n",
    "\n",
    "# To check the results as we consider more layers (variational parameters)\n",
    "lay_max_VQE = 3 # number of maximal layers for VQE hardware efficient ansatz\n",
    "lay_max_QAOA = 2*lay_max_VQE # number of maximal layers for QAOA ansatz\n",
    "\n",
    "# Optimizer options\n",
    "methods = 'COBYLA'\n",
    "grad_methods = '2-point'\n",
    "backend = 'qulacs'\n",
    "lr = 0.01\n",
    "# make sure scipy version is >= 1.5\n",
    "mthd_opt = {'finite_diff_rel_step': 0.0001}\n",
    "\n",
    "# MAX-CUT problem example\n",
    "graph = [[0, 1], [0, 3], [1, 2], [2, 3]]\n",
    "# assuming random parameters\n",
    "#param = [random.uniform(0,1) for i in range(len(graph))]\n",
    "# assuming param=1 (original Max-cut example)\n",
    "param = [1 for i in range(len(graph))]\n",
    "\n",
    "print(\"parameters = \", param)\n",
    "print(\"graph: \", graph)\n",
    "print(\"exact energy = \", exact(HP(param,graph)))\n",
    "#print(\"exact states = \", exact_state(HP(qubits,param,graph)))\n",
    "\n",
    "energ_QAOA = []\n",
    "energ_VQE = []\n",
    "param_QAOA=[]\n",
    "param_VQE=[]\n",
    "\n",
    "for layers in range(1,lay_max_QAOA+1):\n",
    "    \n",
    "    # Objective function\n",
    "    Obj_QAOA = tq.ExpectationValue(H=HP(param,graph), U=QAOA_circuit(qubits,layers,graph))\n",
    "    var_QAOA = Obj_QAOA.extract_variables() # how many variational parameters contain\n",
    "    \n",
    "    # Random initialization of variables\n",
    "    var0 = {key: random.uniform(0, np.pi) for key in var_QAOA}\n",
    "    initial_values = var0\n",
    "    \n",
    "    # classical optimization\n",
    "    QAOA = tq.minimize(objective=Obj_QAOA, adaptive = True, lr=lr, method_options=mthd_opt, method=methods, gradient=grad_methods, samples=None,\n",
    "                          initial_values=initial_values, backend=backend, noise=None, device=None, silent=True)      \n",
    "    \n",
    "    energ_QAOA.append(QAOA.energy) # energy solution\n",
    "    # QA0A.angles # solution (variational parameters)\n",
    "    param_QAOA.append(len(var_QAOA))\n",
    "    \n",
    "    print(f\"QAOA) {layers} layers ({len(var_QAOA)} parameters): {QAOA.energy}\")\n",
    "\n",
    "# Simulate wavefunction final number of layers \n",
    "wfn = tq.simulate(QAOA_circuit(qubits,layers,graph), samples=100, variables=QAOA.angles, backend=backend)\n",
    "print(wfn)\n",
    "\n",
    "\n",
    "# Hardware-efficient VQE\n",
    "\n",
    "for layers in range(1,lay_max_VQE+1):\n",
    "    \n",
    "    Obj_VQE = tq.ExpectationValue(H=HP(param,graph), U=VQE_HE(qubits,layers))\n",
    "    var_VQE = Obj_VQE.extract_variables()\n",
    "    \n",
    "    # Random initialization of variables\n",
    "    var0 = {key: random.uniform(0, np.pi) for key in var_VQE}\n",
    "    initial_values = var0\n",
    "    \n",
    "    VQE = tq.minimize(objective=Obj_VQE, adaptive = True, lr=lr, method_options=mthd_opt, method=methods, gradient=grad_methods, samples=None,\n",
    "                          initial_values=initial_values, backend=backend, noise=None, device=None, silent=True)      \n",
    "    \n",
    "    energ_VQE.append(VQE.energy)\n",
    "    param_VQE.append(len(var_VQE))\n",
    "    \n",
    "    print(f\"VQE) {layers} layers ({len(var_VQE)} parameters): {VQE.energy}\")\n",
    "    \n",
    "wfn = tq.simulate(VQE_HE(qubits,layers), samples=100, variables=VQE.angles, backend=backend)\n",
    "print(wfn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEGCAYAAABLgMOSAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAABBwUlEQVR4nO3dd3gVZfbA8e9JL4Sa0EsglFClBBRiuYqwiCxFdAVRsbL+FIV1XRsuisquylrXdV3boisiroqAIlJsCwgSEOlNioQSeicQkvP7YyYxxDSSW1LO53nuc+87887MmSvek2nnFVXFGGOMOVdBgQ7AGGNM+WQJxBhjTIlYAjHGGFMilkCMMcaUiCUQY4wxJRIS6AD8KTY2VuPj4wMdhjHGlCtLly7dp6pxeadXqgQSHx9PSkpKoMMwxphyRUS25TfdTmEZY4wpEUsgxhhjSsQSiDHGmBKpVNdASupIegb3ffAjf/vdeVSNCA10OMYYH8nIyCA1NZX09PRAhxIQERERNGzYkNDQ4v3OWQIphrlr0pi9Jo15a9MY1KlhoMMxxvhIamoqMTExxMfHIyKBDsevVJX9+/eTmppK06ZNi7VMQE5hicgEEVknIitEZKqIVC+gXx8RWS8im0TkwVzTa4rIHBHZ6L7X8GW8H6Rsd96XpPpyM7+WlQnrZ8E3zzjvWZn+3b4xlUx6ejq1atWqdMkDQESoVavWOR19BeoIZA7wkKqeEZGngYeAB3J3EJFg4B9ALyAVWCIi01V1DfAgME9Vn3ITy4N5ly+NYa8vYsFP+3PaocHOP6aUbQeIf/CznOnJCbWYdPsF3trs2bIy4T+DYEcKnD4BYVHQIAlumApBwb7ZpjGmUiaPbOe67wE5AlHV2ap6xm0uAvI7L9QN2KSqm1X1NPA+MMCdNwB42/38NjDQm/HddVlzIkN/+ZHOyNSz3gEiQ4MZeVkLb272bBvnQOoSOH0cUOd9R4oz3RhTZhxJz2DEOykcSc8IdCh+VxbuwroF+Dyf6Q2A7bnaqe40gDqqugvAfa9d0MpFZISIpIhIyt69e4sVUI+EWN68KemsJJJbZGgwb93Ule4JtYq1vhLZvQIyTp497fQJ2L3Sd9s0xpyz3NdIvSE1NZUBAwbQokULmjVrxsiRIzl16lTO/FGjRtGgQQOysrLOWu61114jMTGRxMREunXrxvz588+a/8MPPyAifPHFF16JE3yYQERkroisyuc1IFefMcAZYFJ+q8hn2jmPfqWqr6lqkqomxcX96kn8AvVIiOXl6zoRHnL2VxQeEsTL13XybfIAqNvBOW2VW1gU1G3v2+0aY86JN6+RqipXXXUVAwcOZOPGjWzcuJGTJ09y//33A5CVlcXUqVNp1KgR3377bc5yn376Kf/617+YP38+69at49VXX+W6665j9+7dOX0mT57MhRdeyOTJk0sdZzafXQNR1csLmy8iw4F+QE/Nf1jEVKBRrnZDYKf7OU1E6qnqLhGpB+zxRsx5HUnPICRIyBAIDhIyMpUsxT+Hqi16Odc8tn8PZ05CcKjTbtHL99s2xhTIl9dIv/zySyIiIrj55psBCA4O5vnnn6dJkyaMHz+exYsX065dO6699lomT56Mx+MB4Omnn2bChAnExsYC0LlzZ4YPH84//vEPnnjiCVSVDz/8kDlz5nDRRReRnp5OREREab4GIEAX0UWkD85F70tU9UQB3ZYALUSkKbADGAJc586bDgwHnnLfp/kizilLtnMiI5M29aryQJ9W3DnpB46dOsO7i7b5/nbeoGDngvnGOfDFQ3B8PwyZbBfQjfGTcTNWs2bnkV9NP3wygyCBLPfP3vyukQYJHDh+mmv/9d1Zy7apX5VHf9u2wG2uXr2aLl26nDWtatWqxMfHs2nTJiZPnszQoUMZMGAADz/8MBkZGYSGhua7XFJSEm+/7VwqXrBgAU2bNiUhIQGPx8PMmTO56qqriv9lFCBQ10BeBmKAOSKyXEReBRCR+iIyE8C9yD4S+AJYC3ygqqvd5Z8CeonIRpy7tJ7yRZAx4aE83Lc1M0ZeyMUtazNj5IUEC2w/cLLohb0hKBha9YF+L8Cpw7Dqv/7ZrjGmQNUiQ2lZJ4agAm5YChJoVSeGqpHn/tCxquZ7J5Sqcvr0aWbOnMnAgQOpWrUq559/PrNnzy7WuiZPnsyQIUMAGDJkiNdOYwXkCERVmxcwfSfQN1d7JjAzn377gZ4+C9D1+vCks9pN46K5t3crJnyxnq/W7eHSxAKv3XtX04uhXkdY8BJ0usGOQozxg8KOFADmrU3jzknLOHXml4vZ4SFBvDKsMz1b1ynRNtu2bctHH3101rQjR46QlpbGrl27OHz4MO3bO9dBT5w4QVRUFFdeeSVt2rRh6dKlXHbZZTnLLVu2jDZt2pCZmclHH33E9OnTGT9+fM4Dg0ePHiUmJqZEcWYrC3dhlSu3X9SMhLhoxk5fxcnTfnqwTwSSR8GBn2DdZ0X3N8b4XPY10iCBiNAggtxrpaW5RtqzZ09OnDjBO++8A0BmZiZ//OMfGTlyJO+//z5vvPEGW7duZevWrWzZsoXZs2dz4sQJ7r//fh544AH273euzSxfvpyJEydy5513MnfuXM477zy2b9/O1q1b2bZtG4MHD+aTTz4p9XdgCeQchYUE8eTA9mw/cJKXv9rovw23GQA14mHBC5DvPQfGGH/Kvkbaul5VXr8xidb1qnIyI7NUd2OJCFOnTuXDDz+kRYsW1KpVi6CgIP7whz/wxRdfcOWVV+b0jY6O5sILL2TGjBn079+fW265hR49epCYmMjtt9/Ou+++S7169Zg8eTKDBg06azuDBw/mvffeK3GcOfHmfwNUxZSUlKTeGlDq3g+WM+PHnXw+6iKa1y7dYWCxLXkDPvsj3DQT4pP9s01jKpG1a9fSunXrYvW9/e0UujWrya3JTQkKEjKzlDfnb2bJ1oO8fmNS0SsohoULFzJ06FA+/vjjX10k95X8vgMRWaqqv9opSyAltO/YKXo++w2JdWN4f8QF/il/kHESnm8HDTrDMLugboy3nUsCqajOJYHYKawSiq0SzgN9Elm85QAfL9vhn42GRsL5v4eNsyFtddH9jTHGhyyBlMKQro3o1Lg6f5m5lkMnTvtno11vg9AoWPh3/2zPGGMKYAmkFIKChPED23PoZAZPz1rvn41G1YTOw2Hlf+Gwn8vLG2NMLpZASqlN/arc3COeyd//zNJtB/2z0e53OndiffeKf7ZnjDH5sATiBaN7taRetQjGTF3JmcysohcoreqNod1gWDoRTvopaRljTB6WQLygSngIj/62Det2H2Xiwq3+2WjyKMg4Dkve9M/2jDE+5/F4flVu/YUXXuDOO+9k9erVXHbZZbRs2ZKEhAQeffTRnJLuEydOJC4ujo4dO+a81qxZ4/N4LYF4yW/a1uWyxNo8N2cDOw/5oVZW3XbQ/HJY/CpkFH8ISmOMF3l52OmhQ4fy/vvvnzXt/fffZ+jQofTv358HH3yQDRs2sHLlSr7//ntefPHFnH7XXnsty5cvz3m1adOmVLEUhyUQLxERxvVvS5Yqj8/wfeYHnKOQ43vhx9I/UWqMOUfZw05/dAt89Rfn/T+DSpVErr76aj799NOcAaS2bt3Kzp072bBhA8nJyfTu3RuAqKgoXn75ZSZMmOCVXSmpQI2JXiE1qhnF3Ze1YMIX65m3Nq3EBdWKLf4iqN/JuaW383ArsmiMN33+YOEjgJ44APvWgbrXPU8fh63/g39e6NwtmZ+67eGKgouH16pVi27dujFr1iwGDBjA+++/z7XXXptvufaEhAROnjzJoUOHAJgyZcpZoxB+9913REZGFmtXS8qOQLzs9oua0bx2FR6dvtr3xRZFIHk0HNgM6z717baMMWc7feyX5JFNs5zppZD7NFb26avCyrxny3sKy9fJA+wIxOvCQoIYP7Ad1762iL9/uZH7+yT6doOtfws1msL8F6B1fyepGGNKr5AjBcC55vHRLc6RR7awaOg7wRnHp4QGDhzIvffey7Jlyzh58iSdO3dm2bJlZw1hC7B582ZiY2OpXr16ibdVWnYE4gPnN6vF4M4Nee3bzWxIO+rbjQUFQ4+7Yecy2Dq/6P7GGO/IHnY6LBoQ590Lw05XqVIFj8fDLbfcwtChQwEYNmwY8+fPZ+7cuQCcPHmSe+65h3HjxpV2L0olIAlERCaIyDoRWSEiU0Wkej59GonIVyKyVkRWi8ioXPMeE5Ed7miGy0Wkb97lA+3hvolEh4fwyCer8HnByo7XQXQcLHix6L7GGO/IHnZ68Ftw6Rjn/YapXrkWOXToUH788cecUQQjIyNzBoRq2bIlsbGxJCcnM2zYsJxlpkyZctZtvAsXLix1HEVSVb+/gN5AiPv5aeDpfPrUAzq7n2OADUAbt/0YcN+5brdLly7qT5MXb9MmD3yq/03Z7vuNffOM6qNVVXet9P22jKmg1qxZE+gQimXq1KnatGlT3bp1q9fXnd93AKRoPr+pATkCUdXZ6ox5DrAIaJhPn12qusz9fBRnXPQG/ouy9H6X1IjObrHFg8d9XGwx6VYIjYaFL/l2O8aYgBs4cCCbN2+mSZMmAY2jLFwDuQX4vLAOIhIPdAIW55o80j0F9paI1Chk2REikiIiKXv37vVKwMUVFCSMH9SewyczeHrWOt9uLKomdBkOKz+EQz/7dlvGGIMPE4iIzBWRVfm8BuTqMwY4A0wqZD1VgI+A0ap6xJ38TyAB6AjsAp4taHlVfU1Vk1Q1KS4urvQ7do5a16vKrRc25f0l21m67YBvN3bBnc5dWIv+6dvtGGMMPkwgqnq5qrbL5zUNQESGA/2AYe45tl8RkVCc5DFJVT/Ote40Vc1U1SzgdaCbr/bDG0b1bEH9ahGMmbqKDF8WW6zeCNpdDUvfdh5yMsYYHwrUXVh9gAeA/qp6ooA+ArwJrFXV5/LMq5erOQhY5atYvSE6PIRH+7dl3e6j/HvBFt9uLPkeK7JojPGLQF0DeRnnzqo57m24rwKISH0Rmen2SQZuAC7L53bdZ0RkpYisAC4F/uDvHThXvdvU4fLWtXlh7kZ2+LLYYp220LyXW2TRD0UdjTGVVqDuwmquqo1UtaP7usOdvlNV+7qf56uqqGqHXP1muvNuUNX27rz+qrorEPtxLkSER3/rFFscN93H45knj4IT+2C5FVk0xjg++eQTr5d4Lwt3YVUajWpGMapnS2avSWPumjTfbSj+QmjQxSmyWMry0saYisESSAVw20VNaVnHKbZ44vSZohcoCRHnKOTgFlg73TfbMMb4xLvvvku3bt3o2LEjv//971m8eDEdOnQgPT2d48eP07ZtW1atWsWxY8fo2bMnnTt3pn379kybNi1nHe+88w4dOnTgvPPO44YbbmDhwoVMnz6dP/3pT3Ts2JGffvrJK7FaAvGz0OAgnhzYnh2HTvLSvE2+21BiP6jZzCmy6OtSKsZUUB6Ph4kTJwKQkZGBx+Ph3XffBeDEiRN4PB6mTJkCwOHDh/F4PHz8sXPD6L59+/B4PMyYMQOA3bt3F7m9tWvXMmXKFBYsWMDy5csJDg5m/fr19O/fn0ceeYT777+f66+/nnbt2hEREcHUqVNZtmwZX331FX/84x9RVVavXs348eP58ssv+fHHH3nxxRfp0aMH/fv3Z8KECSxfvpyEhASvfD9WjTcAujWtyTVdGvLG/zYzqFMDWtWN8f5GsossfvoHZ4yCphd7fxvGGK+aN28eS5cupWvXroBTNLF27dqMHTuWrl27EhERwUsvOdUmVJWHH36Yb7/9lqCgIHbs2EFaWhpffvklV199NbGxsQDUrFnA2CReYAkkQB7q25o5a9N45JOVTBnRnaAgH5RhP+86Z6S0BS9aAjGmBL7++uucz6GhoWe1o6KizmpXq1btrHZsbOxZ7bp16xa5PVVl+PDh/PWvfz1r+u7duzl27BgZGRmkp6cTHR3NpEmT2Lt3L0uXLiU0NJT4+HjS09MLHDvEF+wUVoDUjA7joSsSWbL1IB8uS/XNRkIj4Pw7YNPcwkdWM8aUCT179uTDDz9kz549ABw4cIBt27YxYsQInnjiCYYNG8YDDzwAOKfMateuTWhoKF999RXbtm3LWccHH3zA/v37c9YBEBMTw9Gj3h1ewhJIAF3TpRFJTWrwV18WW+zqFllcYEUWjSnr2rRpw5NPPknv3r3p0KEDvXr14u233yYkJITrrruOBx98kCVLlvDll18ybNgwUlJSSEpKYtKkSSQmOoPXtW3bljFjxnDJJZdw3nnnce+99wIwZMgQJkyYQKdOnbx2EV0KqCJSISUlJWlKSkqgwzjLut1H6PfSfAZ3bsjTV3fwzUZmPew8WDhqOVRv7JttGFMBrF27ltatWwc6jIDK7zsQkaWqmpS3rx2BBFhiXafY4pSU7SzZ6qP6Vd3dIovfveKb9RtjKiVLIGXAqMtb0KB6JI/4qthitYbQ/hpYZkUWjTHeYwmkDIgKC+Gx/m1Zn3aUt+b7qNhij3sg4wQsecM36zemgqhMp/XzOtd9twRSRvRqU4fLW9fhhbkbST2Yb4Hi0qnTBlr8xoosGlOIiIgI9u/fXymTiKqyf/9+IiIiir2MXUQvQ1IPnqDXc9+S3DyWN4b/6npV6W1dABP7wpXPQtfbvL9+Y8q5jIwMUlNTSU9PD3QoAREREUHDhg0JDQ09a3pBF9HtQcIypGGNKEZf3oK/fr6O2at307tt0Q8enZMmPaBBklNksfNNEGz/+Y3JLTQ0lKZNmwY6jHLDTmGVMbdc2JRWdWJ4bPpqjp/ycrHFnCKLW63IojGm1CyBlDGhwUE8OagdOw+n89K8jd7fQOKVUKu5U96kEp2+NMZ4X6CGtJ0gIutEZIWITBWR6gX02+qOPLhcRFJyTa8pInNEZKP7XsNvwftB1/iaXJvUiDfnb2Hd7iPeXXl2kcVdy2HLt95dtzGmUgnUEcgcoJ2qdgA2AA8V0vdSdzTC3BdwHgTmqWoLYJ7brlAevCKRmIgQHpm6iqwsLx8pdBgC0bVhwQveXa8xplIJ1JC2s1U1+wT/IqDhOa5iAPC2+/ltYKCXQiszakSH8VDf1qRsO8h/l2737spDI+CCO+CnL2HXCu+u2xhTaZSFayC3AJ8XME+B2SKyVERG5JpeJ3scdPe9dkErF5ERIpIiIil79+71WtD+cHXnhnSLr8lfP1/HAW8XW0y6FcKqwEIrsmiMKRmfJRARmSsiq/J5DcjVZwxwBphUwGqSVbUzcAVwl4ic86AWqvqaqiapalJcXFyJ9iVQgoKEJwe141j6Gf46c613Vx5ZHbrcBKs+hoPbvLtuY0yl4LMEoqqXq2q7fF7TAERkONAPGKYFPM2oqjvd9z3AVKCbOytNROq566kH7PHVfgRayzox3HZRM/67NJXvt3i5jtUF2UUW/+Hd9RpjKoVA3YXVB3gA6K+q+dbtEJFoEYnJ/gz0Bla5s6cDw93Pw4Fpv15DxXFPz+ZOscVPVnL6jBeLLVZrAO1/B8vegeP7vbdeY0ylEKhrIC8DMcAc9xbdVwFEpL6IzHT71AHmi8iPwPfAZ6o6y533FNBLRDYCvdx2hRUVFsK4/m3ZkHaMN71dbDH5Hjhz0oosGmPOWUBqWahq8wKm7wT6up83A+cV0G8/0NNnAZZBl7epQ+82dXhx3gb6dahHo5pR3llx7dbQsg98/y/n+ZAwL63XGFPhlYW7sEwxPdq/LUEiPDZ9tXerhSaPghP7YXlB9zIYY8yvWQIpRxpUj2T05S2Yt24Ps9ekeW/FjbtDw65OkcVML9ffMsZUWJZAypmbk5uSWDeGcd4stigCyaPh0DZYW6HvRzDGeJElkHImNDiI8W6xxRe9WWyxVV8rsmiMOSeWQMqhLk1qMqSrU2xx7S4vFVsMCnKGvd31I2z+2jvrNMZUaJZAyqkH+iRSLTKURz7xYrHFDtdClTrOUYgxxhTBEkg5VSM6jIf7tmbptoN8kOKlYouhEXDB/8Hmr5wjEWOMKYQlkHJscOcGdGvqFFvcf+yUd1ba5WYIi4EFVmTRGFM4SyDlmIgwfmA7jp86w18/X+edlUZWh6SbYPXHztC3xhhTAEsg5VyLOjGMuLgZHy5NZdFmL9WzuuBOkGArsmiMKZQlkArg7sta0LBGJI98sso7xRar1ncuqC/7jxVZNMYUyBJIBRAZFszjA9qyac8x3pi/2Tsr7XG3W2Txde+szxhT4VgCqSAuS6xDn7Z1eWneRrYfyLdC/rmpnQgtr4DF/4LTx0u/PmNMhWMJpAJ5tH8bgkV41FvFFpNHwckD8IMVWTTG/JolkAqkXrVI/tCrJV+u28MXq71QbLFJd2h0PnxnRRaNMb8WqBEJJ4jIOhFZISJTRaR6Pn1auYNNZb+OiMhod95jIrIj17y+/t6HsuqmHvG0rleVcTNWc8wbxRaTR8Ghn2HNJ6VflzGmQgnUEcgcoJ2qdgA2AA/l7aCq61W1o6p2BLoAJ3DGRc/2fPZ8VZ2Zd/nKKsQttrj7SDovzNlQ+hW2vAJqtYAFL1iRRWPMWQKSQFR1tqpm/3m8CGhYxCI9gZ9UdZtvI6sYOjeuwZCujfn3wq2s2VnKYotBQc6wt7tXOiVOjDHGVRaugdwCfF5EnyHA5DzTRrqnwN4SkRoFLSgiI0QkRURS9u7dW9pYy40H+rSiemQoYz5ZWfpiix2uhSp1rciiMeYsPksgIjJXRFbl8xqQq88Y4AxQ4G0+IhIG9Af+m2vyP4EEoCOwC3i2oOVV9TVVTVLVpLi4uNLtVDlSPSqMMVe25oefD/H+klIWWwwJd4ssfg07l3sjPGNMBeCzBKKql6tqu3xe0wBEZDjQDximhd9zegWwTFVzbitS1TRVzVTVLOB1oJuv9qM8G9SpARc0q8nTs9axr7TFFpOyiyzaUYgxxhGou7D6AA8A/VW1qKfehpLn9JWI1MvVHASs8m6EFYOI8OTAdpw4fYa/zFxbupVFVHOSyJpP4MAWr8RnjCnfAnUN5GUgBpjj3ob7KoCI1BeRnDuqRCQK6AV8nGf5Z0RkpYisAC4F/uCnuMud5rVj+P3FCXy8bAff/VTKulZWZNEYk4sU54ll94f8j0BjVb1dRFoArVT1U18H6E1JSUmakpIS6DD8Lj0jk17Pf0NYcBCfj7qYsJBS/N0w7S5Y+RH8YRVEx3ovSGNMmSUiS1U1Ke/04v6S/Bs4BXR326nAk16KzfhYRGgwj/dvx097j/P6/0pZbLHHPU6Rxe9f805wxphyq7gJJEFVnwEyAFT1JCA+i8p43aWJtenb3im2+PP+UhRbjGsFrfo6CcSKLBpTqRU3gZwWkUhAAUQkAeeIxJQjY/u1JSRIGDt9VemKLSaPhpMH4Yd3vRabMab8KW4CeQyYBTQSkUnAPOB+XwVlfKNutQju7d2Kr9fvZdaq3SVfUePzodEFsPBlK7JoTCVWrASiqrOBq4CbcG6pTVLVr30XlvGV4d2b0KZeVcbNWFO6YovJo+Dwz7B6atF9jTEVUrESiIhMB3oDX6vqp6q6z7dhGV/JLraYdjSd50tTbLFlH4ht5TxYaEUWjamUinsK61ngImCNiPxXRK4WkQgfxmV8qFPjGlzXrTH/XrCFVTsOl2wl2UUW01bCT196N0BjTLlQ3FNY36jqnUAz4DXgd8AeXwZmfOv+3yRSMzqMRz5ZVfJii+2vgZh6Vt7EmEqq2E+UuXdhDQbuALoCb/sqKON71aJCGXNla5ZvP8TkJT+XbCXZRRa3fAM7f/BugMaYMq+410CmAGuBy4B/4DwXcrcvAzO+N7BjA7o3q8XTn69j79ES3pXd5SYIr2pHIcZUQufyJHqCqt6hql+6VXBNOSciPDmoHekZWSUvthhRDZJugTXT4EApn3I3xpQrhSYQEbnM/RgFDBCRq3K/fB+e8bWEuCr8/pJmTP1hBwt/KuHNdeffAUEhVmTRmEqmqCOQS9z33+bz6ufDuIwf3XVpcxrXjOKRT1Zx6kzmua+gaj1n1MIf3oVjlWfUR2Mqu0ITiKo+6n58XFVvzv0CnvB9eMYfIkKDeXxAWzbvPc5r35TwNFTyKDhzyoosGlOJFPcayEf5TPvQm4GYwPK0qs2V7evx8leb2Lb/OEfSMxjxTgpH0jOKt4LYFpB4pZNATh3zbbDGmDKhqGsgiSIyGKiW5/rHTYA9SFjB/LlfG0KDgxg7bTVzVqcxe00a89amFb1gtuRRkH4IfviPz2I0xpQdRR2BtMK51lGds69/dAZuL+lGReQJEVnhjkY4W0TqF9Cvj4isF5FNIvJgruk1RWSOiGx032uUNBbzi7rVIvhj75Z8s2Evr37zEwAfLEkt/goadYPG3Z2L6ZnFPHIxxpRbxR2RsLuqfue1jYpUVdUj7ud7gDaqekeePsHABpwhbVOBJcBQVV0jIs8AB1T1KTex1FDVB4rabmUdkbA4hr2+iAX5DHkbGixkZP7ybyQ5oRaTbr+g4BWtnwWTr4WrXocOv/NFqMYYPyvtiIR3iEj1XCurISJvlTSY7OThisYdZySPbsAmVd2sqqeB94EB7rwB/PIk/NvAwJLGYhx3XdacyNDgX03PnTwiQ4MZeVmLwlfUojfEJVqRRWMqgeImkA6qeii7oaoHgU6l2bCIjBeR7cAwYGw+XRoA23O1U91pAHVUdZcbyy6gdiHbGSEiKSKSsnev3WJakB4Jsbx5U1K+SQSc5PHWTV3pnlCr8BUFBTnD3qatgk3zfBCpMaasKG4CCcp9nUFEagIhhS0gInNFZFU+rwEAqjpGVRsBk4CR+a0in2nn/Cetqr6mqkmqmhQXF3eui1cqPRJiefm6ToSHnP3PIjwkiJev61R08sjW/hqIqQ8LXvB+kMaYMqPQJJDLs8BCEfkQ50f8d8D4whZQ1cuLue73gM+AR/NMTwUa5Wo3BHa6n9NEpJ6q7hKRelhlYK85kp5BSJCQIRAc5Fz/UKX4t/MChIRB9zth9iOwYyk06OK7gI0xAVPccu7v4FTiTQP2Alepaonv1RSR3CfS+wPr8um2BGghIk1FJAwYAkx3500HhrufhwPTShqLOduUJds5kZFJ63pVeXN4EtFhwZzOzOLdRdvObUWdh0N4NVjwkm8CNcYEXLHLuQM1geOq+ndgr4g0LcV2n3JPZ63AGelwFICI1BeRmQCqegbn1NYXOJWAP1DV1dnLA71EZCPOXVpPlSIWk0tMeCgP923NjJEXcnHL2nxyVzJBAjsOnjy3FUVUha63wNrpsP8n3wRrjAmo4t7G+yiQBLRS1Zbucxv/VdVkXwfoTXYbb8k8N2cDL83byKTbzie5eWzxFzy6G15oD52uh37P+y5AY4xPlfY23kE4p5qOA6jqTiDGe+GZsuxOTwJNakXx53MtthhTF84bCj9MgmN2mcqYiqa4CeS0OocqCiAi0b4LyZQ1EaHBPDGgHZv3Hedf51psscfdkHnaiiwaUwEVN4F8ICL/AqqLyO3AXOB134VlypqLW8bRr4NTbHHrvuPFXzCnyOLrVmTRmAqmuHdh/Q2n+u5HOPWxxroX000l8ud+bQgPDuLP01ZRnGtnOZJHO0UWl73jq9CMMQFQ7LuwVHWOqv5JVe9T1Tm+DMqUTXWqOsUW/7dxH5+t3FX8BRt1hSbJVmTRmAqmqHLu8933oyJyJJ/XFhG50z+hmrLghu7xtG9QjcdnrDm3hwuTR8GRVFj1se+CM8b4VVEjEl7ovseoatW8L5xbe0f5I1BTNgQHCeMHtWPvsVM8N3tD8Rds3gviWluRRWMqkGKfwhKRziJyj4jcLSKdAFR1P+DxVXCmbOrQsDo3XNCEd77bysrUw8VbKCgIku+BPath01zfBmiM8YtiJRARGYtTNr0WEAtMFJFHIKcarqlk7vtNK2pVCWfMJyvJzCrmEUW7q6FqA+coxBhT7hX3CGQo0FVVH1XVR4ELcMqwm0qqakQof+7XhhWph5m0uJh1skLC4II7Yev/IHWpbwM0xvhccRPIVs4eAz0csAJHldxvO9TjwuaxTJi1nj1H04u3UJfsIosv+DQ2Y4zvFXUX1t9F5CXgFLBaRCaKyERgFWBPhVVyIsITA9txKjOLJz9dW7yFwmOg662wdoYVWTSmnCvqCCQFWArMwql4uwRYBIwBPvdtaKY8aBobzZ2eBKb/uJP/bSzmiI/n3wHBYbDQnkU1pjwrKoG8B7QFngRuAm5zP7dz5xnDHZckEF8rirHTVpOeUYxiizF1oONQWP6eFVk0phwrKoE8A9QAmqpqF1XtBDQDqgETfB2cKR8iQoN5YmA7tuw7zqvfFPO0VHe3yOLiV30bnDHGZ4pKIP2AEap6NHuCqh4B/g+40peBmfLlohZx9D+vPq989RNbilNsMbY5tO4HS96AU0eL7m+MKXOKSiCq+VTNU9VM3NLuJSEiT4jIChFZLiKz3QGq8vZpJCJfichaEVktIqNyzXtMRHa4yy8Xkb4ljcV4zyP9WhMeEsTY4hZbTB4N6YetyKIx5VRRCWSNiNyYd6KIXE/+45gX1wRV7aCqHYFPgbH59DkD/FFVW+M8d3KXiLTJNf95Ve3ovmaWIhbjJbVjIvhTn1b8b+M+ZqwoxvOlDZOgyYVWZNGYcqqoBHIXzg/31yLyrIj8TUS+Ae7BOY1VIu5psGzR5HM0o6q7VHWZ+/kozrjoDUq6TeMfw85vQoeG1Xji02IWW0weBUd2wMoPfR+cMcariiqmuENVzwcex3mY8GfgcVXtpqo7SrNhERkvIttxnmjP7wgkd994oBOwONfkke5psLdEpEYhy44QkRQRSdm7t5i3mZoSCw4Sxg9sz/5jp3j2i/VFL9CiF9RuY0UWjSmHijug1Jeq+ndVfUlV5xVnGRGZKyKr8nkNcNc5RlUbAZOAkYWspwrOQFajcx25/BNIADoCu4BnC4n9NVVNUtWkuLi44oRuSql9w2rc2D2edxZtY0XqocI7izhHIXvXwkYbZsaY8qTY1XjPlaperqrt8nlNy9P1PWBwfusQkVCc5DFJVXMGklDVNFXNVNUsnKF1u/lqP0zJ3Nu7JXFVwhkzdVXRxRbbDYaqDa3IojHljM8SSGFEpEWuZn/yuSAvIgK8CaxV1efyzKuXqzkIp7SKKUOyiy2u3HGYdxcVUWwxOBS63wnb5kNqin8CNMaUWkASCPCUezprBdAbd1AqEakvItl3VCUDNwCX5XO77jMistJd/lLgD/7eAVO0fh3qcVGLWCZ8sZ60I0UUW+x8I0RYkUVjyhMp1v36FURSUpKmpNhfuP60dd9xer/wLb3b1OHl6zoX3nneE/C/Z2FkivOgoTGmTBCRpaqalHd6oI5ATCURHxvNXZ7mfLpiF99uKOIuuPN/7xRZ/M6KLBpTHlgCMT53h6cZzWKj+fO0VYUXW6xSGzpeB8snw9E0/wVojCkRSyDG58JDnGKL2/af4JWviyi22MOKLBpTXlgCMX6R3DyWAR3r8+rXP7F5byFjkdVKgDb9YcmbVmTRmDLOEojxmzFXtiY8NIg/F1VsMXkUnDoMS9/2X3DGmHNmCcT4Te2YCO7vk8iCTfuZ/uPOgjs26ALxFzlFFs+c9l+AxphzYgnE+NV13RpzXsNqPPHpWg6fLKTYYvIoOLoTVlmRRWPKKksgxq+Cg4Txg9pz4Pgp/lZYscXml0PttrDgJcjK8l+AxphiswRi/K5dg2oM7xHPu4u3sXz7ofw75S6yuMmKLBpTFlkCMQFxb6+W1I4JZ8zUlZzJLOAIo91VTpHF+S/4NTZjTPFYAjEBERMRyth+bVm98wj/KajYYnAodL8Lfl4I27/3b4DGmCJZAjEB07d9XS5pGcezszcUXGyx840QUd1KvRtTBlkCMQEjIjw+oC0ZmVk8/uma/DuFV4Fut8O6z2DfRv8GaIwplCUQE1BNakUz8tLmfLZiF1+v35N/p25ukcWFL/k3OGNMoSyBmIAbcUkzmsVFM3ba6vyLLVaJg07D4Mf34ehu/wdojMlXoEYkfEJEVriDRM0WkfoF9NvqDhy1XERSck2vKSJzRGSj+17Df9EbbwsPCebJAe34+cAJXvlqU/6duo+ErDNWZNGYMiRQRyATVLWDqnYEPgXGFtL3UlXtmGcwkweBearaApjntk051qN5LIM6NeCf3/zEpj35FFuslQCt+8OStyD9iP8DNMb8SkASiKrm/gWIBs51WMQBQHalvbeBgV4IywTYw31bExkazJ8/KaDYYk6RxYl+j80Y82sBuwYiIuNFZDswjIKPQBSYLSJLRWRErul1VHUXgPteu5DtjBCRFBFJ2bu3iBHxTEDFxYRzf59Evtu8n2nL8ym22KAzNL0YFr1iRRaNKQN8lkBEZK6IrMrnNQBAVceoaiNgEjCygNUkq2pn4ArgLhG5+FzjUNXXVDVJVZPi4uJKvD/GP67r1piOjarz5GdrOHwin2KLyaPg6C5Y+V//B2eMOYvPEoiqXq6q7fJ5TcvT9T1gcAHr2Om+7wGmAt3cWWkiUg/AfS/g/k9T3gQFCeMHtePA8dM888W6X3dI6Al12jm39FqRRWMCKlB3YbXI1ewP/OqXQkSiRSQm+zPQG1jlzp4ODHc/DwfyJiVTjrWtX42bejTlve9/5oefD549M6fI4jrY+EVgAjTGAIG7BvKUezprBU5iGAUgIvVFZKbbpw4wX0R+BL4HPlPVWdnLA71EZCPQy22bCuTe3i2pExPBmKmrfl1sse0gqNbYypsYE2AhgdioqhZ2yqqv+3kzcF4B/fYDPX0WoAm4KuEhPPrbNvzfpGW8/d02br2w6S8zs4ssznoAfl4Mjc8PXKDGVGL2JLops/q0q4unVRzPzV7P7sN5ii12vgEia1h5E2MCyBKIKbNEhMf7t+NMlvL4p6vPnhkWDV3dIot7NwQmQGMqOUsgpkxrXCuKe3q2YObK3XyVt9ji+b+HkHA7CjEmQCyBmDLv9ouakRAXzdhpq84uthgdC52uhxVT4MiuwAVoTCVlCcSUeWEhQTw5sD3bD5zk5S/zFFvsfpcVWTQmQCyBmHKhe0ItrurcgH99+xOb9hz9ZUbNZtBmAKS8BemHAxegMZWQJRBTbjzctzVRYSE8krfYYvIoOHXEiiwa42eWQEy5EVslnAf6JLJo8wGm/rDjlxn1O0HTS2DRP+HMqcAFaEwlYwnElCtDujaiU+PqjP9sLYdO5KrIa0UWjfE7SyCmXAkKEsYPbM+hkxk888X6X2YkXAZ12zvlTazIojF+YQnElDtt6lfl5h7xvLf4Z5Zuc4stikDyaNi3ATbMKnR5Y4x3WAIx5dLoXi2pVy2CRz7JVWyxzUArsmiMH1kCMeVSdrHFtbuOMHHhVmdicAj0GAnbF8HPiwIanzGVgSUQU279pm1dLkuszXNzNrDz0ElnYqfrIbKmHYUY4weWQEy5JSKM69+WLFUen7HGmRgWDd1GwPqZsHd94SswxpSKJRBTrjWqGcXdl7Vg1urdfLkuzZnY7XYIibQii8b4WKCGtH1CRFaIyHIRmS0i9fPp08qdn/06IiKj3XmPiciOXPP6+n0nTJlx+0XNaF67CmOnrebk6cxfiiz+OAWO7Ax0eMZUWIE6Apmgqh1UtSPwKTA2bwdVXa+qHd0+XYATwNRcXZ7Pnq+qM/MubyqPsJAgxg9sR+rBk/z9y43OxO53gWY6T6cbY3wiIAlEVY/kakYDWlBfV0/gJ1Xd5ruoTHl2frNaDO7ckNf/t5mNaUehZlNn7PSUf1uRRWN8JGDXQERkvIhsB4aRzxFIHkOAyXmmjXRPg70lIjUK2c4IEUkRkZS9e/eWMmpTlj3cN5GosBDGZBdb7HEPnD7qJBFjjNf5LIGIyFwRWZXPawCAqo5R1UbAJGBkIesJA/oDuYsc/RNIADoCu4BnC1peVV9T1SRVTYqLiyv9jpkyq1aVcB66IpHvtxzgo2U7oH5HaOaxIovG+IjPEoiqXq6q7fJ5TcvT9T1gcCGrugJYpqppudadpqqZqpoFvA508/4emPLod0mN6Ny4On+Z6RZbTB4Fx3Y7oxYaY7wqUHdhtcjV7A+sK6T7UPKcvhKRermag4BV3ovOlGdBQcL4Qe05fDKDp2etg2aXQt0OsOAlK7JojJcF6hrIU+7prBVAb2AUgIjUF5GcO6pEJAroBXycZ/lnRGSlu/ylwB/8FLcpB1rXq8qtFzZl8vfbWfrzQecoZP9G2PB5oEMzpkKRs0Z2q+CSkpI0JSUl0GEYPzh+6gy9nvuGqpGhzLjrAkL/0QWq1IVbZzuVe40xxSYiS1U1Ke90exLdVEjR4SE82r8t63YfZeJ3qdD9bkj93oosGuNFlkBMhdW7TR0ub12b5+duYFezwVZk0RgvswRiKiwR4dHfOsUWH5u1Bc7/vXMdZE9h92wYY4rLEoip0BrVjGJUz5Z8sTqNb6oNsCKLxniRJRBT4d12UVNa1qnCw1/sIuO8YbDiAzi8I9BhGVPuWQIxFV5ocBBPDmzPjkMneSPzStAsWGxFFo0pLUsgplLo1rQm13RpyLPfp3MkoR+kTISThwIdljHlmiUQU2k81Lc1VSJCePJgL7fI4luBDsmYcs0SiKk0akaH8dAViXywoya7Y7vD4lchIz3QYRlTblkCMZXKNV0akdSkBo/u7wXH0qzIojGlYAnEVCpBQcKTg9ox71QiqREtnVt6rciiMSViCcRUOol1q3Lrhc146shvYP8mWP9ZoEMyplyyBGIqpVGXt2BFzCXskjpkzX8BKlFRUWO8xRKIqZSiwkL484AO/ON0X4J2pMDP3wU6JGPKHUsgptLq1aYOB5pfzQGN4eTXz3EkPYMR76RwJD0j0KH5je1z5eCrfbYEYiq1hwd2ZpL2IXLLHL5fNJ/Za9KYtzat6AUriLlr0myfKwFf7XNAE4iI3CciKiKxBczvIyLrRWSTiDyYa3pNEZkjIhvd9xr+i9pUJA1rRBFz8f9xQsOo+b+x3B38MRv/9yFkZQY6NL/4IGW7874k1X8bzcqE9bPgm2ecdz9/1wHZ5wDz1T4HbERCEWkEvAEkAl1UdV+e+cHABpwhbVOBJcBQVV0jIs8AB1T1KTex1FDVB4rapo1IaHIb9voiFvy0nyCymB92D/XkAAqcJJwfsppzY8ZDZBFEckItJt1+QaDD9Yrsfc4WGixkZGrOezaf7XNWJvxnEOxIgdMnICwKGiTBDVMhKNj726MM7HMAeHufy+KIhM8D9wMFZbBuwCZV3ayqp4H3gQHuvAHA2+7nt4GBxdng+vXrmThxIgAZGRl4PB7effddAE6cOIHH42HKFOfBssOHD+PxePj4Y2c49n379uHxeJgxYwYAu3fvxuPxMGvWLAC2b9+Ox+Nh7ty5AGzevBmPx8M333yTs22Px8PChQsBWLVqFR6PhyVLlgCwfPlyPB4Py5cvB2DJkiV4PB5WrVoFwMKFC/F4PKxfvx6Ab775Bo/Hw+bNmwGYO3cuHo+H7dudvzRmzZqFx+Nh9+7dAMyYMQOPx8O+fU6e/vjjj/F4PBw+fBiAKVOm4PF4OHHiBADvvvsuHo+HjAznnOnEiRPxeDw53+Xrr7/O5ZdfntN+5ZVXuOKKK3LaL774Iv37989p/+1vf2Pw4ME57aeeeoohQ4bktJ944gmuv/76nPbYsWO5+eabc9oPPfQQI0aMyGnfd9993HXXXTnt0aNHM3r06Jz2XXfdxX333ZfTHjFiBA899FBO++abbyZsxYdEhgbjCVrOvZ/s4slvTxEkEC2n+PvUFIYsuZF5YX/kLztv4Yo2VRnzmzpsf7wN2x9vQ6/EGB67oi6pj7cm9fHWXNqiCuP71WXH44nseDyRixKieaZ/XXaMa8WOca3oHh/F8wPrsnNcK7aNbUmP+Cj+flU9do1ryU9jmtMjPopXr6nHrnEtWP9QM3rER/LmtfXYPa4Fq+932u9cV5+0cc1ZeV88PeIjmXx9ffaMS+CHe5uQHB/Jf2+sx55xCXw/qjHJ8ZF8cnN99o5rxnd3NyQ5PpLPbq3PCzuGMvHw7dR/fzCv7ruN+SH/x5sHbyNu0lW8efA2vg+/k3/tu5Wt4z3Mv7M+B8Y1YdatdbkwPoLFd9fn4LjGfHpzbS6Mj+CHUU576o1Oe+W99Tk0rhEfXh/HhfHhrL+vPofHNeT962pxYXw4W+6vx4nHG/DRzLl4XtvD4fQsOH2cd2Z8Q5PmLfH8dRY9n/2atteNoUZCRy59Zi49n/2aNkMepEZCR3o++zU9n/2axGvuo2bLLjntVleNplbi+TntlgNHEtc2Oae9YOq/2Tv1Lzn/7fct+IC9057O+SE9tGAy+2b8jS37jtPz2a9p2utG6nW9Imf5+J7DaHBBv5x2Y8+1NEwemNNudPHVNLr46px2w+SBNPZcm9NucEE/4nsOy2nX63oFTXvdmNOu26UXzfrcktOu0+kyEq4ckdOu3eFiWvz2/3LacW2TaTlwZE67VuL5tLpqdE67ZssuLPr8A7IHbt793oMcXD4HgNOnM9j93oMcW/0VkaHB3Na9QbF+9woSUuAcHxKR/sAOVf1RCh6fugGwPVc7FTjf/VxHVXcBqOouEaldyLZGACMAwsPDSxu6qUAa1Yji9uFJpLz9EYvJIvffU4JyijDSoloQHhLEqeBDHA+tQVpUCwBOBR/kaFgsu6IS3PYBjoTVYWdUM6cddIDD4fXYUaUpAKeDD3AwoiGpVZpwJjOLU8EH2R/RkJ+rNCb99BnSgw+zL6IxP1dpxLHgDNKDj7I3silbqzTgkJwiPfgoaRHN2FylPvuz0kkPXsTuyOb8VKUuaWdOcDLke3ZGtWJTTB12nD7GyZAl7IhKZGNMbbadPMKJkBS2R7WletVYfjp8gEO6lJSslhzMrMHWzMMc0pUszkxkb1Y19kcc50TICjZHd4KY6myO2sfxkOVsik7idExVtkTt4XjIj2yI6caxqjFsid7F8ZCVrK9yAQerRrM5aifHQlazOiaZ3VUj2RKVyrGQNayIuZgOIVuBlWf9dwglg6rBp2hZN4bQ8EikeiSHwkNoVTeG4JBQqBbB4fAQEutVBSCrWgTHw35pn6kWQXqu+aeqRpAR8Uv7ZI1IzkSHESyQmc+fq0ECsVXC6NTEORN+JCacY6dDc5Y/GB3GKQ3Lae+PDiPztOa090WHAeS090SFERz2S//dUWGER//S3hkVSpWY8Jz29shQauRq/xwRSmyu9taIUOKq/tLeHBFC7aoROe1N4SHUrfZLe0NYCI1rRlG7WU1SthzI999+WHAQb93UlfPqReY7v7h8dgpLROYCdfOZNQZ4GOitqodFZCuQlM8prGuA36jqbW77BqCbqt4tIodUtXquvgdVtcjrIHYKy+Rn+bz3afHtKKLll7pYxzWcjRe/RMeeQwpZsvyatzaNOyct49SZX57CDw8J4pVhnenZuo7vNrx+Fnx0C5w+/su0sGgY/Ba06uO77RLAfQ4gb+2z309hqerlqtou7wvYDDQFfnSTR0NgmYjkTTapQKNc7YbATvdzmojUA3Df9/hqP0zFt6VGd1bQnOMaTqYKxzWcFbRgS43ugQ7NZ46kZxASJAQJRIQGESQQHCS+v7W1RS/nmkdYNCDOe4MkZ7qPBWyfA8jX++z3ayCqulJVa6tqvKrG4ySKzqq6O0/XJUALEWkqImHAEGC6O286MNz9PByY5ofQTQU1JWUn159+kAkxD7D9vNFMiHmA608/yAcpuwIdms9MWbKdExmZtK5XlddvTKJ1vaqczMj0/Z1JQcHOBfPBb8GlY5x3H15Azy1g+xxAvt7ngFwDKYiI1AfeUNW+qnpGREYCXwDBwFuqutrt+hTwgYjcCvwMXBOYiE1FEBMeyoN923Jrcj+CgoQ/Zyn1529mydaDgQ7NZ2LCQ3m4b2tuTW5KUJDQIyGWN/21z0HBzukqH5+yyiug+xwgvt7ngN3GGwh2DcQYY85dWbyN1xhjTDlmCcQYY0yJWAIxxhhTIpZAjDHGlEiluoguInuBbSVcPBbYV2SvisX2uXKwfa4cSrPPTVQ1Lu/ESpVASkNEUvK7C6Eis32uHGyfKwdf7LOdwjLGGFMilkCMMcaUiCWQ4nst0AEEgO1z5WD7XDl4fZ/tGogxxpgSsSMQY4wxJWIJxBhjTIlYAimCiDQSka9EZK2IrBaRUYGOyR9EJFhEfhCRTwMdiz+ISHUR+VBE1rn/rSvuYCAuEfmD+296lYhMFpGIQMfkbSLylojsEZFVuabVFJE5IrLRfS9yMLrypIB9nuD+214hIlNFpLo3tmUJpGhngD+qamvgAuAuEWkT4Jj8YRSwNtBB+NGLwCxVTQTOo4Lvu4g0AO7BGQ20Hc6QCRVx+MWJQN668Q8C81S1BTDPbVckE/n1Ps8B2qlqB2AD8JA3NmQJpAiquktVl7mfj+L8sDQIbFS+JSINgSuBNwIdiz+ISFXgYuBNAFU9raqHAhqUf4QAkSISAkTxy4ifFYaqfgvkHRh8APC2+/ltYKA/Y/K1/PZZVWer6hm3uQhnhNdSswRyDkQkHugELA5wKL72AnA/kFVEv4qiGbAX+Ld72u4NEYkOdFC+pKo7gL/hDMi2CzisqrMDG5Xf1FHVXeD8gQjUDnA8/nYL8Lk3VmQJpJhEpArwETBaVY8EOh5fEZF+wB5VXRroWPwoBOgM/FNVOwHHqXinNc7invcfADQF6gPRInJ9YKMyviYiY3BOy0/yxvosgRSDiITiJI9JqvpxoOPxsWSgv4hsBd4HLhORdwMbks+lAqmqmn1k+SFOQqnILge2qOpeVc0APgZ6BDgmf0kTkXoA7vueAMfjFyIyHOgHDFMvPQBoCaQIIiI458bXqupzgY7H11T1IVVtqKrxOBdVv1TVCv2XqaruBraLSCt3Uk9gTQBD8oefgQtEJMr9N96TCn7jQC7TgeHu5+HAtADG4hci0gd4AOivqie8tV5LIEVLBm7A+Ut8ufvqG+igjNfdDUwSkRVAR+AvgQ3Ht9yjrQ+BZcBKnN+CClfeQ0QmA98BrUQkVURuBZ4CeonIRqCX264wCtjnl4EYYI77G/aqV7ZlpUyMMcaUhB2BGGOMKRFLIMYYY0rEEogxxpgSsQRijDGmRCyBGGOMKRFLIKbME5GvReQ3eaaNFpFXzmEdj4vI5UX0uUlE6udqv+HtwpnuNl725jp9ya1SfGeg4zBlkyUQUx5M5teVYoe404skIsGqOlZV5xbR9Sacsh4AqOptqlrmHygUh6/+X64OnFMC8XE8pgyx/8imPPgQ6Cci4ZBT1LI+MF9E/ikiKe64FuOyFxCRrSIyVkTmA9eIyEQRudqdN1ZElrjjYLzm/uBdDSThPEy4XEQi3SOfJHeZoSKy0l3m6VzbOSYi40XkRxFZJCJ13Om/FZHFbnHGudnTCyIij4nIf0TkS3ecitvd6VVEZJ6ILHO3PyD7OxBn3JJXcB4GbFTEd/EXEfnOnd9ZRL4QkZ9E5I5c/f7kfi8rci3/FJDgficTCupXQDwT3e9rpYj84Rz/m5vyQFXtZa8y/wI+Awa4nx8EJrifa7rvwcDXQAe3vRW4P9fyE4Grcy/jfv4P8Fv389c442OQu42TrH4G4nAKL34JDHT7aK7lnwEecT/X4JcHdW8DnnU/3wS8nM/+PQb8CEQCscB2d7shQFW3TyywCRAgHqda8gW51lHYd/F/7ufngRU4TyXH4RTOBOiN8yS64Pxh+SlOift4YFWubRTWLyceoAswJ9dy1QP9b8he3n/ZEYgpL3Kfxsp9+up3IrIM+AFoC+S+ZjGlgHVd6h4drAQuc5crTFfga3UKD2ZXMr3YnXca50cUYCnODyk44y184W7jT8XYBsA0VT2pqvuAr4BuOD/Uf3FLrMzFGYsm+2hmm6ouyrV8Yd/FdPd9JbBYVY+q6l4gXZzR6Xq7rx9wjiASgRb5xFhYv9zxbAaaicjfxanDVGErWFdmIYEOwJhi+gR4TkQ6A5GqukxEmgL3AV1V9aCITARyD8t6PO9KxBm29RWcI43tIvJYnmXyI4XMy1DV7HpAmfzy/9TfgedUdbqIeHCOMIqSt66QAsNwjhS6qGqGOFWSs+PN2b9ifBen3PesXJ+z2yE4+/hXVf1X7gDc04VnTSqkX048bgznAb8B7gJ+hzMOhalA7AjElAuqegzntMxb/HL0URXnR+uwe43himKsKvtHdZ84Y7xcnWveUZxTO3ktBi4RkVgRCQaGAt8UsZ1qwA738/DCOuYyQEQiRKQW4AGWuOvZ4yaPS4EmBSxbku8ity+AW9zvBBFpICK1+fV3UlC/s4hILBCkqh8Bf6bil8evlOwIxJQnk3HGrRgCoKo/isgPwGqcUyYLilqBqh4SkddxTuVsxfmRzjYReFVETgLdcy2zS0QewjmtJMBMVS2qBPhjwH9FZAfOEKJNi7F/3+Nc62kMPKGqO0VkEjBDRFKA5cC6AvbrnL+LPMvPFpHWwHciAnAMuF5VfxKRBSKyCvhcVf+UXz+co6/cGuCM8Jj9R6pXxuA2ZYtV4zWmDHBPpR1T1b8FOhZjistOYRljjCkROwIxxhhTInYEYowxpkQsgRhjjCkRSyDGGGNKxBKIMcaYErEEYowxpkT+H1dIeMjio+/YAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot the results\n",
    "plt.plot(param_QAOA,energ_QAOA, marker=\"*\", markersize=10,label=\"QAOA\")\n",
    "plt.plot(param_VQE,energ_VQE,marker=\".\",markersize=10,label=\"VQE\")\n",
    "plt.axhline(exact(HP(param,graph)), color=\"black\", linestyle=\":\", label=\"exact\")\n",
    "plt.xlabel(\"Variational parameters\")\n",
    "plt.ylabel(\"Objective\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Some questions and suggestions\n",
    "\n",
    "- VQE vs QAOA not clear which one is better (this is a very simple example). Explore the advantages of QAOA (precision of the results, number of variational parameters, circuit depth, etc).\n",
    "- Almost all QAOA tutorials use the \"square\" MAX-CUT problem. Other problems where combinatorial optimization is required? (e.g. Travelling Salesman Problem)\n",
    "- Write a beautiful Tequila tutorial about QAOA, upload it to the [Tequila tutorials Github](https://github.com/aspuru-guzik-group/tequila-tutorials) and become a Tequila contributor!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:tqenv2] *",
   "language": "python",
   "name": "conda-env-tqenv2-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}