tsp_qaoa.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "entitled-aquarium",
   "metadata": {},
   "source": [
    "## This notebook is part of self learning Quantum Computing Algorithms"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c538bacd",
   "metadata": {},
   "source": [
    "## Solve TSP using QAOA "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "6099ddc8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# useful additional packages\n",
    "import networkx as nx\n",
    "import numpy as np\n",
    "\n",
    "from qiskit import Aer\n",
    "from qiskit_optimization.applications import Tsp\n",
    "from qiskit.algorithms import QAOA, NumPyMinimumEigensolver\n",
    "from qiskit.algorithms.optimizers import SPSA\n",
    "from qiskit.utils import algorithm_globals, QuantumInstance\n",
    "from qiskit_optimization.algorithms import MinimumEigenOptimizer\n",
    "from qiskit_optimization.problems import QuadraticProgram\n",
    "from qiskit_optimization.converters import QuadraticProgramToQubo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "c146598c",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "distance\n",
      " [[ 0. 48. 91. 33.]\n",
      " [48.  0. 63. 71.]\n",
      " [91. 63.  0. 92.]\n",
      " [33. 71. 92.  0.]]\n",
      "\n",
      "Exact:\n",
      " optimal function value: 236.0\n",
      "optimal value: [0. 0. 0. 1. 1. 0. 0. 0. 0. 1. 0. 0. 0. 0. 1. 0.]\n",
      "status: SUCCESS\n",
      "\n",
      "Route:\n",
      " [1, 2, 3, 0]\n",
      "\n",
      "QAOA:\n",
      " optimal function value: 236.0\n",
      "optimal value: [0. 1. 0. 0. 0. 0. 1. 0. 0. 0. 0. 1. 1. 0. 0. 0.]\n",
      "status: SUCCESS\n",
      "\n",
      "Route:\n",
      " [1, 2, 3, 0]\n"
     ]
    }
   ],
   "source": [
    "# Generating a TSP instance of n cities\n",
    "n = 4\n",
    "num_qubits = n ** 2\n",
    "ins = Tsp.create_random_instance(n, seed=123)\n",
    "adj_matrix = nx.to_numpy_matrix(ins.graph)\n",
    "print('distance\\n', adj_matrix)\n",
    "\n",
    "# Create a random TSP Hamiltonian\n",
    "qp = ins.to_quadratic_program()\n",
    "\n",
    "\n",
    "algorithm_globals.random_seed = 10598\n",
    "quantum_instance = QuantumInstance(Aer.get_backend('qasm_simulator'), seed_simulator=algorithm_globals.random_seed, seed_transpiler=algorithm_globals.random_seed)\n",
    "qaoa_mes = QAOA(quantum_instance=quantum_instance, initial_point=[0., 0.])\n",
    "exact_mes = NumPyMinimumEigensolver()\n",
    "\n",
    "def route_x(x):\n",
    "    n = int(np.sqrt(len(x)))\n",
    "    route = []\n",
    "    for p__ in range(n):\n",
    "        for i in range(n):\n",
    "            if x[i * n + p__]:\n",
    "                route.append(i)\n",
    "    return route\n",
    "            \n",
    "# solving Quadratic Program using exact classical eigensolver\n",
    "exact = MinimumEigenOptimizer(exact_mes)\n",
    "exact_result = exact.solve(qp)\n",
    "print(\"\\nExact:\\n\", exact_result)\n",
    "print(\"\\nRoute:\\n\", route_x(exact_result.x))\n",
    "\n",
    "# solving the Problem using QAOA\n",
    "qaoa = MinimumEigenOptimizer(qaoa_mes) \n",
    "qaoa_result = qaoa.solve(qp)\n",
    "print(\"\\nQAOA:\\n\", qaoa_result)\n",
    "qaoa_result = np.asarray([int(y) for y in reversed(list(qaoa_result))])\n",
    "print(\"\\nRoute:\\n\", route_x(qaoa_result))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8583779b",
   "metadata": {},
   "source": [
    "## Weight Matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "72cd1ff7",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "array([[68.23518632,  5.38210188],\n",
    "       [22.03598728, 18.43718107],\n",
    "       [17.59059011, 81.20945067],\n",
    "       [92.3344998 , 27.65743978]])\n",
    "ref: https://numpy.org/doc/stable/reference/random/generator.html\n",
    "\"\"\"\n",
    "coordinates = np.random.default_rng(123).uniform(low=0, high=100, size=(4, 2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "57787e31",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "{0: (68.23518632481435, 5.382101880222267), 1: (22.03598727726114, 18.43718106986697), \n",
    "2: (17.590590108503037, 81.20945066557736), 3: (92.33449980270564, 27.657439779710625)}\n",
    "\"\"\"\n",
    "pos = dict()\n",
    "for i, coordinate in enumerate(coordinates):\n",
    "    pos[i] = (coordinate[0], coordinate[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "7b6f0980",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "ref: https://networkx.org/documentation/latest/reference/generated/networkx.generators.geometric.random_geometric_graph.html\n",
    "\"\"\"\n",
    "high = 100\n",
    "low = 0\n",
    "graph = nx.random_geometric_graph(n=4, radius=np.sqrt((high - low)**2 + (high - low)**2) + 1, pos=pos)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "b1207a00",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "graph.edges = EdgeView([(0, 1), (0, 2), (0, 3), (1, 2), (1, 3), (2, 3)])\n",
    "graph.nodes = NodeView((0, 1, 2, 3))\n",
    "ex: graph.nodes[0] = {'pos': (68.23518632481435, 5.382101880222267)}\n",
    "rint (Round elements of the array to the nearest integer)\n",
    "ref : https://numpy.org/doc/stable/reference/generated/numpy.rint.html\n",
    "weights or edges = sqrt((x1 - x2)**2 + (y1 -y2)**2)\n",
    "ref: https://en.wikipedia.org/wiki/Euclidean_distance\n",
    "\"\"\"\n",
    "for w, v in graph.edges:\n",
    "    delta = []\n",
    "    for i in range(2):\n",
    "        delta.append(graph.nodes[w][\"pos\"][i] - graph.nodes[v][\"pos\"][i])\n",
    "    graph.edges[w, v][\"weight\"] = np.rint(np.sqrt(delta[0]**2 + delta[1]**2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "00a57d1b",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0. 48. 91. 33.]\n",
      " [48.  0. 63. 71.]\n",
      " [91. 63.  0. 92.]\n",
      " [33. 71. 92.  0.]]\n"
     ]
    }
   ],
   "source": [
    "\"\"\"\n",
    "graph.adjacency() = holds weights of the graph\n",
    "ex: 0 {1: {'weight': 48.0}, 2: {'weight': 91.0}, 3: {'weight': 33.0}}\n",
    "\n",
    "np.full((4, 4), np.nan) = \n",
    "array([[nan, nan, nan, nan],\n",
    "       [nan, nan, nan, nan],\n",
    "       [nan, nan, nan, nan],\n",
    "       [nan, nan, nan, nan]])\n",
    "index = dict(zip([0, 1, 2, 3], [0, 1, 2, 3]))\n",
    "\"\"\"\n",
    "index = dict(zip(list(graph), range(n)))\n",
    "A = np.full((n, n), np.nan)\n",
    "for u, wdict in graph.adjacency():\n",
    "    for v, d in wdict.items():\n",
    "        A[index[u], index[v]] = d.get(\"weight\", 1)\n",
    "\n",
    "A[np.isnan(A)] = 0.0\n",
    "A = np.asarray(A)\n",
    "M = np.asmatrix(A)\n",
    "print(M)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "d12ee2e4",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "3edda678",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAADnCAYAAAC9roUQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAABBUUlEQVR4nO3dd1zVZf/H8ddZ7C0IoqKiYlqag9ttmrnKxJUjU8xUcJOaMzNH7pEjcaa34MqJYN6ljZ8NM3dDU1Fyo4Age51zvr8/vh0SB6IezsDr+XjwMOBwzucQvLnO53sNhSRJCIIgCKahNHcBgiAIzxMRuoIgCCYkQlcQBMGEROgKgiCYkAhdQRAEE1IX9UlPT0+pcuXKJipFEAShdDhx4kSSJEleD/tckaFbuXJljh8/XjJVCYIglFIKheLKoz4n2guCIAgmJEJXEATBhEToCoIgmJAIXUEQBBMq8kKaSWi1EBcH587BrVuQnw/OzlCtGrzwAnh4mLtCQRAEozFf6F69Cps3y2/5+SBJkJcn/6tSgY2NHMgvvQShofDaa/LHBEEQrJjpQzc3F5YuhbVrQa8HR0d5ZPswkgR//QWjRkH58rBsGdSpY9p6BUEQjMi0Pd3r1+H112H1ajls3d2LHr0qFHIgu7rCzZvQrRuEh8thLAiCYIVMN9K9cUMOzZQUOWyflIuL3G5YtAiys2HsWOPXKAiCUMJMM9LNy4PgYDlwXV2f/n7Uajl8w8MhOtp49QmCIJiIaUJ36VK4fPnZAtdApQJ7e/jwQ0hIePb7EwRBMKGSD91r1+SLZo+6WPY07OwgKwvmzzfefQqCIJhAyYfu1q2g08mtAWNycYGYGLllIQiCYCVKNnT1eoiMlGcqFCFPr2fGzZu8efEir5w/T5+4OA5nZBR93yqVHOb79xuxYEEQhJJVsqF7+bI8L/cxixp0gI9Gwxo/P/4vIIChXl5MvHGDm3l5Rd+/QgE//mi0cgVBEEpayU4ZO3dODsbHsFcqCfH6d7/fFs7O+Go0nMvJwbeowLazg9OnjVCoIAiCaZRs6N6+LU8XewJarZY04GpeHv62tkXfWKOBxMSnr08QBMHESjZ0dbonWj2Wm5vLH2fP8qlaTVsXFyo9LnQVCvkxBEEQrETJ9nQdHeULXsWlVLJWpUKh09EzJ4e4uDjupqaif1Rw63Tg4GCcWgVBEEygZEO3enW5BVAMkiQxJzGRFEnis4AANAoFjo6OpKamcvHiRZLu3EGn1xf+otxcqFGjBAoXBEEoGSUbui+8IO+XUIwWw5xbt7icl8cUe3scNRoqVapEeloaLi4u+FWsSG5ODhcvXuR2QgL5Wq38RTk50KhRiT4FQRAEYyrZnq6TEzRsCMeOFbkEOD4/n91372KjUDAwNxeb2FiUKhXjPD2pmZSEh4cH5cuXJy8/n+TkZOLi4nB2dMTT1hab9u1L9CkIgiAYU8nvMjZokBy6kvTI6WPlNBqO16wJwPXr13F2ccHVxQWAPBcXrl69iqTX4+npiY+3N16enqTfusWZtDT+u2EDwXo99evXR1GM6WmCIAjmVPLLgFu2BH9/SEsr1s3VGg1aQ/sAsPmn1ZCalkZCYiISoFIocHNw4KWICFq2asWsWbPo378/33zzDfr7+76CIAgWpORDV6WSdxkD+Viex9BoNOTfdzuNWk2lSpXIyMjg9u3bSKmp0K4dmo4d6datGzt37mTAgAFs2bKFbt26sWPHDnJyckri2QiCIDwT02zt+MILMGUKZGTIF9aKoFGr0T4knNUqFZX8/CA1les6HfqZMwvaFUqlkldffZX169czffp0fvnlF4KCglizZg13794tiWckCILwVEx3XE///jBuHKSnyyc/PIJao/l3dsK99HpUaWl41a7NsmbNmLZsGbqHLIx4+eWXWbx4MatXr+b27dt069aN+fPnc+PGDWM+G0EQhKdi2jPShgyRz0fTaOQtGR8Srg+MdCVJHiGnpUGHDqhiYpj5+eckJyczadKkB1oRBlWqVOGjjz5i+/btODg4EBwczKRJkzh79mxJPTtBEITHUkhFzKENDAyUjh8/bvxHTU6GTz+FnTv/ncdrZwcaDRIQe+EC1StWRJGbK+/D6+8P48fLx7D/Iy8vj0mTJqHT6Zg/fz42j9nJLCsri6ioKDZv3kzFihUJDg6mSZMmYsaDIAhGp1AoTkiSFPjQz5kldA3S0+HAAXl7xpMn5Q1y9Hr+vnWLcq++il3LltC2LdSu/dDpZlqtlqlTp3L37l0WLVqEvb39Yx9Sq9Vy4MABIiMj0ev1BAcH065dOzTFXDknCILwOJYbuo/w7rvvMmbMGOrUqfPY2+r1ej755BOuXLnC0qVLcXJyKtZjSJLEkSNHiIyM5PLly/Tp04euXbvi+JgN1wVBEB6nqNA1bU+3mMqWLUtCMQ+dVCqVTJkyhRo1ajBs2DDSijkfWKFQ0KRJE8LDw1m8eDFnz56lU6dOLFu2jESxXaQgCCXEIkPX29ub27dvF/v2SqWScePGERgYSEhICMnJyU/0eC+88AKzZ89m06ZN5Obm0rNnT2bMmEFcXNyTli4IglAkiwzdsmXLPlHogjxyHTlyJK1bt2bw4MHFHinfy9fXl3HjxhEVFYWvry9Dhgxh9OjRnDp1iqLaMIIgCMVlkaHr7e39VKGpUCgICQmhc+fODB48mJs3bz7V47u6ujJo0CBiYmJo3rw5M2bMYMCAAXz33XdimbEgCM/EYkP3SUe69woODuadd95h8ODBXL169anvx9bWlu7du7Nr1y6Cg4OJiIgoeD83N/ep71cQhOeXRYbu07QX7tezZ09CQ0MJDQ3l4sWLz3RfSqWS1q1bs2HDBqZOncpPP/1Ep06dWLduHampqc9034IgPF9KfmvHp+Dl5UVycjI6nQ7Vkxz3c5+goCDs7OwYNmwYS5YsoVatWs9Ul0KhoF69etSrV4+4uDg2bdpE165def3113nnnXfw9fV9pvsXBKH0s8iRrlqtxs3NjTt37jzzfbVr144PP/yQsLAwfvvtNyNUJ/P392fq1Kl88cUX2NnZ0bdvXyZPnsy5c+eM9hiCIJQ+Fhm6YJwWg0HLli2ZMWMGY8eO5dixY0a5TwMvLy9GjhxJTEwMtWrVYsyYMQwbNoxffvlFzHgQBOEBFhu6Pj4+TzWD4VGaNGnC/PnzmTx5Mj///LPR7tfA0dGRvn37snfvXjp27MiSJUvo06cP+/fvL7QpuyAIzzeLDV1jjnQN6tevz+LFi5k+fTrfffedUe/bQKPR0LFjR7Zt28bIkSOJjo6mc+fObNmyhaysrBJ5TEEQrMdzFboAtWvXZvny5cybN4///e9/Rr9/A4VCQdOmTVm1ahULFizgjz/+oFOnTqxYscIovWpBEKyTxYausdsL96pRowYrV65k+fLl7Nmzp0Qe4161atVizpw5REREkJmZSY8ePfjkk0+4fPlyiT+2IAiWxWJDt6RGugb+/v6sWrWKzz//nK1bt5bY49yrfPnyjB8/nt27d1O2bFlCQkIYO3asUWdVCIJg2Z7b0AXw8/Nj7dq1fPHFF2zYsKFEH+tebm5uhISEEB0dTePGjZk6dSrvvfce33//vVhmLAilnMWGrpeXF3fu3CnxECpXrhxr165l//79rFy50qTTvOzs7OjRowd79uyhT58+bNiwgbfeeovdu3eTl5dnsjoEQTAdiw1djUaDq6srSUlJJf5YXl5erFmzhh9//JFPP/3U5PNrlUolbdq0YePGjXz44YccOnSITp06sX79+mLvDywIgnWw2NCFp99t7Gm4u7uzatUqfv/9d+bMmWOWl/kKhYIGDRqwdOlSwsPDuXr1Kl26dGHRokXEx8ebvB5BEIzPokP3SU6QMAYXFxfCw8O5fPky06ZNe+gR76ZStWpVpk2bxrZt29BoNLzzzjtMmTKFCxcumK0mQRCenUWHrre3N7du3TLpYzo4OLBs2TKSk5OZPHnyI494N5WyZcsyatQooqOjCQgIICwsjOHDh/Prr7+KZcaCYIUsPnRNOdI1sLOzY/HixeTn5zNu3DiLuKjl5OREcHAw0dHRdOjQgUWLFvHOO+/w1VdfmXVELgjCk7Ho0DV1e+FeNjY2zJ8/HwcHB95//32ys7PNUsf9NBoNnTp1Ytu2bQwbNozdu3fTpUsXtm7dKpYZC4IVsOjQfdYTJJ6VWq3mk08+wdvbm5EjR5KZmWm2Wu6nVCpp3rw5a9asYe7cuZw+fZqgoCDCw8Of+GBOQRBMR4TuYyiVSj766COqV6/O0KFDLXIK14svvsi8efPYsGEDaWlpdO/enVmzZnHlyhVzlyYIwn0sOnQ9PT1JSkoy+yotpVLJ+PHjadCgAaGhoRY7kqxYsSITJ05k165deHp6MmjQID744AN+//13c5cmCMI/LDp0bWxscHFxsYiQUygUjBo1ildfffWpj3g3FQ8PD0JDQ4mOjqZhw4ZMmTKFgQMH8sMPP5j9D5ggPO8sOnTBMloMBsY64t1U7O3t6dmzJ3v27KF3796sXbuWnj17snfvXouYkSEIzyOLD11TbHzzpIx1xLupqFQq2rZtS0REBBMnTuTbb78lKCiooAcsCILpWORpwPcy11zdx+nZsyd2dnaEhoby2WefUbVqVXOX9FgKhYLAwEACAwOJjY0lMjKSLl268Oabb9KnTx98fHzMXaIglHoWP9K1pPbC/YKCgnj//fcZOnSo1Z0CXL16dWbMmMHWrVtRKBT06dOHqVOnEhsba+7SBKFUs/jQtcT2wr3at2/P5MmTGTVqlFXOEvD29mb06NHs3bsXf39/Ro4cyciRIzl69KhYZiwIJcDiQ9dS2wv3atWqFdOnTy+RI95NxdnZmXfffZfo6Gjatm3L/Pnz6devHwcOHBDLjAXBiKwidC15pGvQpEkT5s6dy+TJkzl8+LC5y3lqNjY2BAUFsX37dkJCQtixYwddu3bliy++sJil0IJgzSw+dL28vEhMTLSK+aUNGjRg8eLFTJs2rcSOeDcVpVLJK6+8wtq1a5k1axYnTpygU6dOrFq1yiLmTQuCtbL40LWxscHZ2ZmUlBRzl1IstWvXZtmyZcybN4+vvvrK3OUYRe3atZk/fz7r168nOTmZ7t27M2fOHK5du2bu0gTB6lh86IL1tBgMXnjhBVauXMnSpUuJiooydzlG4+fnx+TJk9m1axdubm4MGDCA8ePH8+eff5q7NEGwGlYRupY+g+Fh/P39Wb16NevWrWPbtm3mLseoPDw8GDp0KNHR0dSvX59JkyYREhLCjz/+aBVtIEEwJ6sIXWuYwfAwhiPet23bxn//+19zl2N0Dg4O9O7dm6ioKLp3786qVavo1asX0dHRYpmxIDyC1YSutY10DcqVK8eaNWvYt28fq1atKpVzX1UqFe3bt2fTpk2MGzeOAwcO0LlzZzZu3Eh6erq5yxMEi2IVoWuN7YV7lS1bljVr1vDDDz+wZMmSUhm8IC8zbtiwIZ999hlLly7l4sWLdO7cmSVLlljlKxVBKAlWEbrW2l64l4eHB6tWreLUqVPMnTu31Pc+AwICmDlzJlu2bEGv19O7d28+/vhjLl68aO7SBMGsrCZ0rXmka+Di4sLKlSuJi4tjxowZz8VKLx8fH8aMGUNUVBSVK1dm+PDhjBo1iuPHj5faEb8gFMUqQrds2bJWs0DicRwdHVm2bBmJiYlMmTIFrVZr7pJMwsXFhQEDBhATE0Pr1q2ZPXs2/fv35+DBg0X+8YmPj+fIkSMmrFQQSpZVhK6NjQ2Ojo7cvXvX3KUYhb29PZ9++im5ubmMHz/+ubrSb2NjQ5cuXdi5cycDBw5k69atnDhx4qF/UP/44w/69OnDJ598QvXq1dm/f78ZKhYE47KK0IXS02IwMBzxbmtra1FHvJuKUqmkZcuWrF+/nsDAQJTKwj+K2dnZfPHFFzRr1ox9+/axZs0aoqOjAfj++++5deuWOcoWhGcmQteM1Go1s2bNomzZshZ3xLsp3R+4ho/t27ePbt26AZCTk4ODgwPXrl1j9+7dtGzZktDQULEUWbA6VhO6H3/8Mc2bNzd3GUanVCqZOnUq1apVs9gj3s1BoVDQsmVLfv75Z+7cucOWLVtQKpVoNBqWL1/O+fPn0ev1bN682dylCsITsZrQdXV1Ra22+NOFnopSqWTChAnUr1/foo94NyUbGxtCQ0P59ddfGTp0KDdv3iQgIAAfH5+CC28dOnTg77//FgswBKtiNaFb2ikUCsLCwmjVqhUhISFWPy/5Wdy4cYPz589Tq1YtNm3aRKNGjXjppZdo0aIFIK+A++WXXzhy5Aiurq44Ozs/cB9iOppgqUrn0NFKKRQKQkNDsbOzIyQkhPDwcHx9fc1dlsn99NNPTJkyhZYtW1KmTBlOnTrFokWLqFmzJufPn2fbtm189dVX9O3bly5dugCg0+lQqVQAxMXFsWrVKvbv309QUBDTpk3DxsbGjM9IEP4lRroWqH///rz99tuEhIRYxRHvxtarVy+OHz+On58fDRo0ICYmBoARI0Ywe/ZskpKSWL16NcOHD6d8+fIABYF7/fp1xo8fj1qt5vvvv+f8+fN8+eWXZnsugnA/qxzpSpJEXl4etra25i6lxPTq1avgiPcVK1bg7+9v7pJMytXVlalTpwLy9LH+/fvz22+/ERkZSZ8+fQD550ChUBT8m5WVxddff42zszMff/wxtra2uLq6FrRqDLcTBHOyqpHuyZMn+fTTT+nRowcjR46kVatWhIaGcufOHXOXViI6d+5MWFgYQ4YMsboj3o3J3t6ekydPsm/fPtavX0+rVq04dOhQQYAa/r127Rq//fYbbdu2xdbWloSEBKpXr17QWhCBK1gCqwndkydPMm/ePE6dOgXA5cuXGTFiBE2aNGHIkCFcunTJzBWWjA4dOjBp0iSrPeLdmF5//XW++eYb5s6dW/Aq59SpU+Tk5ACQnp7OtWvXeP311wG4dOkS169fx8fHBxAX1wTLYDWhO23aNN555x0iIiLYuXMn1apVo3z58rz77rtotdpSHUivvvoq06ZNY+zYsZw4ccLc5Zhd48aNady4MQD79+/nhx9+AODnn39GqVTi7u5OUlISp06dQq1W07p1a3OWKwiFWE3oVqxYkfPnz5Ofn8/NmzeJj4/n77//BuQTg0v7yqSmTZsyZ84cJk6caNVHvBvbhx9+SLt27QBo0qQJOp2OGzdusHjxYv788086deqEra1twXWAiIgIUlNTzVy18DyzmtB98803OX/+PNOmTWPatGnUq1ePN954A4BRo0bx5ptvmrnCkhcYGMiiRYuYNm0a33//vbnLsTiBgYHUrVuXzp07c/PmTT744APatGkDyP3c3Nxcrly5QpcuXZg/fz43b940c8XC80hRVJ8rMDBQOn78uAnLKdoPP/zAzz//jKenJ126dMHLywuAI0eO8PvvvxMSEmLmCk3j3LlzjBo1ijFjxtChQwdzl2ORMjIycHJyeuiMhaSkJLZt28aePXto1KgR/fr1o2bNmmaqVCiNFArFCUmSAh/6SUmSHvnWoEEDyVIcPXpUCgoKKng/IyNDiomJKfhcmzZtzFWaWVy8eFHq0KGDFBUVZe5SrFZmZqa0adMm6fXXX5eGDBki/fzzz5Jerzd3WUIpAByXHpGrVtNe8PT0LLQ0Nj8/nzlz5gBQpUqV5279fdWqVVm9ejVr167liy++MHc5VsnBwYF33nmH6OhoOnXqxLJly3j77bfZv3//c7O5vGB6VhO65cqVIzs7m6SkJACOHTvG2bNn2bp1K19++SWurq5mrtD0/Pz8WLNmDVu2bGHjxo3mLsdqqdVq3njjDbZu3UpYWBgxMTEEBQWxefNmsrKyzF2eUMpYVU937NixZGRk4O3tTWJiIh06dODLL78kPT2dGTNmUL16dXOXaBYJCQkMGzaMtm3bEhISIhYBGMFff/1FZGQkv/76K126dOHtt9/G09PT3GUJVqKonq5VhS7A2rVr+euvv2jfvj3t27cnISEBGxsb3NzczF2aWSUnJzN8+HAaNWpEWFiYCF4juXnzJps2beKrr76iVatW9OvXjypVqpi7LMHClZrQvXcnKeFBaWlpjBgxglq1ajF+/PiHnsggPJ3U1FR27NjB9u3beemllwgODubll18Wf9yEhyo1ofvLL7/g6OhInTp1yMvLQ61Wo1AoxA/+PTIzMwkLC6NixYp89NFHIniNLDc3l3379hEZGYm7uzv9+vWjVatW4vssFFJU6FrVT8qtW7cKtjq0sbFBqVSKwL2Po6Mjy5cv5/bt24SHh4v9BozM1taW7t27s3v3bvr168fGjRt566232L17N7m5ueYuT7ACVjXSvdfp06eJj48nJyeH6tWr89JLL5m7JIuSn5+PVqvF3t7+kbeJi4tDqVRSuXJl0xVWykiSxKlTp4iMjOTMmTP06tWLHj164OLiYu7SBDMqNe0FkHeO2rRpE6mpqQWthdjYWG7evMmsWbMK1uELj94/9vr166xevZoDBw6gUqkYOHAgAwcONEOFpUtcXBybNm3i//7v/3jjjTfo06fPc3nyh1CKQtfQr9TpdLz11ltUqlQJZ2dnHBwcuHLlCsHBwZw9e9bcZVq8efPmceHCBWbMmIGDgwNvvvkmGzZsICAgwNyllQoJCQls27aNqKgomjRpQnBwMDVq1DB3WYIJFRW6VndyRExMDLdv337g487OzqK/Wwz5+fn88ccfDBo0iPLly5OVlYWTkxO2trbo9XpxQcgIypYty6hRoxg4cCC7d+9m9OjRVK5cmeDgYBo1aiR+Tp9zVhW6hpkL8+bNw9/fH1tbW9LT04mNjWXz5s188MEH4kiWx9BoNHh6erJz504CAgKIiIigadOmeHh4sGXLFvLz8xkwYIC5yywVHB0d6devH7179+brr79m8eLFqNVq+vXrR9u2bVGrrerXTzASq2ovgLxSaOHChdjb22NnZ4dKpcLT05OOHTtSq1Ytc5dnsWJjY/n6668ZMWIEIG+HefPmTW7cuMHQoUPp0aMHly5dok6dOixcuJAxY8aYueLSR5IkDh8+TEREBDdu3KBPnz506dIFBwcHc5cmGFmpai/UrFmTzz//nDt37pCUlISjoyNubm44OTmZuzSLZmtry9q1a8nLyyM0NJSKFSuiUCho0aIFwcHBgLxlZIsWLZ67QzBNRaFQ0KxZM5o1a8bZs2eJiIjg888/p3v37vTq1YsyZcqYu0TBBKxupAviVNendfXqVQYNGoSLiwt+fn60a9eO1157DY1GwxdffMHKlSv54IMPaN++PRqNxtzlPheuX7/O5s2b+frrr3nttdfo27cvlSpVMndZwjMqNbMXDAyr0cRFnyeXn59fsJLPYPfu3SxfvpwRI0YQFBQkAtcM7t69y/bt29mxYwd16tShf//+1KlTx9xlCU+p1IXut99+yy+//MKUKVPMXYrVOnz4ML/99hv/+c9/mDhxIu+++y49evQoOGVXMI+cnBxiYmLYtGkTZcqUoX///rRo0UIMMKxMqVkGbODs7Mz169fNXYZVa9CgAStXrqRZs2YMHTqUvn37PjRw9Xq9Gap7ftnZ2dGjRw/27NlDnz59+Pzzz+nRowdRUVHk5eWZuzzBCKwydMuWLVvoFAnhydna2nLq1CmaN2/O7t27H3k7pVLJ0aNHmTlzJjNmzDBhhc83pVJJmzZt2LhxI5MnT+b777+nU6dOrF+/nrS0NHOXJzwDqw3d27dvi81cnpFKpeLbb7/F39+/0HFHhtHttWvXWLduHWFhYTg5OXH+/Hm6dOlipmqfTwqFggYNGrB06VLCw8O5evUqXbp0YfHixdy6dcvc5QlPwSpD18HBARsbG/EX30hmzpyJs7MzP/74IyCPspKTk/nyyy/55ZdfmDlzJqNHj2bz5s1kZ2dz48YNM1f8fKpatSrTpk1j27ZtqFQq+vTpw0cffcSFCxfMXZrwBKwydOHf0a5gPF9++SWGC6e///4733//Pb1796ZNmzbodDpWrlyJj48P5cuXN3Olz7eyZcsSFhZGdHQ01atXJywsjBEjRnD06FHx6s8KWG3oent7i76ukc2dO5fAwEDy8vIYP348LVq0oG3btgAcOXKEK1eu8Nprr5m5SsHAycmJ4OBg9u7dS7t27Vi4cCF9+/bl66+/RqfTmbs84RGsOnTFSLdk5OXlUaVKFYYMGQLAN998w759+5AkiaCgIDNXJ9zPxsaGoKAgtm3bxtChQ9m1axddu3Zl27ZtZGdnm7s84T5WtwzYQLQXSo6TkxNeXl40b96cqlWrolQqqVGjBqGhobi5uYkVgRZKqVTSvHlzmjdvzp9//klERATr1q0rWGbs4eFh7hIFrDh0vb29OXnypLnLKLU+++wzdu7cSXp6Oq+//jqurq7Y29uLwLUSL730EvPnz+fq1ats3ryZt956izZt2tC3b1/8/PzMXd5zzWpDV4x0S95bb731wMe0Wu0Dy4gFy+Xn58ekSZMIDQ1l+/btDBw4kLp16xIcHEzt2rXNXd5zSfR0hWKTJImlS5cyf/58sVLNynh4eDBkyBCio6MJDAzkww8/ZNCgQfzwww/i/6WJWXXoJiQkiCkyJqRQKBgyZAgXLlxg5syZ4pfVCtnb29OrVy/27NlDz549WbNmDT179iQ6OlosMzYRqw1dBwcHVCpVoZVUQslzcnLis88+4/bt20yZMgWtVmvukoSnoFKpaNeuHZGRkUyYMIFvvvmGoKAg/vvf/4rfqRJmtaELosVgLvb29ixZsoSsrCwmTJggRkhWTKFQ8J///Idly5axfPly4uLi6Ny5M59++qn43SohVh+6YoGEedjY2LBgwQLUajVjxowhJyfH3CUJz6h69erMmDGDrVu3AvD2228zdepULl68aObKSherDl0xg8G8NBoNs2fPxsPDg1GjRpGVlWXukgQj8Pb2ZvTo0ezduxd/f39GjBjBqFGjOH78uLiGYgRWHbqivWB+KpWKadOmUblyZYYNGyY2ISpFnJ2deffdd4mOjua1115j7ty5BAcHc+DAAbHM+BlYfeiK9oL5KZVKJk2aRJ06dRgyZAgpKSnmLkkwIhsbGzp37sz27dsZPHgw27dvp2vXrmzfvl0sM34KVh26or1gORQKBaNHj6ZFixaEhISQmJho7pIEI1MqlbzyyiusW7eOTz75hGPHjhEUFMTq1avFH9onYNWhK9oLlkWhUDB06FDeeOMNBg8eTHx8vLlLEkpInTp1WLBgAZ9//jl37tyhe/fuzJ07l2vXrpm7NItn9aErFkhYngEDBtCrVy9CQkLEL2Ep5+fnx+TJk9m5cyeurq4MGDCACRMmcObMGXOXZrGsOnQdHR1RKBRkZGSYuxThPm+//TbvvfceoaGhxMXFmbscoYR5eHgwdOhQoqOjqVevHhMnTiQkJISffvpJrFy8j1WHLoCPj49oMViorl27MmLECIYOHcr58+fNXY5gAg4ODvTu3ZuoqCi6d+9OeHg4vXv3JiYmhvz8fHOXZxGsPnTFycCW7Y033mDChAmMHDmSP//809zlCCaiUqlo3749mzdv5oMPPuCrr76ic+fOREREPPevTK12a0cDMYPB8rVu3RpbW1tGjx7NvHnzqF+/vrlLEkxEoVDQsGFDGjZsyIULF4iMjCQoKIguXbrQu3dvypYta/wHvXkTzp6FM2cgMREUCvDxgRdfhFq1oCQe8wlYfeiK9oJ1aNasGbNmzWL8+PHMnDmTJk2amLskwcQCAgKYOXMm8fHxbN26ld69e/PKK6/Qr18/qlat+mx3np8PBw/C6tVy2KpUkJcHyn9ezOv1YGsLWi00bAiDB0PLlv9+3oREe0EwmYYNG7Jo0SKmTp3KoUOHzF2OYCblypVjzJgxREVF4efnx7BhwwgLC+PEiRNPNxPp/Hno1AlGjYJz58DVFZydoUwZcHeX38qUAScncHGBY8cgJAR69wYzzK4pFaErRrrW4+WXX2bp0qXMmjWLAwcOmLscwYxcXFx47733iImJoVWrVsyaNYv+/fvzzTffFH/Gw44dcuBeuvRv2BZ1qolSKd/OxQVOnoR27eC774zzhIrJ6kNXLJCwPrVq1SI8PJzFixcTExNj7nIEM7OxsaFr167s3LmT9957jy1bttCtWzd27NhR9O51W7fCxIlgZycH6ZMcIaVQgJub3IYIDYVvv33m51FcpSZ0xQIJ61KtWjVWr17NypUr2bFjh7nLESyAUqmkVatWrF+/nunTp3PkyBGCgoJYs2YNd+/eLXzjkydh6lRwdAQbm6d/UDs7udc7fDj8/fcz1V9cVh+6hgUSmZmZ5i5FeEKVKlVi3bp1REZGEhkZae5yBAvy8ssvs2jRItasWUNCQgLdunVj/vz53LhxA7KzISxMbhU8S+Aa2NnJF9hGj5YvuJUwqw9dhUIh+rpWzNfXl3Xr1hEVFcWaNWvEKxahkMqVKzNlyhS2b9+Oo6MjwcHB7OrWjdwrV+T+rbG4usKff8JXXxnvPh/B6kMXRF/X2pUtW5a1a9fy3XffsXz5chG8wgM8PT0ZPnw4MXv30jI2lvjkZK5cuUJGRgaP+mn56MYN2sfG8sr583S7dImo+1sU91Io5P7uqlUlUX4hpSZ0xbQx6+bh4cHq1as5fvw4CxYsEOv1hYdyOHsWT5WKSjVq4ObmRkJCAnFxcdxNTUV/3x/rAZ6exFStyg81arC4QgXCExL4q6j9f52c4K+/Sry3WypCV7QXSgdXV1fCw8M5f/48n3zyiQhe4UGnT0N+PgqFAldXV6r4++Pt7U1qaiqXLl7kTnIyun9+bvxtbbH5Z/GDArkVeb2o/R8UCvntjz9K9CmUitAV7YXSw3DEe3x8vDjiXXjQ0aOg/nchrQJwcnSkkp8fFSpWJCc7m4sXL3I7IYF8rZa5t27R7Nw5usfF4alW08zJqej7z8+Xg70ElZrQFe2F0sNwxHtmZqY44l0o7OZN0Gge+il7Ozt8y5fHz8+P/Px8YmNj6ZmdzaGAANZVqkRrZ2dsHjeXV6OB69dLoPB/lYrQFe2F0sfW1paFCxeiUqkYO3asOOL9OSVJEmlpaVy9epXTp09zJymJ1LQ0EpOSuHXrFtdv3ODK1atciovjwoULnPvrL65euSKf3SZJpKeloQDqOjhwOz+fncU5VqiED920+g1vQLQXSiuNRsOcOXOYPn06o0aNYsmSJTg4OJi7LOEZSJJEZmYmKSkpBW/JyckPvH/37t2Cf21tbfHw8MDd3Z2RqalUzM5GcnREY2ODvVqNSqVCrVKh+ue/U+/eJTExEQ8PD7y8vFD909fVAdcf96pJp5P3aShBpSJ0nZycCv5nOjo6mrscwYgMR7zPmTOH4cOHs2zZMpyNOT9TeCaSJJGdnf3Y8Lz34xqNBnd3d9zd3fHw8MDNzQ0PDw98fHyoWbNmwcfd3d1xc3PD5t4FEAsWQHj4Q4MxKzuba9euoVQqcfL15Xe9nuYKBXaSxK+ZmXydmsrs8uUf/6Tq1jXeN+ghSkXo3rtAwt/f39zlCEamVCqZPHkyv/76K6dPn6Zx48ZoHtHXE56dIUQfF56G95VKZaHwNASql5cXNWrUKHjf8GZra/v0xb388gOr0LRaLbcTEsjKzKSstzcuLi7c1WrZeeMGs+Pj0QPlNBrGenvzyuP+YKvV8p67JahUhC7822IQoVs6KRQKGjVqhOJJNjURAMjNzS32KDQ5ORlJkgqFp2Hk6eHhQdWqVQuNRN3d3bGzszPdk2nWTF7+m5+PpFaTnJJCUlISbm5u+FetWtBKcFerWVOp0pPdd06OvDKtdu0SKPxfpSp0xQyG0k0EriwvL6/Yo9CUlBR0Ot0Do1BDcFauXPmBkai9vb3lfq8dHaFHD3I+/5wb2dlo1GoqV66MrTH2YDDs6aBSPft9FaHUhK6YwSAUh06nY/v27Zw8eZLmzZvz5ptvoirhX7LHyc/PJyUl5bHhaXg/Ly/vgfA0/HfFihUfeKnv4OBguSH6hOLj4/k8KYn+aWl4e3nh6OGBUZ5ZVpY8yu3Z0xj3VqRSE7re3t6cOXPG3GUIJpCZmUlCQgJJSUkkJCTg6upK2bJl2bVrF/v372fixIl07Nix0NdIkoRCoWDbtm1ERUXx0ksvsXz5cmxtbenQoYNR69Nqtdy9e7fYV+izs7MfORJ98cUXH/i4k5NTqQnR4srLyyMiIoItW7bw9ttv49OkCZpp0+RdwZ71yB2dDnJz5Qt0rq5GqbcopSp0v//+e3OXIZSwnJwcQkNDOXPmDDVq1GDHjh2EhobSt29fypQpQ1JSEufPn+f1119Hec8voyF0IyMjeffdd+nduzdLly7lu+++o3bt2pQv4qq2TqcjNTW12FfoMzMzcXV1fegV+hdeeOGBEaqTk1OhWoV/SZLEjz/+yKJFi6hRowabNm3C19dXDttff4X//U8Oyqf9/ul0kJoK770HrVoZtfZHsf7QTUuDY8eo/u23vBUdLa9YUauhcmUIDIR69aBaNXNXKRiJnZ0dmzZtKvR+/fr1adq0KU2bNiUxMZHs7Gz0en2hIFMqldy6dYsqVaoUTCts0aIFGzdu5MKFCwWhm5eXx549e/j2228LgjQ9PR0XF5eHXqGvXr36AyNRFxcXEaJGcPXqVRYuXMjNmzeZPHkyjRo1+veTSiUsXgySJG/H+DSbmefkyH3c4GCYPNm4xRfBekP3779h7VrYvRskCc+cHKokJMhHL0sSHD8Ou3bJG1jUrAlDhkCHDmY5/VMwLsOo9fz58+Tn5xc6xtvBwYFr1649dLOcnJwcHBwcyM3NBeQzutRqNSn3rFJSqVTUq1ePqlWrFoSsm5ubCFETysrKYv369ezZs6fgVclDpwja2MDy5RARAXPmQEaGPOp9XI9eq5UHa46OsGgRdO78ZEf9PCPrC12dDj7/XP5mabXyRsZqNQpnZ7Lu3EFnb18wbQSQA/j8efmk0MBA+euKM0FaMBtJkkhPTyclJQUfH58H5nUaQve3335DkiSaNm1a8DkvLy/OnDmDVqstPKke+ZQRe3v7glkukiRhY2NTKKBVKhUBAQEl+OyER5EkiYMHD7JkyRICAwP54osv8PT0LPqLlEp491149VVYuhT27ZNbD5IknwhhCGCtVh7ZqlTyK+F33oGRI+GeP9imYl2hm5kpn1f/66/y3pf3/PVTAGqNBm1+Pqp7f0kVCjmYJQlOnID27WH9emjY0PT1P6cMqwWLM73J8GZnZ4eHhweLFi16YO614SJSfHw8mZmZeHp6kpeXh42NDe7u7qSmppKdnf3AkmEvLy98fX05dOgQw4YNQ6vV8vvvvzNixAiTfS+Eh4uNjWXBggVkZGQwe/Zs6j7pqrBKleR2w5Qp8P338ivdEyfg7l05Azw84D//kQderVvL+WEm1hO6ublys/v4cfkUz4e8HNBoNORrtQ9f8aJQyC89srLkHs6WLVC/fsnXXQoZln4WZ3qT4Sq+Yennvcs7i7308z4KhYK0tDRu375N48aNAQqmfdWrV4/IyEg8PDwKbq/T6UhKSsLb25tXX32Vffv2MWXKFC5fvsyLL75I+fLlC0bPgmmlp6ezatUqDhw4QGhoKN26dXu2Vo6HB3TvLr9ZKOsJ3fnz4dgxcHd/ZP9Fo1ajLWqTYgAHB3nEPHCg/BfRzc34tVoZQ4gWNQq9/3MqleqB6U1ubm6ULVu20NJPQ5AWFaJPIyEhgR9++IF169YVPIdvvvmGadOmcfjwYU6ePEmbNm0IDw/n2LFjrFy5ko0bN1KjRg3mzZvHhg0baN68OV26dAHEwgtT0+v1xMTEsGLFClq1asWOHTtwe05+FxVFnUcVGBgoHT9+3ITlPMLJk/KkZWfnIpvkCYmJKBQKvB7XBwL5ZUfHjnIfqBTKyckp9jzRlJSUgqWfD7tCf/8EfJMv/XyIuLg4Bg0axHfffQfIoXvq1Cl+//13/P39KVOmDL6+vri7u5u1TuFBZ86cYd68eajVasaNG0fNmjXNXZLRKRSKE5IkBT70c1YRul26yLMSXFyKvFlKSgrZOTn4liv3+PvU6+UrmF9+CTVqGKfOEmRY+lncNfRFLf182AjVopd+CqVCcnIyK1as4Oeff2bkyJEPzKUuTYoKXctvL5w/Lx+N/JiVItuTk9l95w4XsrIIkiSm+foWfb9KpXxxLSICZs0yYsHFY1j6Wdw19EUt/fTz8yvVSz8F66bT6dixYwfr1q3jjTfeYOfOnTiZ8UKWuVl+6O7dK49KHxMgnmo173l4cOBJjnZxdpbn8s6Y8cybXBiWfhb3Cn1RSz99fX3F0k8j0+v1ZGdnY2Njg1qtFt9LE7l27Rrjxo3D3d2dNWvWiF0AsYbQPXIEirH/ZmsXF7Q6HT/Hxxf/vtVqebR7+TJUrVroUzqd7oH180W9vM/IyMDNzU0s/bRQSqUSSZIYN24cfn5+vP/+++L7bwQnTpzAx8fnkcuoPT09GT58OM2bNxd/6P5h2T1dSZJXk9nbF2skKgEf//EHeHgw474fAgnQabXodDq0Oh06rRatTocyI4MDHTpw3Nu7UJBmZGQ8cunnw3qjYumndcjIyGDUqFFUrlyZKVOmiP9nTykrK4tBgwZx7NgxAgICqFu3Lp988skDwfq8TsWz3p5ufr48P7eY/R/D2fYpKSlc1mpBktDp9Wi1WvQ6Hcp7zlIy/Oug1+Pn6op727aFRqhi6WfpZDjifezYsXz00UdMnz4dtdqyfw0s0alTp0hPTyc2Npbr168TEhLC5s2b6dKlS6F+7fMYuI9j2anyFP/DXP+54JaTk0NmZiZ6nQ57e3u8vLwo7+uLX6VKVKpUiQoVKlDOxwdXFxcaN25M27ZtCQwMxN/fHw8PDxG4pZiDgwNLliwhIyODiRMniiPen4K7uztZWVkkJiZSoUIF+vfvz08//cTvv/9u7tIsnmUni1otr59+giOR7f45OfSFGjWoWasWFStWxMXFBZ1Ox53kZP7+56jmK1eucOv2bdJzcriallawCYrwfDAc8a5UKvnggw/EEe9PKD8/nzp16nD27FkAevXqhY2NDSdOnADktoLwcJYduoYdworxC6GTJPL0evSAHsjT65GQt/5zc3XF29ubSn5+BAQEUNXfnzKenmjUanJzcli8fz+vvvoq3bt3Z9KkSaxfv54ffviB+Ph48cNTihmOeHd1dSUsLIysrCxzl2RRJElizZo1nDp1irt37wLyBWaAl156CbVazbFjx7h+/ToAb7zxRsEKQdFWeDTLvpAGsHChvKP7PWvpH2ZNYiJrkpIKfSzE05MQL69Hf5FWK/eMz5xBq9dz+fJlYmNjC73l5ORQvXr1Qm9Vq1bF3t7eGM9OsAB6vZ7Zs2dz6dIlccT7P3bs2MHcuXOpUqUKbm5ulC9fno8//hilUlmwV/Hhw4fZvHkzL774IsOGDUOSJHr06MG6deuemyW9j2LdK9IuXIA33pAXRxj7r2dKCvTpAzNnFnGTFGJjY7l48SKxsbFcuHCBy5cv4+3tTbVq1ahevToBAQFUr16dcuXKib/wVkqSJBYtWsSpU6dYsWLFcxsaOp2OtLQ0xo4dy4ABA2jRogWRkZGcO3eOWbNmFbzyM/ycf/vtt3z22Wc4ODhw5MgRunbtyoIFC5773wPrDl2Abt3gjz+Me36RXg/p6fIy4CfcP1Wn03HlyhUuXrzIhQsXCkbFmZmZBUFseKtWrdoDWwwKlkmSJMLDwzl06BDh4eGP38u1FMnJyeHDDz8kJyeHPn36UK9ePRwcHEhPT6dHjx5oNBomT55MzZo1H/iDlJyczMGDBylXrhyvvPKKeZ6AhbH+0D11Cnr0eOyGN0/k7l0ICpL34DSS1NTUQiPiixcvEhcXR5kyZQpGxIZQLl++vJghYaHWr1/Pvn37WLlyJd7e3uYup8RFRUUxZcoU2rVrx4svvsiyZcvYvn07NWrU4MMPP0ShUPDCCy9w8uRJ9Ho9S5YsAWD79u0FM36Ewqx3nq5BvXryXrrr1hW5tWOxGY71+Phj49T3D1dXVxo0aECDBg0KPqbX67l69WrBaDgmJobY2FjS0tKoWrXqA6Pi53lNuqV47733sLOzY/DgwYSHh1OhQgVzl1SiypUrR3h4eMEodd++fWRmZgLw0UcfFewo5+fnx8qVK0lOTiYzMxONRiMfEik8EesY6QLk5UH//nD06CM3MS+WzEz5azdvlsPcTNLT0wtGxYa3S5cu4ebmViiIAwICqFChghgVm8Hu3btZt24dK1asoEqVKuYup8QlJyfTuXNnbt26Rf369Rk6dCgNGzYsaI9t3LiRo0ePsmLFCjNXavmsv71gkJUFoaFw+PCTn/4pSfJRy/b28N//wj2jUUuh1+u5fv36A73i5ORk/P39Cy7YGUbFLo/Z6vJJHlepVJKYmEhOTg4VK1Y0yv2WBvv372fp0qUsX7681J+dptPpWL16NcOGDSMyMpLjx48TFBSEr68va9eu5bvvvmP69Ol07tzZ3KVavNITuiAvlNi4EebNk6d83XdW2gMkSb5gptdD48byCRRW9pIoMzPzgVHxxYsXcXZ2fqBX7OfnV3B0zZM+xoQJE4iIiCAmJoaWLVsSERFBSkoKjRs3pn79+g8/kfU58M033zB//nw+/fRTXnzxRXOXUyIetkdC+/btmTp1KgcPHiQhIYGFCxeKi8LFVLpC1+DqVbnHu2OHHKh5eXL4Gl6G5+X9e9Gtdm35CPa2bU161HJJ0uv1xMfHF4yIDaPjxMRE/P39C0bDhtGx6yNmfhh+2YYPH46bmxuHDx8mKiqK6dOnk5eXh0Kh4MiRI4waNYp+/fqZ+Flajh9//JEZM2Ywf/586pmxLWUs+fn5HD9+nLp16z50zvmZM2cYMWIE4eHhpfJkh5JWOkPXIDNTPjvtzz/h9Gl5VKtWg7+/fPBk3brwHPTjDLKysrh06dIDizwcHByoUaMGixYtemAkPG3aNDIyMujevTuzZ88mJiam4NTc6tWrc/DgQdauXVtoGpUhrKOiotiwYQNKpZLZs2eX6l/Qo0ePMnnyZGbNmkWjRo3MXc5T+/XXX1mwYAHly5dnwYIFBefXGdpbkyZN4sKFC4SFhdG3b18zV2udigpdJEl65FuDBg0kwfrp9Xrpxo0b0tGjR6Xc3NxCH9+zZ480fvx4SZIkad++fVLv3r0lSZKk9evXSx07dpRWrFghhYWFSdOmTXvgfs+cOSPVrVtXOnjwoLRw4UIpJCSk0P2XRqdOnZLatGkjHTp0yNylPLGbN29K48aNkzp16iQdOnRI0uv1D9zm1q1b0urVqyWtVmuGCksP4Lj0iFy1/pGu8NTS0tJ46623uHPnDgMHDuS3337D3d2dDz/8kAkTJqBUKrly5QoqlYqoqKhCX5uens7s2bPx9PRk7Nix5ObmUqtWLc6dO1fqe79nz57l/fffZ9y4cbRt29bc5TxWbm4uERERbN26lT59+tCvXz9si3EwgPD0ihrpinlIzzFnZ2eWL1/O0KFD+fXXX9m7dy8HDx6kZcuWeHp68tlnnxETE4OPj0/BRiaGDU+uXLlCcnIyDRs2BODy5cs0atSIS5cume35mEqtWrVYsWIFixYtYt++feYu55EkSeLQoUP06NGD2NhYNm/ezKBBg0Tgmpl1LI4QSoRCoaBGjRrUqFGD4OBg+vbty9y5cxk9ejRlypQpuF1mZmbB1of5+fmoVCquXbuGg4NDwe2uXr2Ku7s7er3eLM/F1KpXr86qVasYNmwYubm5dO/e3dwlFXLlyhUWLlxIfHw8U6ZMKfjjKJifCN3nnKHPdPPmTby8vHB0dGTEiBEsXLiQ7OxsMjMz0el0tGrVCqBgdZKDgwN6vb7gotyZM2dwd3cvFNYAsbGxJCYmUq1aNby8vErVRiiVK1dmzZo1DBs2jJycHN555x1zl0RWVhbr1q1j7969DBgwgF69epX6do+1ET1d4QE6nY6vvvqKo0ePkpCQwKRJk/Dz8+Pnn38mNjaWt99+G4BXXnmFXbt2UaFCBZo2bcrYsWMLjfi0Wi2HDh1i165dxMbGotfrH9gm09/f3+pf7t6+fZuhQ4fSsWNH3nvvvUf/YblzB376Sd5L5PRpeTm6Wg0VK0KjRvJsm3r1nmpaoyRJfP311yxbtozAwEBGjRr1XG3YY2lK95QxwWQSExP566+/CtboL1y4kI0bN+Lt7U3FihVZu3btI88bkySJ5OTkQpsBxcbGcuXKFXx9fQst8AgICKBs2bJWNSq+c+cOQ4cO5ZVXXmH48OGFaz93DlasgK++kt/X6eQTrlUqefFOfr680EetBh8feU55jx5FL/q5R2xsLAsWLCAjI4MJEybw8ssvl8AzFJ6ECF3BqKR7Vi/FxcURFxdHkyZNcHR0fOL7ys/PL9g8/t7lz3l5eQ8dFVvy5vF3795lxIgR1K1blzFjxqDUauGzz2DlSnkBj4tL0bvkSRJkZ8sb6wcEwNKlRW47mpaWxqpVqzh48CBDhgyha9euYo8OCyFCV7A6ycnJhZY+G8618/b2LrQZULVq1Sxq8/j09HTCwsJ4wdeXcX/9heL0aXlL0ic5cdiwdF2hkEfIbdoU+rReryc6Oprw8HBeffVVhg0b9sgVh4J5iNAVSgWtVsvVq1cLbQYUGxtLVlZWoSXP1apVM+vm8VlJSVz4z3+oePcuHlWqoHja0WdOjjzqXbcOWrYE4M8//2T+/Pmo1WrGjx/PCy+8YMTKBWMRoSuUanfv3n1gQ6C4uDi8vLweaFH4+vqW/Evwjz9GiojgWmYmCqVS3rD+aUfi2dmgUJCycyfLt2/n8OHDjBw5ktdff120EiyYCF3huaPT6Qo2j7+3V5yenl4wEr53ZPw0/eiHOnYMevcGFxf0SiU3btxAkiR5T+SnmZUAZFy/zmGlkrMffEBIaKjxahVKjAhdQfhHWlraQ0fF7u7uD4yKn2rz+C5d4OxZ+aIZcmjevHEDrVZLhYoVUT3B/WVmZXHr1i3UKhXlHB2x2bPHrBvvC8UnQlcQimDYXctwwc4wOk5JSSm0ebxhStsjj2g/fx46dnzg5GoJiI+PJy83l4p+fqiUStJ0OmbEx3MkIwM3tZoRXl50+OdiWL5Wy+3bt8nOzsbb2xtnZ2cUKSnyfS9bZoLviPCsrP+MNEEoQUqlEj8/P/z8/HjttdcKPp6RkVFoVPy///2PS5cu4eLi8sCo2M/PD+W+ffLUsPvaCArkc8hu37rFlStX8PPzY+6tW2gUCg4EBHAhJ4ewa9eoZmuLW0YGd5KTcXd3l/vPhvtycYH//U/eJ/pJTkwRLI4Y6QrCE9Dr9dy8ebPQAo8LFy5w584dFl+7RpWsLJTOztjZ2WFna1to72IJSExIIDE9nffy8thRtSp+/wTohMuXsc3KYrCzM94+Ptg8bGFERgbs3g21apno2QpPS4x0BcFIlEolFSpUoEKFCrRu3brg41lZWShr1ybH0ZHsvDzS09LIyc1FpVRi+08A29rZ4eLqyiWdDm16OuUUCvLy8rh9+zZe2dlctrMr+nw6nU5e3SZC16qJ0BUEI3CwswOtFjtPT9z+aQlIyCvucnNyyMnNJT0tjcTcXOJzcrCTJM5fuIBGo6FMmTJUcXLiQlpa0Q+i1cLjbiNYPBG6gmAMD5kOpgBsNBpsNJpCF9/ys7LQX75MxfLlcXB0RKNWk33nDo7Fmdkg5uZaPfF/UBCMQaEANzd5NPoYle3skJRKUu3t0fyzPPhCbi7+j9ttTaMBd3cjFCuYkwhdQTCW2rXlpbuPYa9U0trZmVWJiWTr9fyWlcWh9HQ6Pm7/BIUCxLJfqydCVxCMpWnTYo10ASb6+JCr19P2wgUm37jBJB+foke6Wq0cuv7+RipWMBfR0xUEY+nUCRYvlufqPqb36qJSsaiomQr3S0+HPn2K3hpSsApipCsIxlKxojzaNfYMA0OI9+tn3PsVzEKEriAY05QpckDm5xvvPlNT5T0ditjQXLAeInQFwZgCAiAsTF49VsRqz2JLT4cyZeCjj579vgSLIEJXEIxtyBBo2xbu3pVbA08rI0OeJrZ+fcGuZYL1E6ErCMamUsm7gXXsKLcGijGNrBC9Xg5sOzvYsgVeeqlEyhTMQ4SuIJQEGxv5YMlFi+T3U1Lko3eKYgjbtDRo3Rq++Qbq1CnxUgXTElPGBKGkKBTyBbBmzeQR64YNco/WMJfXcAS7ViuHtF4v33bwYPlfCzlsUzAusbWjIJhKXh78+ae8U9jp0/Ko1sYGqlSRR7R16oC3t7mrFIxAbO0oCJbAxgbq15ff+vQxdzWCmYieriAIggmJ0BUEQTAhEbqCIAgmJEJXEATBhEToCoIgmFCRU8YUCkUicMV05QiCIJQKlSRJ8nrYJ4oMXUEQBMG4RHtBEATBhEToCoIgmJAIXUEQBBMSoSsIgmBCInQFQRBM6P8BISffgv7p8boAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def draw_graph(G, colors, pos):\n",
    "    default_axes = plt.axes(frameon=True)\n",
    "    nx.draw_networkx(G, node_color=colors, node_size=600, alpha=.8, ax=default_axes, pos=pos)\n",
    "    edge_labels = nx.get_edge_attributes(G, 'weight')\n",
    "    nx.draw_networkx_edge_labels(G, pos=pos, edge_labels=edge_labels)\n",
    "    \n",
    "# Draw the graph\n",
    "colors = ['r' for node in graph.nodes]\n",
    "pos = [graph.nodes[node]['pos']  for node in graph.nodes]\n",
    "draw_graph(graph, colors, pos)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "167b03d6",
   "metadata": {},
   "source": [
    "## Quadratic Program"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "10d676b7",
   "metadata": {},
   "outputs": [],
   "source": [
    "from docplex.mp.model import Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "6f3f047b",
   "metadata": {},
   "outputs": [],
   "source": [
    "mdl = Model(name=\"TSP\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "1a92edf5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{(0, 0): docplex.mp.Var(type=B,name='x_0_0'), (0, 1): docplex.mp.Var(type=B,name='x_0_1'), (0, 2): docplex.mp.Var(type=B,name='x_0_2'), (0, 3): docplex.mp.Var(type=B,name='x_0_3'), (1, 0): docplex.mp.Var(type=B,name='x_1_0'), (1, 1): docplex.mp.Var(type=B,name='x_1_1'), (1, 2): docplex.mp.Var(type=B,name='x_1_2'), (1, 3): docplex.mp.Var(type=B,name='x_1_3'), (2, 0): docplex.mp.Var(type=B,name='x_2_0'), (2, 1): docplex.mp.Var(type=B,name='x_2_1'), (2, 2): docplex.mp.Var(type=B,name='x_2_2'), (2, 3): docplex.mp.Var(type=B,name='x_2_3'), (3, 0): docplex.mp.Var(type=B,name='x_3_0'), (3, 1): docplex.mp.Var(type=B,name='x_3_1'), (3, 2): docplex.mp.Var(type=B,name='x_3_2'), (3, 3): docplex.mp.Var(type=B,name='x_3_3')}\n"
     ]
    }
   ],
   "source": [
    "x = dict()\n",
    "for i in range(n):\n",
    "    for j in range(n):\n",
    "        x[(i, j)] = mdl.binary_var(name=\"x_{0}_{1}\".format(i, j))\n",
    "print(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1d86883d",
   "metadata": {},
   "source": [
    "## $C\\left(x\\right) = \\sum\\limits_{i=1,j = 1\\\\ i \\ne j}^{n} w_{ij} \\sum\\limits_{p = 1}^{n} x_{i,p} x_{j,p+1}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "5b357ec0",
   "metadata": {},
   "outputs": [],
   "source": [
    "C_x = mdl.sum(\n",
    "            M[i, j] * x[(i, k)] * x[(j, (k + 1) % n)]\n",
    "            for i in range(n)\n",
    "            for j in range(n)\n",
    "            for k in range(n)\n",
    "            if i != j\n",
    "        )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "7162621d",
   "metadata": {},
   "outputs": [],
   "source": [
    "mdl.minimize(C_x)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "65aeb3e0",
   "metadata": {},
   "source": [
    "## $Constraints \\\\ \\sum\\limits_{i=1}^{n} x_{i,p} = 1,  \\ \\forall p \\\\ \\sum\\limits_{p=1}^{n} x_{i,p} = 1, \\ \\forall i$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "94f17a06",
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(n):\n",
    "    mdl.add_constraint(mdl.sum(x[(i, p)] for p in range(n)) == 1)\n",
    "for p in range(n):\n",
    "    mdl.add_constraint(mdl.sum(x[(i, p)] for i in range(n)) == 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "66e052ca",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\\ This file has been generated by DOcplex\n",
      "\\ ENCODING=ISO-8859-1\n",
      "\\Problem name: TSP\n",
      "\n",
      "Minimize\n",
      " obj: [ 96 x_0_0*x_1_1 + 96 x_0_0*x_1_3 + 182 x_0_0*x_2_1 + 182 x_0_0*x_2_3\n",
      "      + 66 x_0_0*x_3_1 + 66 x_0_0*x_3_3 + 96 x_0_1*x_1_0 + 96 x_0_1*x_1_2\n",
      "      + 182 x_0_1*x_2_0 + 182 x_0_1*x_2_2 + 66 x_0_1*x_3_0 + 66 x_0_1*x_3_2\n",
      "      + 96 x_0_2*x_1_1 + 96 x_0_2*x_1_3 + 182 x_0_2*x_2_1 + 182 x_0_2*x_2_3\n",
      "      + 66 x_0_2*x_3_1 + 66 x_0_2*x_3_3 + 96 x_0_3*x_1_0 + 96 x_0_3*x_1_2\n",
      "      + 182 x_0_3*x_2_0 + 182 x_0_3*x_2_2 + 66 x_0_3*x_3_0 + 66 x_0_3*x_3_2\n",
      "      + 126 x_1_0*x_2_1 + 126 x_1_0*x_2_3 + 142 x_1_0*x_3_1 + 142 x_1_0*x_3_3\n",
      "      + 126 x_1_1*x_2_0 + 126 x_1_1*x_2_2 + 142 x_1_1*x_3_0 + 142 x_1_1*x_3_2\n",
      "      + 126 x_1_2*x_2_1 + 126 x_1_2*x_2_3 + 142 x_1_2*x_3_1 + 142 x_1_2*x_3_3\n",
      "      + 126 x_1_3*x_2_0 + 126 x_1_3*x_2_2 + 142 x_1_3*x_3_0 + 142 x_1_3*x_3_2\n",
      "      + 184 x_2_0*x_3_1 + 184 x_2_0*x_3_3 + 184 x_2_1*x_3_0 + 184 x_2_1*x_3_2\n",
      "      + 184 x_2_2*x_3_1 + 184 x_2_2*x_3_3 + 184 x_2_3*x_3_0 + 184 x_2_3*x_3_2\n",
      "      ]/2\n",
      "Subject To\n",
      " c0: x_0_0 + x_0_1 + x_0_2 + x_0_3 = 1\n",
      " c1: x_1_0 + x_1_1 + x_1_2 + x_1_3 = 1\n",
      " c2: x_2_0 + x_2_1 + x_2_2 + x_2_3 = 1\n",
      " c3: x_3_0 + x_3_1 + x_3_2 + x_3_3 = 1\n",
      " c4: x_0_0 + x_1_0 + x_2_0 + x_3_0 = 1\n",
      " c5: x_0_1 + x_1_1 + x_2_1 + x_3_1 = 1\n",
      " c6: x_0_2 + x_1_2 + x_2_2 + x_3_2 = 1\n",
      " c7: x_0_3 + x_1_3 + x_2_3 + x_3_3 = 1\n",
      "\n",
      "Bounds\n",
      " 0 <= x_0_0 <= 1\n",
      " 0 <= x_0_1 <= 1\n",
      " 0 <= x_0_2 <= 1\n",
      " 0 <= x_0_3 <= 1\n",
      " 0 <= x_1_0 <= 1\n",
      " 0 <= x_1_1 <= 1\n",
      " 0 <= x_1_2 <= 1\n",
      " 0 <= x_1_3 <= 1\n",
      " 0 <= x_2_0 <= 1\n",
      " 0 <= x_2_1 <= 1\n",
      " 0 <= x_2_2 <= 1\n",
      " 0 <= x_2_3 <= 1\n",
      " 0 <= x_3_0 <= 1\n",
      " 0 <= x_3_1 <= 1\n",
      " 0 <= x_3_2 <= 1\n",
      " 0 <= x_3_3 <= 1\n",
      "\n",
      "Binaries\n",
      " x_0_0 x_0_1 x_0_2 x_0_3 x_1_0 x_1_1 x_1_2 x_1_3 x_2_0 x_2_1 x_2_2 x_2_3 x_3_0\n",
      " x_3_1 x_3_2 x_3_3\n",
      "End\n",
      "\n"
     ]
    }
   ],
   "source": [
    "\"\"\"\n",
    "DoCPLEX to Quadratic\n",
    "# Convert inequality constraints into equality constraints by adding slack variables\n",
    "# Map integer variables to binary variables\n",
    "# Penalize linear equality constraints with only binary variables\n",
    "\"\"\"\n",
    "mod = QuadraticProgram()\n",
    "mod.from_docplex(mdl)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "f78e62a3",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Convert binary to integer for the paths\n",
    "def route_x(x):\n",
    "    n = int(np.sqrt(len(x)))\n",
    "    route = []\n",
    "    for p__ in range(n):\n",
    "        for i in range(n):\n",
    "            if x[i * n + p__]:\n",
    "                route.append(i)\n",
    "    return route"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "22331bf2",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "QAOA:\n",
      " optimal function value: 236.0\n",
      "optimal value: [0. 1. 0. 0. 0. 0. 1. 0. 0. 0. 0. 1. 1. 0. 0. 0.]\n",
      "status: SUCCESS\n",
      "\n",
      "Route:\n",
      " [1, 2, 3, 0]\n"
     ]
    }
   ],
   "source": [
    "algorithm_globals.random_seed = 10598\n",
    "quantum_instance = QuantumInstance(Aer.get_backend('qasm_simulator'), seed_simulator=algorithm_globals.random_seed, seed_transpiler=algorithm_globals.random_seed)\n",
    "qaoa_mes = QAOA(quantum_instance=quantum_instance, initial_point=[0., 0.])\n",
    "# solving the Problem using QAOA\n",
    "qaoa = MinimumEigenOptimizer(qaoa_mes) \n",
    "qaoa_result = qaoa.solve(mod)\n",
    "print(\"\\nQAOA:\\n\", qaoa_result)\n",
    "qaoa_result = np.asarray([int(y) for y in reversed(list(qaoa_result))])\n",
    "print(\"\\nRoute:\\n\", route_x(qaoa_result))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "detailed-hometown",
   "metadata": {},
   "source": [
    "## If you have any queries on this notebook please reach to me bala.na@hcl.com"
   ]
  }
 ],
 "metadata": {
  "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.9.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}