samples/traveling-salesperson/traveling-salesperson.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Solving the Traveling Salesperson Problem with Azure Quantum QIO\n",
    "Hello and welcome! In this notebook we will walk you through how you (can) solve the traveling salesperson problem (also known as the traveling salesman problem) with the Azure Quantum quantum inspired optimization (QIO) service.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "The traveling salesperson problem is a well-known optimization problem in which the aim is to find a (near) optimal route through a network of nodes. As you will see later on, it is not a straightforward problem to solve as the complexity (difficulty) grows exponentially with the number of nodes. Additionally, due to the rugged (non-convex) optimization space, it is difficult or even impossible to find an optimal route (global minimum/optimal solution) through a large network. Common solvers for these rugged problems are based on searches, which you will implement in this tutorial with Azure QIO!\n",
    "\n",
    "\n",
    "Imagine you have to calculate an optimal route for the salesperson over $N$ nodes (addresses). Your manager wants you to minimize the traveling cost (time/distance/money/etc.) such that the salesperson is more profitable. There are a number of constraints which have to be satisfied by the salesperson:\n",
    "\n",
    "1. The salesperson is required to visit each node. \n",
    "2. The salesperson may visit each node only once!\n",
    "3. The salesperson starts and finishes in the starting node (headquarters). \n",
    "\n",
    "For more information regarding the traveling salesperson problem, check out these links:\n",
    "\n",
    "- [Introduction to binary optimization - Azure Quantum docs](https://docs.microsoft.com/en-us/azure/quantum/optimization-binary-optimization) \n",
    "- [Traveling salesperson problem - Wikipedia](https://en.wikipedia.org/wiki/Travelling_salesperson_problem)\n",
    "\n",
    "> Please note this sample is intended to demonstrate how to formulate the cost function for a well-understood problem mathematically and then map it to a binary QUBO/PUBO format. The Traveling Salesperson Problem is not a good example of a problem that scales well in this format, as detailed in [this paper](https://arxiv.org/abs/1702.06248).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup\n",
    "\n",
    "First, you need to import some dependencies and connect to your Azure Quantum Workspace.\n",
    "\n",
    "> If you do not have an Azure Quantum workspace, or an Azure subscription, please take a look at the [Get started with Azure Quantum](https://docs.microsoft.com/learn/modules/get-started-azure-quantum/) Learn module, or alternatively the [Azure Quantum docs site](https://docs.microsoft.com/azure/quantum/how-to-create-quantum-workspaces-with-the-azure-portal)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import dependencies\n",
    "import numpy as np\n",
    "import os\n",
    "import time\n",
    "import math\n",
    "import requests\n",
    "import json\n",
    "import datetime  \n",
    "\n",
    "from azure.quantum.optimization import Problem, ProblemType, Term, HardwarePlatform, Solver\n",
    "from azure.quantum.optimization import SimulatedAnnealing, ParallelTempering, HardwarePlatform, Tabu, QuantumMonteCarlo\n",
    "from typing import List"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This allows you to connect to the Workspace you've previously deployed in Azure.\n",
    "# Be sure to fill in the settings below which can be retrieved by running 'az quantum workspace show' in the terminal.\n",
    "from azure.quantum import Workspace\n",
    "\n",
    "# Copy the settings for your workspace below\n",
    "workspace = Workspace (\n",
    "    subscription_id = \"\",\n",
    "    resource_group = \"\",\n",
    "    name = \"\",\n",
    "    location = \"\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Defining a cost function: Minimizing the travel cost\n",
    "The generated route should minimize the travel cost for the salesperson. It is time to define this mathematically. You will first need to define the travel costs between the nodes and have a suitable mapping for the location of the salesperson. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Defining the travel cost matrix\n",
    "\n",
    "Consider a single trip for the salesperson, from one node to another node. They start at node $i$ and travel to node $j$, which requires $c(i,j)$ in travel costs. Here, $i$ denotes the origin node and $j$ denotes the destination node. To keep matters simple for now, both $i$ and $j$ can be any node in the set of $N$ nodes. \n",
    "\n",
    "$$ \\text{The origin node } (\\text{node } i) \\text{ with } i \\in \\{0,N-1\\}.$$\n",
    "$$ \\text{The destination node } (\\text{node } j) \\text{ with } j \\in \\{0,N-1\\}.$$\n",
    "$$ \\text{Time traveling from } \\text{node } i \\text{ to } \\text{node } j \\text{ is } c_{i,j}.$$\n",
    "\n",
    "As you may have already noticed, the travel cost between two nodes can be written out for every $i$ and $j$, which results in a travel cost matrix $C$ (linear algebra), shown below:\n",
    "\n",
    "$$C = \\begin{bmatrix} c_{0,0} & c_{0,1} & \\dots & c_{0,N-1} \\\\ c_{1,0} & c_{1,1} & \\dots & c_{1,N-1} \\\\ \\vdots & \\ddots & \\ddots & \\vdots \\\\ c_{N-1,0} & c_{N-1,1} & \\dots & c_{N-1,N-1} \\end{bmatrix}. $$\n",
    "\n",
    "Here, we define the rows to be origin nodes and the columns to be destination nodes. \n",
    "For example, traveling from node 0 to node 1 is simply described by:\n",
    "\n",
    "$$ C(0,1) = c_{0,1}. $$\n",
    "\n",
    "The unit for the travel cost is arbitrary: it can be time, distance, money, or a combination of these and/or other factors."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Defining the location vectors\n",
    "\n",
    "Now that the travel cost between two nodes has been formulated, a representation of the origin and destination nodes of the salesperson must be defined to specify which element of the matrix gives the associated travel cost. Remember that this is still for a single trip! \n",
    "\n",
    "This is the same as saying \"we need a way to select a row-column pair in the cost matrix\". This can be done by multiplying the cost matrix with a vector from the left, and a vector from the right! The left vector will specify the origin node and the right vector the destination node. For brevity they'll be named the origin vector (left) and the destination vector (right).\n",
    "\n",
    "Consider the example where you tell the salesperson to travel from node 1 to node 2:\n",
    "\n",
    "$$ \\text{Travel cost node 0 to node 1 }=  \\begin{bmatrix} 1 & 0 & \\dots & 0 \\end{bmatrix} \\begin{bmatrix} c_{0,0} & c_{0,1} & \\dots & c_{0,N-1} \\\\ c_{1,0} & c_{1,1} & \\dots & c_{1,N-1} \\\\ \\vdots & \\ddots & \\ddots & \\vdots \\\\ c_{N-1,0} & c_{N-1,1} & \\dots & c_{N-1,N-1} \\end{bmatrix} \\begin{bmatrix} 0 \\\\ 1 \\\\ \\vdots \\\\ 0 \\end{bmatrix} = c_{0,1}.$$\n",
    "\n",
    "> Please note that the salesperson can only visit one node at a time (no magic allowed!) and that there is only one salesperson, thus the sum elements in the origin and destination vector must equal 1!\n",
    "\n",
    "Fantastic, you now know how to the express the travel cost for a single trip. However, it is necessary to express these ideas in mathematical format that the solver understands. In the previous example the trip was hard-coded from node 1 to node 2, let's generalize this from any origin node to any destination node:\n",
    "\n",
    "$$ x_k \\in \\{0,1\\} \\text{ for } k \\in \\{0,2N-1\\}, $$\n",
    "\n",
    "$$ \\text{Travel cost for single trip }=  \\begin{bmatrix} x_0 & x_1 & \\dots & x_{N-1} \\end{bmatrix} \\begin{bmatrix} c_{0,0} & c_{0,1} & \\dots & c_{0,N-1} \\\\ c_{1,0} & c_{1,1} & \\dots & c_{1,N-1} \\\\ \\vdots & \\ddots & \\ddots & \\vdots \\\\ c_{N-1,0} & c_{N-1,1} & \\dots & c_{N-1,N-1} \\end{bmatrix} \\begin{bmatrix} x_{N} \\\\ x_{N+1}\\\\ \\vdots \\\\ x_{2N-1} \\end{bmatrix}. $$\n",
    "\n",
    "The solver will determine which $x_k$ are given a value of 1 or 0 (this is the binary variable you are optimizing). If the value is 1, it means the salesperson is travelling between the corresponding origin and destination nodes. Correspondingly, if the value is 0 it means the salesperson is not originating from or traveling to that location.\n",
    "\n",
    "> You may be wondering why the destination vector is indexed from $N+1$ to $2N$ (seperate variables). The reason for this is that otherwise, the solver would consider the left vector and the right vector equal. As as consequence, there would be an origin vector on both sides of the cost matrix, meaning the salesperson would remain at the origin node, which isn't allowed (it would count as visiting the same node more than once, which isn't allowed according to the constraints set out at the start). \n",
    "\n",
    "To summarize: for the salesperson to be at one node at a time, the sum of the vector elements for both the origin and destination vectors has to be 1.\n",
    "\n",
    "$$ \\text{Sum of the origin vector elements: }\\sum_{k = 0}^{N-1} x_k = 1, $$ \n",
    "$$ \\text{Sum of the destination vector elements: }\\sum_{k = N}^{2N-1} x_k = 1.$$ \n",
    "\n",
    "Later in this notebook, constraints will be designed to disable magic/superposition for the salesperson.\n",
    "\n",
    "Now that the basic mathematical formulations are covered, the scope can be expanded. Let's take a look at how two trips can be modeled!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. Defining the travel costs for a route\n",
    "\n",
    "To derive the cost function for a route through a network, you'll need a way to describe the 'total travel cost'. As you might expect, the total cost of a route through the network is the sum of the travel costs between the nodes (sum of the trips). Say you have a route ($R$) from node 1 to node 3 to node 2. The total cost of the route would then be: \n",
    "\n",
    "$$ \\text{Cost of route: } R_{1-3-2} = c_{1,3} + c_{3,2}. $$\n",
    "\n",
    "Note that for the second trip, the origin node is the same as the destination node of the previous trip. Knowing that the last destination equals the new origin is a useful property when reducing the number of variables the solver has to optimize for. Therefore these vectors can also be called the 'location' vectors. \n",
    "\n",
    "Recall that the costs can be expressed with linear algebra. Then the total cost of two trips is:\n",
    "\n",
    "$$ \\text{Cost of route } = \\begin{bmatrix} x_0 & x_1 & \\dots & x_{N-1} \\end{bmatrix} \\begin{bmatrix} c_{0,0} & c_{0,1} & \\dots & c_{0,N-1} \\\\ c_{1,0} & c_{1,1} & \\dots & c_{1,N-1} \\\\ \\vdots & \\ddots & \\ddots & \\vdots \\\\ c_{N-1,0} & c_{N-1,1} & \\dots & c_{N-1,N-1} \\end{bmatrix} \\begin{bmatrix} x_{N} \\\\ x_{N+1}\\\\ \\vdots \\\\ x_{2N-1} \\end{bmatrix} + \\begin{bmatrix} x_{N} & x_{N+1} & \\dots & x_{2N-1} \\end{bmatrix} \\begin{bmatrix} c_{0,0} & c_{0,1} & \\dots & c_{0,N-1} \\\\ c_{1,0} & c_{1,1} & \\dots & c_{1,N-1} \\\\ \\vdots & \\ddots & \\ddots & \\vdots \\\\ c_{N-1,0} & c_{N-1,1} & \\dots & c_{N-1,N-1} \\end{bmatrix} \\begin{bmatrix} x_{2N} \\\\ x_{2N+1}\\\\ \\vdots \\\\ x_{3N-1} \\end{bmatrix}.$$\n",
    "\n",
    "Generalizing the small example to a route in which the salesperson visits all $N$ nodes and returns back to the starting location gives\n",
    "\n",
    "$$\\text{Travel cost of route } = \\sum_{k=0}^{N-1} \\left(  \\begin{bmatrix} x_{Nk} & x_{Nk+1} & \\dots & x_{Nk+N-1} \\end{bmatrix} \\begin{bmatrix} c_{0,0} & c_{0,1} & \\dots & c_{0,N-1} \\\\ c_{1,0} & c_{1,1} & \\dots & c_{1,N-1} \\\\ \\vdots & \\ddots & \\ddots & \\vdots \\\\ c_{N-1,0} & c_{N-1,1} & \\dots & c_{N-1,N-1} \\end{bmatrix} \\begin{bmatrix} x_{N(k+1)} \\\\ x_{N(k+1)+1}\\\\ \\vdots \\\\ x_{N(k+1)+N-1} \\end{bmatrix} \\right),$$\n",
    "\n",
    "which can equivalently be written as:\n",
    " \n",
    "$$\\text{Travel cost of route} = \\sum_{k=0}^{N-1}\\sum_{i=0}^{N-1}\\sum_{j=0}^{N-1} \\left( x_{Nk+i}\\cdot x_{N(k+1)+j}\\cdot c_{i,j} \\right).$$\n",
    "\n",
    "\n",
    "Fantastic! A cost function to optimize for the salesperson's route has been found! Because you want to minimize (denoted by the 'min') the total travel cost with respect to the variables $x_k$ (given below the 'min'), and make mathematicians happy, you'll want to slightly adjust the model: \n",
    "\n",
    "$$\\text{Travel cost of route} := \\underset{x_0, x_1,\\dots,x_{(N^2+2N)}}{min}\\sum_{k=0}^{N-1}\\sum_{i=0}^{N-1}\\sum_{j=0}^{N-1} \\left( x_{Nk+i}\\cdot x_{N(k+1)+j}\\cdot c_{i,j} \\right).$$\n",
    "\n",
    "Time to write this out in code."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4. Coding the cost function\n",
    "\n",
    "For the solver to find a suitable route, you'll need to specify how it calculates the travel cost for that route. The solver requires you to define a cost term for each possible trip-origin-destination combination given by the variables $k,i,j$, respectively. As described by the cost function, this weighting term is simply the $c(i,j)$ element of the cost matrix. The solver will optimize for the $x$ variables of the location vectors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Define variables\n",
    "\n",
    "# The number of nodes\n",
    "NumNodes = 5\n",
    "\n",
    "# Max cost between nodes \n",
    "maxCost = 10\n",
    "\n",
    "# Node names, to interpret the solution later on\n",
    "NodeName = {0:'A', 1:'B', 2:'C', 3:'D', 4:'E', 5:'F', 6:'G', 7:'H', 8:'I', 9:'J', 10:'K', \n",
    "            11:'L', 12:'M', 13:'N', 14:'O', 15:'P', 16:'Q', 17:'R', 18:'S', 19:'T',\n",
    "            20:'U', 21:'V', 22:'W', 23:'X', 24:'Y', 25:'Z'}\n",
    "\n",
    "# Cost to travel between nodes -- note this matrix is not symmetric (traveling A->B is not same as B->A!)\n",
    "CostMatrix = np.array([[1, 4, 7, 4, 3], [3, 3, 3, 1, 2], [2, 5, 2, 3, 1], [7, 8, 1, 3, 5], [3, 2, 1, 9, 8]])    # If you want to rerun with the same matrix\n",
    "#CostMatrix = np.random.randint(maxCost, size=(NumNodes,NumNodes))                                          # If you want to run with a new cost matrix\n",
    "\n",
    "############################################################################################\n",
    "##### Define the optimization problem for the QIO Solver\n",
    "def OptProblem(CostMatrix) -> Problem:\n",
    "    \n",
    "    # 'terms' will contain the cost function terms for the trips\n",
    "    terms = []\n",
    "\n",
    "    ############################################################################################\n",
    "    ##### Cost of traveling between nodes  \n",
    "    for k in range(0, len(CostMatrix)):  # For each trip\n",
    "        for i in range(0, len(CostMatrix)): # For each origin node\n",
    "            for j in range(0, len(CostMatrix)): # For each destination node\n",
    "                \n",
    "                # Assign a weight to every possible trip from node i to node j for each trip \n",
    "                terms.append(\n",
    "                    Term(\n",
    "                        c = CostMatrix.item((i,j)) ,\n",
    "                        # Plus one to denote dependence on next location\n",
    "                        indices = [i + (len(CostMatrix) * k), j + (len(CostMatrix) * (k + 1))]   \n",
    "                    )\n",
    "                )\n",
    "                ##----- Uncomment one of the below statements if you want to see how the weights are assigned! -------------------------------------------------------------------------------------------------\n",
    "                #print(f'{i + (len(CostMatrix) * k)}, {j + (len(CostMatrix) * (k + 1))}')                                                                 # Combinations of origin and destination nodes \n",
    "                #print(f'For x_{i + (len(CostMatrix) * k)}, to x_{j + (len(CostMatrix) * (k + 1))} in trip number {k} costs: {CostMatrix.item((i,j))}')   # In a format for the solver (as formulated in the cost function)\n",
    "                #print(f'For node {i}, to node {j} in trip number {k} costs: {CostMatrix.item((i,j))}')                                         # In a format that is easier to read for a human\n",
    "\n",
    "    return Problem(name=\"Traveling Salesperson\", problem_type=ProblemType.pubo, terms=terms)\n",
    "\n",
    "OptimizationProblem = OptProblem(CostMatrix)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Defining optimization constraints: Penalizing invalid routes\n",
    "The modeled cost function will allow you to find the cheapest route for the salesperson, however it does not include any information on invalid routes! It is now time to integrate constraints/penalties for routes which the salesperson should not travel."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Constraint 1: The salesperson may not be at more than one node at a time (no magic)\n",
    "\n",
    "The salesperson can only be at one node at a time. In the defined cost function this constraint is not enforced, and thus the solver can return origin and location vectors for which the sum is larger than one. Such vectors represent invalid solutions. \n",
    "\n",
    "To avoid these solutions, the invalid solution has to be penalized. This is done by modifying the cost function. You can imagine the penalization of the cost function as reshaping/redesigning the rugged optimization landscape such that for invalid solutions the solver cannot find (local) minima. Invalid solutions should have such high cost that the solver is very unlikely to select them, but not so high as to block the solver from moving around that area of the cost function to find other, valid solutions. \n",
    "\n",
    "To ensure that the salesperson is only ever at a single location, before or after a trip, we must require that only one element in the origin or destination vector is equal to 1, with the rest being 0. One way of designing the constraint would be to look at the sum of elements of each location vector (the origin and destination vectors):\n",
    "\n",
    "$$ \n",
    "\\begin{align}\n",
    "    \\text{Location vector 0 (HQ)} &: \\text{ } \\hspace{0.5cm}  x_0 + x_1 + \\dots + x_{N-1} = 1, \\\\\n",
    "    \\text{Location vector 1} &: \\text{ } \\hspace{0.5cm} x_{N} + x_{N+1} + \\dots + x_{2N-1} = 1, \\\\\n",
    "    &\\vdots \\\\\n",
    "    \\text{Location vector N (HQ)} &: \\text{ } \\hspace{0.5cm} x_{N^2} + x_{N^2+1} + \\dots + x_{N^2 + N-1} = 1. \n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Enforcing the constraint over all trips would then yield ($N+1$ because the salesperson returns to starting node):\n",
    "\n",
    "$$ \\text{For all locations: } \\hspace{0.5cm}  x_0 + x_1 + \\dots +  x_{N^2 + N-1} = N+1. $$\n",
    " \n",
    "The equation above is a valid way to model the constraint. However, there is a downside to it as well: in this formulation individual locations are penalized but there is no penalty for being in two locations at once! To see this, take the following example:\n",
    "\n",
    "$$\n",
    "\\text{ If the first } N + 1 \\text{ values are all 1:}\\\\\n",
    "x_0=1, x_{1}=1, \\dots, x_{N}=1, \\\\\n",
    "\\text{ } \\\\\n",
    "\\text{ and all other } x \\text{ values are 0:}\\\\\n",
    "x_{N+1}=0, x_{N+2}=0, \\dots, x_{N^2+N-1}=0, \n",
    "$$  \n",
    "\n",
    "the constraint is satisfied but the salesperson is still at all nodes at once. Thus the derived equation is not specific enough to model the location constraint. Let's rethink. \n",
    "\n",
    "Consider three nodes (a length-3 location vector), then if the salesperson is at node 1, they cannot be at node 0 or node 2. If the salesperson is at node 0, they cannot be at node 1 or node 2. Instead of using a sum to express this, an equally valid way would be to use products. The product of elements in a location vector must always be zero, regardless of where the salesperson resides, because only one of the three $x$ values can take value 1 (denoting the salesperson being at that location). Therefore the constraint can also be expressed as: \n",
    "$$ x_0 \\cdot x_1 = 0,$$\n",
    "$$ x_0 \\cdot x_2 = 0,$$ \n",
    "$$ x_1 \\cdot x_2 = 0.$$\n",
    "\n",
    "In this format the constraint is much more specific and stricter for the solver. As a result, the solver will provide solutions that do not violate this constraint. Note that we do not want to count combinations more than once, as this would lead to (assymmetries) inbalances in the cost function. We therefore exclude the reverse combinations:\n",
    "\n",
    "$$ x_{1}\\cdot x_{0}, $$\n",
    "$$ x_{2}\\cdot x_{0}, $$\n",
    "$$ x_{2}\\cdot x_{1}. $$\n",
    "\n",
    "Generalizing the location constraint for a salesperson who passes through all nodes and returns back to the starting node ($N+1$ nodes in total, iterating over $l$):\n",
    "\n",
    "$$ \\sum_{l=0}^{N} \\sum_{i=0}^{N-1} \\sum_{j=0}^{N-1} x_{(i+Nl)} \\cdot x_{(j+Nl)} = 0 \\text{ with } \\{ i,j | i<j \\}  $$\n",
    "\n",
    "Great! Now we have an accurate description of this constraint, lets add it to the code.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Define variables\n",
    "\n",
    "# The number of nodes\n",
    "NumNodes = 5\n",
    "\n",
    "# Max cost between nodes \n",
    "maxCost = 10\n",
    "\n",
    "# Node names, to interpret the solution later on\n",
    "NodeName = {0:'A', 1:'B', 2:'C', 3:'D', 4:'E', 5:'F', 6:'G', 7:'H', 8:'I', 9:'J', 10:'K', \n",
    "            11:'L', 12:'M', 13:'N', 14:'O', 15:'P', 16:'Q', 17:'R', 18:'S', 19:'T',\n",
    "            20:'U', 21:'V', 22:'W', 23:'X', 24:'Y', 25:'Z'}\n",
    "\n",
    "# Cost to travel between nodes -- note this matrix is not symmetric (traveling A->B is not same as B->A!)\n",
    "CostMatrix = np.array([[1, 4, 7, 4, 3], [3, 3, 3, 1, 2], [2, 5, 2, 3, 1], [7, 8, 1, 3, 5], [3, 2, 1, 9, 8]])    # If you want to rerun with the same matrix\n",
    "#CostMatrix = np.random.randint(maxCost, size=(NumNodes,NumNodes))                                          # If you want to run with a new cost matrix\n",
    " \n",
    "############################################################################################\n",
    "##### Define the optimization problem for the Azure Quantum Solver\n",
    "def OptProblem(CostMatrix) -> Problem:\n",
    "    \n",
    "    # 'terms' will contain the cost function terms for the trips\n",
    "    terms = []\n",
    "\n",
    "    ############################################################################################\n",
    "    ##### Cost of traveling between nodes  \n",
    "    for k in range(0, len(CostMatrix)):                          # For each trip\n",
    "        for i in range(0, len(CostMatrix)):                      # For each origin node\n",
    "            for j in range(0, len(CostMatrix)):                  # For each destination node\n",
    "                \n",
    "                #Assign a weight to every possible trip from node i to node j for each trip \n",
    "                terms.append(\n",
    "                    Term(\n",
    "                        c = CostMatrix.item((i,j)),                                     # Element of the cost matrix\n",
    "                        indices = [i + (len(CostMatrix) * k), j + (len(CostMatrix) * (k + 1))]    # +1 to denote dependence on next location\n",
    "                    )\n",
    "                )\n",
    "                ##----- Uncomment one of the below statements if you want to see how the weights are assigned! -------------------------------------------------------------------------------------------------\n",
    "                #print(f'{i + (len(CostMatrix) * k)}, {j + (len(CostMatrix) * (k + 1))}')                                                                   # Combinations between the origin and destination nodes \n",
    "                #print(f'For x_{i + (len(CostMatrix) * k)}, to x_{j + (len(CostMatrix) * (k + 1))} in trip number {k} costs: {CostMatrix.item((i,j))}')     # In a format for the solver (as formulated in the cost function)\n",
    "                #print(f'For node_{i}, to node_{j} in trip number {k} costs: {CostMatrix.item((i,j))}')                                                     # In a format that is easier to read for a human\n",
    "    \n",
    "    ############################################################################################\n",
    "    ##### Constraint: Location constraint - salesperson can only be at 1 node at a time.\n",
    "    for l in range(0,len(CostMatrix)+1):                # The total number of nodes that are visited over the route (+1 because returning to starting node)\n",
    "        for i in range(0,len(CostMatrix)):              # For each origin node\n",
    "            for j in range(0,len(CostMatrix)):          # For each destination node\n",
    "                if i!=j and i<j:                        # i<j because we don't want to penalize twice // i==j is forbidden (above)\n",
    "                    terms.append(\n",
    "                        Term(\n",
    "                            c = int(2 * np.max(CostMatrix)),                                     # Assign a weight penalty dependent on maximum distance from the cost matrix elements\n",
    "                            indices = [i + (len(CostMatrix) * l), j + (len(CostMatrix) * l)]                   \n",
    "                        )\n",
    "                    )\n",
    "                    ##----- Uncomment one of the below statements if you want to see how the weights are assigned! -------------------------------------------------------------------------------------------------\n",
    "                    #print(f'{i + (len(CostMatrix) * k)}, {j + (len(CostMatrix) * (k))}')\n",
    "                    #print(f'Location constraint 1: x_{i + (len(CostMatrix) * l)} - x_{j + (len(CostMatrix) * (l + 1))} (trip {l}) assigned weight: {int(2 * np.max(CostMatrix))}')  # In a format for the solver (as formulated in the cost function)\n",
    "\n",
    "    return Problem(name=\"Traveling Salesperson\", problem_type=ProblemType.pubo, terms=terms)\n",
    "\n",
    "OptimizationProblem = OptProblem(CostMatrix)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Constraint 2: The salesperson must be somewhere, they can't disappear\n",
    "\n",
    "Due to the first constraint, the salesperson is penalized for being in multiple nodes at once. But due do the formulation of this constraint it is possible that the solver puts all all $x_k$ in a location vector equal to zero, meaning that the salesperson could be 'nowhere' after some trip.\n",
    "\n",
    "> By looking at the code and previous sections you can come to the conclusion that the minimal cost is obtained by setting all $x_k$ to 0.\n",
    "\n",
    "To encourage the salesperson to not disappear, it is necessary to reward them for being somewhere. Rewards can be assigned by incorporating negatively weighted terms to the cost function that decrease the cost function value for valid solutions of the optimization problem (remember we are minimizing the cost here, so reducing cost makes it more likely that a solution will be chosen). To demonstrate this point, take the following optimization problem:\n",
    "\n",
    "$$\\text{f(x)} := \\underset{x_0, x_1, x_2}{min} x_0 + x_1 - x_2.$$ \n",
    "$$ \\text{with } x_0, x_1, x_2 \\in \\{0,1\\}$$\n",
    "\n",
    "The minimum value for the example is attained for the solution $x_0$ = 0, $x_1$ = 0, $x_2$ = 1, with the optimal function value equal to -1. Here, the negatively weighted third term encourages $x_2$ to take a value 1 rather than 0, unlike $x_0$ and $x_1$. With this idea in mind, you can stop the salesperson from dissapearing! If the salesperson has to visit $N+1$ nodes over a route then $N+1$ of the $x_k$ variables must be assigned the value 1. Written as an equation: \n",
    "\n",
    "$$ \\sum_{k=0}^{N(N+1)-1} x_k = N+1,$$\n",
    "\n",
    "with $(N(N+1))$ equal to the number of variables used to represent the origin and destination nodes for each trip. You could split the equation for each location vector seperately, but if you keep the constraint linear in $x_k$ the resulting cost function will be the same. To assign a reward to the salesperson for being at a node, the $x_k$ terms are moved to the right side of the equation:\n",
    "\n",
    "$$ 0 = (N+1) -\\left( \\sum_{k=0}^{N(N+1)-1} x_k \\right).$$\n",
    "\n",
    "As you may be aware, there is no guarantee which particular $x_k$ the solver will assign the value 1 to in this equation. However, in the previous constraint the salesperson is already enforced to be in a maximum of one node before/after each trip. Therefore with the cost function weights properly tuned, it can be assumed that the salesperson will not be in two or more nodes at once. \n",
    "\n",
    "In other words, the previous constraint penalizes the salesperson for being at more than one node at a once, while this constraint rewards them for being at as many nodes before/after any trip (also being in multiple nodes at once). The weights of the constraints will effectively determine how well they are satisfied, a balance between the two needs to be found such that both are adhered to. \n",
    "\n",
    "Incorporating this constraint in the code can be done by assigning a negative term to each $x_k$. The $N+1$ value can be ignored, since it results in a linear shift of the optimization landscape and does not have an effect on solutions of the the minimization problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Define variables\n",
    "\n",
    "# The number of nodes\n",
    "NumNodes = 5\n",
    "\n",
    "# Max cost between nodes \n",
    "maxCost = 10\n",
    "\n",
    "# Node names, to interpret the solution later on\n",
    "NodeName = {0:'A', 1:'B', 2:'C', 3:'D', 4:'E', 5:'F', 6:'G', 7:'H', 8:'I', 9:'J', 10:'K', \n",
    "            11:'L', 12:'M', 13:'N', 14:'O', 15:'P', 16:'Q', 17:'R', 18:'S', 19:'T',\n",
    "            20:'U', 21:'V', 22:'W', 23:'X', 24:'Y', 25:'Z'}\n",
    "\n",
    "# Cost to travel between nodes -- note this matrix is not symmetric (traveling A->B is not same as B->A!)\n",
    "CostMatrix = np.array([[1, 4, 7, 4, 3], [3, 3, 3, 1, 2], [2, 5, 2, 3, 1], [7, 8, 1, 3, 5], [3, 2, 1, 9, 8]])    # If you want to rerun with the same matrix\n",
    "#CostMatrix = np.random.randint(maxCost, size=(NumNodes,NumNodes))                                          # If you want to run with a new cost matrix\n",
    " \n",
    "############################################################################################\n",
    "##### Define the optimization problem for the Quantum Inspired Solver\n",
    "def OptProblem(CostMatrix) -> Problem:\n",
    "    \n",
    "    #'terms' will contain the weighting terms for the trips!\n",
    "    terms = []\n",
    "\n",
    "    ############################################################################################\n",
    "    ##### Cost of traveling between nodes  \n",
    "    for k in range(0, len(CostMatrix)):                          # For each trip\n",
    "        for i in range(0, len(CostMatrix)):                      # For each origin node\n",
    "            for j in range(0, len(CostMatrix)):                  # For each destination node\n",
    "                \n",
    "                #Assign a weight to every possible trip from node i to node j for each trip \n",
    "                terms.append(\n",
    "                    Term(\n",
    "                        c = CostMatrix.item((i,j)),                                     # Element of the cost matrix\n",
    "                        indices = [i + (len(CostMatrix) *k ), j + (len(CostMatrix) * (k + 1))]    # +1 to denote dependence on next location\n",
    "                    )\n",
    "                )\n",
    "                ##----- Uncomment one of the below statements if you want to see how the weights are assigned! -------------------------------------------------------------------------------------------------\n",
    "                #print(f'{i + (len(CostMatrix) * k)}, {j + (len(CostMatrix) * (k + 1))}')                                                                   # Combinations between the origin and destination nodes \n",
    "                #print(f'For x_{i + (len(CostMatrix) * k)}, to x_{j + (len(CostMatrix) * (k + 1))} in trip number {k} costs: {CostMatrix.item((i,j))}')     # In a format for the solver (as formulated in the cost function)\n",
    "                #print(f'For node_{i}, to node_{j} in trip number {k} costs: {CostMatrix.item((i,j))}')                                                     # In a format that is easier to read for a human\n",
    "    \n",
    "    ############################################################################################\n",
    "    ##### Constraint: Location constraint - salesperson can only be at 1 node at a time.\n",
    "    for l in range(0,len(CostMatrix)+1):                # The total number of nodes that are visited over the route (+1 because returning to starting node)\n",
    "        for i in range(0,len(CostMatrix)):              # For each origin node\n",
    "            for j in range(0,len(CostMatrix)):          # For each destination node\n",
    "                if i!=j and i<j:                        # i<j because we don't want to penalize twice // i==j is forbidden (above)\n",
    "                    terms.append(\n",
    "                        Term(\n",
    "                            c = int(2 * np.max(CostMatrix)),                                     # Assign a weight penalty dependent on maximum distance from the cost matrix elements\n",
    "                            indices = [i + (len(CostMatrix) * l), j + (len(CostMatrix) * l)]                   \n",
    "                        )\n",
    "                    )\n",
    "                    ##----- Uncomment one of the below statements if you want to see how the weights are assigned! -------------------------------------------------------------------------------------------------\n",
    "                    #print(f'{i + (len(CostMatrix) * k)}, {j + (len(CostMatrix) * (k))}')\n",
    "                    #print(f'Location constraint 1: x_{i + (len(CostMatrix) * l)} - x_{j + (len(CostMatrix) * (l + 1))} (trip {l}) assigned weight: {int(2 * np.max(CostMatrix))}')  # In a format for the solver (as formulated in the cost function)\n",
    "    \n",
    "    ############################################################################################\n",
    "    ##### Constraint: Location constraint - encourage the salesperson to be 'somewhere' otherwise all x_k might be 0 (for example).\n",
    "    for v in range(0, len(CostMatrix) + len(CostMatrix) * (len(CostMatrix))):    # Select variable (v represents a node before/after any trip)\n",
    "        terms.append(\n",
    "            Term(\n",
    "                c = int(-1.65 * np.max(CostMatrix)),          # Assign a weight penalty dependent on maximum distance from the cost matrix elements\n",
    "                indices = [v]   \n",
    "            )\n",
    "        )\n",
    "        ##----- Uncomment one of the below statements if you want to see how the weights are assigned! -------------------------------------------------------------------------------------------------\n",
    "        #print(v)\n",
    "        #print(f'Location constraint 2: x_{v} assigned weight: {int(-1.65 * np.max(CostMatrix))}')                                                      # In a format for the solver (as formulated in the cost function)\n",
    "        #print(f'Location constraint 2: node_{v % NumNodes} after {np.floor(v / NumNodes)} trips assigned weight: {int(-1.65 * np.max(CostMatrix))}')   # In a format that is easier to read for a human\n",
    "\n",
    "    return Problem(name=\"Traveling Salesperson\", problem_type=ProblemType.pubo, terms=terms)\n",
    "\n",
    "OptimizationProblem = OptProblem(CostMatrix)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Constraint 3: Same node constraint - can't travel to the same node more than once (except the starting node)\n",
    "\n",
    "The salesperson may only visit each node once, meaning that routes containing revisits of a node (other than the starting node) must be penalized. As an example, consider node 3 for each trip (these $x$ describe the same node): \n",
    "\n",
    "$$ \\text{Node 3:   } x_3, x_{3+N}, x_{3+2N}, \\dots $$\n",
    "\n",
    "If the salesperson has been in node 3 after any trip, they may not pass through the same node again. This means that if $x_3$ is 1 for example, $x_{3+N}$ and $x_{3+2N}$ have to be 0. As done similarly with the location constraint, the product of these variables is one way of designing a constraint. As you know that the product between variables representing the same node has to be zero, the following can be derived:\n",
    "\n",
    "$$ x_3 \\cdot x_{3+N} \\cdot x_{3+2N} = 0. $$\n",
    "\n",
    "Even though the equation seems to represent the constraint correctly, it is not stringent enough. With such an equation, if $x_3$ and $x_{3+N}$ are both 1 and $x_{3+2N}$ is 0, the constraint is satisfied with the salesperson having been in the same node twice. Therefore, similarly to the location constraint, more specificity is required. Luckily, the constraint can be split into smaller products:\n",
    "\n",
    "$$ x_3 \\cdot x_{3+N} =0$$\n",
    "$$ x_3 \\cdot x_{3+2N} = 0$$\n",
    "$$ x_{3+N} \\cdot x_{3+2N} = 0$$\n",
    "\n",
    "Continuing this for all N trips and N nodes yields:\n",
    "\n",
    "$$ \\large{\\sum}_{p=0}^{N^2+N-1}\\hspace{0.25cm} \\large{\\sum}_{f=p+N,\\hspace{0.2cm} \\text{stepsize: } N}^{N^2-1} \\hspace{0.35cm} (x_p \\cdot x_f) = 0,$$\n",
    "\n",
    "in which the first summation assigns a reference node $x_p$, and the second summation the same node but after some number of trips $x_f$ (multiple of N: stepsize). \n",
    "\n",
    "This constraint penalizes routes in which nodes are visited more than once. It does not include the last trip back to the starting node (headquarters). The last trip does not need to penalized since the salesperson has already visited each node. Therefore, all remaning travels would result in a violation of the constraint. Incorporating it would only make the cost function larger without adding any value to the optimization problem. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Define variables\n",
    "\n",
    "# The number of nodes\n",
    "NumNodes = 5\n",
    "\n",
    "# Max cost between nodes \n",
    "maxCost = 10\n",
    "\n",
    "# Node names, to interpret the solution later on\n",
    "NodeName = {0:'A', 1:'B', 2:'C', 3:'D', 4:'E', 5:'F', 6:'G', 7:'H', 8:'I', 9:'J', 10:'K', \n",
    "            11:'L', 12:'M', 13:'N', 14:'O', 15:'P', 16:'Q', 17:'R', 18:'S', 19:'T',\n",
    "            20:'U', 21:'V', 22:'W', 23:'X', 24:'Y', 25:'Z'}\n",
    "\n",
    "# Cost to travel between nodes -- note this matrix is not symmetric (traveling A->B is not same as B->A!)\n",
    "CostMatrix = np.array([[1, 4, 7, 4, 3], [3, 3, 3, 1, 2], [2, 5, 2, 3, 1], [7, 8, 1, 3, 5], [3, 2, 1, 9, 8]])    # If you want to rerun with the same matrix\n",
    "#CostMatrix = np.random.randint(maxCost, size=(NumNodes,NumNodes))                                          # If you want to run with a new cost matrix\n",
    " \n",
    "############################################################################################\n",
    "##### Define the optimization problem for the Quantum Inspired Solver\n",
    "def OptProblem(CostMatrix) -> Problem:\n",
    "    \n",
    "    #'terms' will contain the weighting terms for the trips!\n",
    "    terms = []\n",
    "\n",
    "    ############################################################################################\n",
    "    ##### Cost of traveling between nodes  \n",
    "    for k in range(0, len(CostMatrix)):                          # For each trip\n",
    "        for i in range(0, len(CostMatrix)):                      # For each origin node\n",
    "            for j in range(0, len(CostMatrix)):                  # For each destination node\n",
    "                \n",
    "                #Assign a weight to every possible trip from node i to node j for each trip \n",
    "                terms.append(\n",
    "                    Term(\n",
    "                        c = CostMatrix.item((i,j)),                                     # Element of the cost matrix\n",
    "                        indices = [i + (len(CostMatrix) *k ), j + (len(CostMatrix) * (k + 1))]    # +1 to denote dependence on next location\n",
    "                    )\n",
    "                )\n",
    "                ##----- Uncomment one of the below statements if you want to see how the weights are assigned! -------------------------------------------------------------------------------------------------\n",
    "                #print(f'{i + (len(CostMatrix) * k)}, {j + (len(CostMatrix) * (k + 1))}')                                                                   # Combinations between the origin and destination nodes \n",
    "                #print(f'For x_{i + (len(CostMatrix) * k)}, to x_{j + (len(CostMatrix) * (k + 1))} in trip number {k} costs: {CostMatrix.item((i,j))}')     # In a format for the solver (as formulated in the cost function)\n",
    "                #print(f'For node_{i}, to node_{j} in trip number {k} costs: {CostMatrix.item((i,j))}')                                                     # In a format that is easier to read for a human\n",
    "    \n",
    "    ############################################################################################\n",
    "    ##### Constraint: Location constraint - salesperson can only be at 1 node at a time.\n",
    "    for l in range(0,len(CostMatrix)+1):                # The total number of nodes that are visited over the route (+1 because returning to starting node)\n",
    "        for i in range(0,len(CostMatrix)):              # For each origin node\n",
    "            for j in range(0,len(CostMatrix)):          # For each destination node\n",
    "                if i!=j and i<j:                        # i<j because we don't want to penalize twice // i==j is forbidden (above)\n",
    "                    terms.append(\n",
    "                        Term(\n",
    "                            c = int(2 * np.max(CostMatrix)),                                     # Assign a weight penalty dependent on maximum distance from the cost matrix elements\n",
    "                            indices = [i + (len(CostMatrix) * l), j + (len(CostMatrix) * l)]                   \n",
    "                        )\n",
    "                    )\n",
    "                    ##----- Uncomment one of the below statements if you want to see how the weights are assigned! -------------------------------------------------------------------------------------------------\n",
    "                    #print(f'{i + (len(CostMatrix) * k)}, {j + (len(CostMatrix) * (k))}')\n",
    "                    #print(f'Location constraint 1: x_{i + (len(CostMatrix) * l)} - x_{j + (len(CostMatrix) * (l + 1))} (trip {l}) assigned weight: {int(2 * np.max(CostMatrix))}')  # In a format for the solver (as formulated in the cost function)\n",
    "    \n",
    "    ############################################################################################\n",
    "    ##### Constraint: Location constraint - encourage the salesperson to be 'somewhere' otherwise all x_k might be 0 (for example).\n",
    "    for v in range(0, len(CostMatrix) + len(CostMatrix) * (len(CostMatrix))):    # Select variable (v represents a node before/after any trip)\n",
    "        terms.append(\n",
    "            Term(\n",
    "                c = int(-1.65 * np.max(CostMatrix)),          # Assign a weight penalty dependent on maximum distance from the cost matrix elements\n",
    "                indices = [v]   \n",
    "            )\n",
    "        )\n",
    "        ##----- Uncomment one of the below statements if you want to see how the weights are assigned! -------------------------------------------------------------------------------------------------\n",
    "        #print(v)\n",
    "        #print(f'Location constraint 2: x_{v} assigned weight: {int(-1.65 * np.max(CostMatrix))}')                                                      # In a format for the solver (as formulated in the cost function)\n",
    "        #print(f'Location constraint 2: node_{v % NumNodes} after {np.floor(v / NumNodes)} trips assigned weight: {int(-1.65 * np.max(CostMatrix))}')   # In a format that is easier to read for a human\n",
    "\n",
    "    ############################################################################################                        \n",
    "    ##### Penalty for traveling to the same node again --- (in the last step we can travel without penalties (this is to make it easier to specify an end node =) ))\n",
    "    for p in range(0, len(CostMatrix) + len(CostMatrix) * (len(CostMatrix))):                                  # This selects a present node x: 'p' for present    \n",
    "        for f in range(p + len(CostMatrix), len(CostMatrix) * (len(CostMatrix)), len(CostMatrix)):              # This selects the same node x but after upcoming trips: 'f' for future\n",
    "            terms.append(\n",
    "                Term(\n",
    "                    c = int(2 * np.max(CostMatrix)),                               # assign a weight penalty dependent on maximum distance from the cost matrix elements\n",
    "                    indices = [p,f]   \n",
    "                )\n",
    "            )     \n",
    "            ##----- Uncomment one of the below statements if you want to see how the weights are assigned! -------------------------------------------------------------------------------------------------\n",
    "            #print(f'x_{p}, x_{f}')                                                                                                                                             # Just variable numbers \n",
    "            #print(f'Visit once constraint: x_{p} - x_{f}  assigned weight: {int(2 * np.max(CostMatrix))}')                                                                     # In a format for the solver (as formulated in the cost function)\n",
    "            #print(f' Visit once constraint: node_{p % NumNodes} - node_{(p + f) % NumNodes} after {(f - p) / NumNodes} trips assigned weight: {int(2 * np.max(CostMatrix))}')  # In a format that is easier to read for a human\n",
    "\n",
    "    return Problem(name=\"Traveling Salesperson\", problem_type=ProblemType.pubo, terms=terms)\n",
    "\n",
    "OptimizationProblem = OptProblem(CostMatrix)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Constraint 4 & 5: Beginning and ending at a specific node \n",
    "\n",
    "The salesperson starts and finishes at the headquarters: the starting node from which they depart on their journey. Similarly to constraint 2, the salesperson can be rewarded for starting/finishing at a specific node. \n",
    "\n",
    "For example, if you want the salesperson to start/end at a particular node, you can assign negative weights to the respective $x_k$ term in the first and last location vector. For node 0 that would result in negatively weighting $x_0$ and $x_{N^2}$, for example. This constraint is the most flexible one, and you can also expand it to encourage the salesperson to visit a specific node, or nodes, at pre-determined trip number $k$. Alternatively, you may see this constraint as a way to integrate prior-knowledge on the set of nodes into the optimization problem. Let's say you want the salesperson to visit node 1 after the second trip ($k=2$, it is the third node they visit), then by negatively weighting $x_{2N+1}$ the cost function will (likely) obtain a lower/smaller optimal (minimal) value if it gives this variable the value 1. \n",
    "\n",
    "In this code, the salesperson is hard-coded to start and finish in node 0. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Define variables\n",
    "\n",
    "# The number of nodes\n",
    "NumNodes = 5\n",
    "\n",
    "# Max cost between nodes \n",
    "maxCost = 10\n",
    "\n",
    "# Node names, to interpret the solution later on\n",
    "NodeName = {0:'A', 1:'B', 2:'C', 3:'D', 4:'E', 5:'F', 6:'G', 7:'H', 8:'I', 9:'J', 10:'K', \n",
    "            11:'L', 12:'M', 13:'N', 14:'O', 15:'P', 16:'Q', 17:'R', 18:'S', 19:'T',\n",
    "            20:'U', 21:'V', 22:'W', 23:'X', 24:'Y', 25:'Z'}\n",
    "\n",
    "# Cost to travel between nodes -- note this matrix is not symmetric (traveling A->B is not same as B->A!)\n",
    "CostMatrix = np.array([[1, 4, 7, 4, 3], [3, 3, 3, 1, 2], [2, 5, 2, 3, 1], [7, 8, 1, 3, 5], [3, 2, 1, 9, 8]])    # If you want to rerun with the same matrix\n",
    "#CostMatrix = np.random.randint(maxCost, size=(NumNodes,NumNodes))                                          # If you want to run with a new cost matrix\n",
    " \n",
    "############################################################################################\n",
    "##### Define the optimization problem for the Quantum Inspired Solver\n",
    "def OptProblem(CostMatrix) -> Problem:\n",
    "    \n",
    "    #'terms' will contain the weighting terms for the trips!\n",
    "    terms = []\n",
    "\n",
    "    ############################################################################################\n",
    "    ##### Cost of traveling between nodes  \n",
    "    for k in range(0, len(CostMatrix)):                          # For each trip\n",
    "        for i in range(0, len(CostMatrix)):                      # For each origin node\n",
    "            for j in range(0, len(CostMatrix)):                  # For each destination node\n",
    "                \n",
    "                #Assign a weight to every possible trip from node i to node j for each trip \n",
    "                terms.append(\n",
    "                    Term(\n",
    "                        c = CostMatrix.item((i,j)),                                     # Element of the cost matrix\n",
    "                        indices = [i + (len(CostMatrix) *k ), j + (len(CostMatrix) * (k + 1))]    # +1 to denote dependence on next location\n",
    "                    )\n",
    "                )\n",
    "                ##----- Uncomment one of the below statements if you want to see how the weights are assigned! -------------------------------------------------------------------------------------------------\n",
    "                #print(f'{i + (len(CostMatrix) * k)}, {j + (len(CostMatrix) * (k + 1))}')                                                                   # Combinations between the origin and destination nodes \n",
    "                #print(f'For x_{i + (len(CostMatrix) * k)}, to x_{j + (len(CostMatrix) * (k + 1))} in trip number {k} costs: {CostMatrix.item((i,j))}')     # In a format for the solver (as formulated in the cost function)\n",
    "                #print(f'For node_{i}, to node_{j} in trip number {k} costs: {CostMatrix.item((i,j))}')                                                     # In a format that is easier to read for a human\n",
    "    \n",
    "    ############################################################################################\n",
    "    ##### Constraint: Location constraint - salesperson can only be at 1 node at a time.\n",
    "    for l in range(0,len(CostMatrix)+1):                # The total number of nodes that are visited over the route (+1 because returning to starting node)\n",
    "        for i in range(0,len(CostMatrix)):              # For each origin node\n",
    "            for j in range(0,len(CostMatrix)):          # For each destination node\n",
    "                if i!=j and i<j:                        # i<j because we don't want to penalize twice // i==j is forbidden (above)\n",
    "                    terms.append(\n",
    "                        Term(\n",
    "                            c = int(2 * np.max(CostMatrix)),                                     # Assign a weight penalty dependent on maximum distance from the cost matrix elements\n",
    "                            indices = [i + (len(CostMatrix) * l), j + (len(CostMatrix) * l)]                   \n",
    "                        )\n",
    "                    )\n",
    "                    ##----- Uncomment one of the below statements if you want to see how the weights are assigned! -------------------------------------------------------------------------------------------------\n",
    "                    #print(f'{i + (len(CostMatrix) * k)}, {j + (len(CostMatrix) * (k))}')\n",
    "                    #print(f'Location constraint 1: x_{i + (len(CostMatrix) * l)} - x_{j + (len(CostMatrix) * (l + 1))} (trip {l}) assigned weight: {int(2 * np.max(CostMatrix))}')  # In a format for the solver (as formulated in the cost function)\n",
    "    \n",
    "    ############################################################################################\n",
    "    ##### Constraint: Location constraint - encourage the salesperson to be 'somewhere' otherwise all x_k might be 0 (for example).\n",
    "    for v in range(0, len(CostMatrix) + len(CostMatrix) * (len(CostMatrix))):    # Select variable (v represents a node before/after any trip)\n",
    "        terms.append(\n",
    "            Term(\n",
    "                c = int(-1.65 * np.max(CostMatrix)),          # Assign a weight penalty dependent on maximum distance from the cost matrix elements\n",
    "                indices = [v]   \n",
    "            )\n",
    "        )\n",
    "        ##----- Uncomment one of the below statements if you want to see how the weights are assigned! -------------------------------------------------------------------------------------------------\n",
    "        #print(v)\n",
    "        #print(f'Location constraint 2: x_{v} assigned weight: {int(-1.65 * np.max(CostMatrix))}')                                                      # In a format for the solver (as formulated in the cost function)\n",
    "        #print(f'Location constraint 2: node_{v % NumNodes} after {np.floor(v / NumNodes)} trips assigned weight: {int(-1.65 * np.max(CostMatrix))}')   # In a format that is easier to read for a human\n",
    "\n",
    "    ############################################################################################                        \n",
    "    ##### Penalty for traveling to the same node again --- (in the last step we can travel without penalties (this is to make it easier to specify an end node =) ))\n",
    "    for p in range(0, len(CostMatrix) + len(CostMatrix) * (len(CostMatrix))):                                  # This selects a present node x: 'p' for present    \n",
    "        for f in range(p + len(CostMatrix), len(CostMatrix) * (len(CostMatrix)), len(CostMatrix)):              # This selects the same node x but after upcoming trips: 'f' for future\n",
    "            terms.append(\n",
    "                Term(\n",
    "                    c = int(2 * np.max(CostMatrix)),                               # assign a weight penalty dependent on maximum distance from the cost matrix elements\n",
    "                    indices = [p,f]   \n",
    "                )\n",
    "            )     \n",
    "            ##----- Uncomment one of the below statements if you want to see how the weights are assigned! -------------------------------------------------------------------------------------------------\n",
    "            #print(f'x_{p}, x_{f}')                                                                                                                                             # Just variable numbers \n",
    "            #print(f'Visit once constraint: x_{p} - x_{f}  assigned weight: {int(2 * np.max(CostMatrix))}')                                                                     # In a format for the solver (as formulated in the cost function)\n",
    "            #print(f' Visit once constraint: node_{p % NumNodes} - node_{(p + f) % NumNodes} after {(f - p) / NumNodes} trips assigned weight: {int(2 * np.max(CostMatrix))}')  # In a format that is easier to read for a human\n",
    "\n",
    "\n",
    "    #############################################################################################                        \n",
    "    ##### Begin at x0\n",
    "    terms.append(\n",
    "        Term(\n",
    "            c = int(-10 * np.max(CostMatrix)),                   # Assign a weight penalty dependent on maximum distance from the cost matrix elements\n",
    "            indices = [0]   \n",
    "        )\n",
    "    )\n",
    "\n",
    "    ############################################################################################                        \n",
    "    ##### End at x0\n",
    "    terms.append(\n",
    "        Term(\n",
    "            c = int(-10 * np.max(CostMatrix)),                   # Assign a weight penalty dependent on maximum distance from the cost matrix elements\n",
    "            indices = [len(CostMatrix) * (len(CostMatrix))]   \n",
    "        )    \n",
    "    )\n",
    "\n",
    "    return Problem(name=\"Traveling Salesperson\", problem_type=ProblemType.pubo, terms=terms)\n",
    "\n",
    "\n",
    "OptimizationProblem = OptProblem(CostMatrix)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Submitting the optimization problem\n",
    "\n",
    "The optimization problem subject to the necessary constraints have been defined! Now it is time hand the problem over to the Azure QIO solvers and analyze the routes that are returned to us. \n",
    "\n",
    "Some essential specifications/explanations:\n",
    "\n",
    "1. The problem type is a [PUBO (Polynomial Unconstrained Binary Optimization)](https://docs.microsoft.com/azure/quantum/optimization-binary-optimization) - the variables ($x_k$) that are optimized for can take a value of 0, or 1.\n",
    "\n",
    "2. We submit the problem to an Azure QIO solver, which one is up to you. Each has benefits/drawbacks, and [selecting the best one requires some experimentation](https://docs.microsoft.com/azure/quantum/optimization-which-solver-should-you-use). To read more about the available solvers, please refer to the [Microsoft QIO solver overview page](https://docs.microsoft.com/azure/quantum/provider-microsoft-qio) on the Azure Quantum docs site.\n",
    "\n",
    "3. These optimizers are heuristics, which means they aren't guaranteed to find the optimal solution (or even a valid one, depending on how well you have encoded your cost function). Due to this, it is important to validate the solution returned. It is also recommended to run the solver several times to see if it returns the same solution. Additionally, having longer optimization times (larger timeout) can return better solutions.  \n",
    "\n",
    "4. The solution returned by the solver is heavily influenced by the constraint weights. Feel free to play around with these (called tuning), and discover if you can make the optimization more efficient (speed vs solution quality trade-off). A suggestion would be to make the weights dependent on the cost matrix norm (Frobenius,1,2,etc.). \n",
    "\n",
    "5. Here the hardware implementation is defaulted to CPU. For further information on available hardware, please refer to the [Microsoft QIO solver overview page](https://docs.microsoft.com/azure/quantum/provider-microsoft-qio) on the Azure Quantum docs site."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "############################################################################################\n",
    "##### Choose the solver and parameters --- uncomment if you wish to use a different one  --- timeout = 120 seconds\n",
    "\n",
    "solver = SimulatedAnnealing(workspace, timeout = 120)   \n",
    "#solver = ParallelTempering(workspace, timeout = 120)\n",
    "#solver = Tabu(workspace, timeout = 120)\n",
    "#solver = QuantumMonteCarlo(workspace, sweeps = 2, trotter_number = 10, restarts = 72, seed = 22, beta_start = 0.1, transverse_field_start = 10, transverse_field_stop = 0.1) # QMC is not available parameter-free yet\n",
    "\n",
    "route = solver.optimize(OptimizationProblem)                                        # Synchronously submit the optimization problem to the service -- wait until done.\n",
    "print(route)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Parse the results\n",
    "\n",
    "Below are defined some utility functions which are needed to read and analyze the result returned by the solver."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "############################################################################################\n",
    "##### Read the results returned by the solver - need to make the solution readable\n",
    "def ReadResults(Config: dict, NodeName, CostMatrix, NumNodes):  \n",
    "\n",
    "    #############################################################################################\n",
    "    ##### Read the return result (dictionary) from the solver and sort it\n",
    "    PathChoice = Config.items()\n",
    "    PathChoice = [(int(k), v) for k, v in Config.items()] \n",
    "    PathChoice.sort(key=lambda tup: tup[0]) \n",
    "\n",
    "    #############################################################################################\n",
    "    ##### Initialize variables to understand the routing    \n",
    "    TimeStep=[]                                                     # This will contain an array of times/trips - each node is represented during/for each time/trip interval\n",
    "    Node = []                                                       # This will contain an array of node names \n",
    "    Location = []                                                   # This will contain the locations the salesperson is for each time/trip\n",
    "    RouteMatrixElements = []                                        # This will contain the indices of the cost matrix representing where the salesperson has traveled (to determine total cost)\n",
    "\n",
    "    #############################################################################################\n",
    "    ##### Go through nodes during each timestep/trip to see where the salesperson has been\n",
    "    for Index in PathChoice:\n",
    "        TimeStep.append(math.floor(Index[0] / len(CostMatrix)))         # Time step/trip = the k-th is floor of the index divided by the number of nodes\n",
    "        Node.append(NodeName[(Index[0] % len(CostMatrix))])             # Append node names for each time step\n",
    "        Location.append(Index[1])                                       # Append locations for each time step\n",
    "        if Index[1] == 1:                                               # Save selected node where the salesperson travels to in that trip (if the variable == 1, the salesperson goes to that node)\n",
    "            RouteMatrixElements.append(Index[0] % len(CostMatrix))      # Save the indices (this returns the row index)\n",
    "    SimulationResult = np.array([TimeStep, Node, Location])             # Save all the route data (also where the salesperson did not go during a turn/trip/timestep)\n",
    " \n",
    "    #############################################################################################\n",
    "    ##### Create the route dictionary \n",
    "    k=0                                                                                                             \n",
    "    PathDict = {}                                                                                                                                              \n",
    "    PathDict['Route'] = {}\n",
    "    Path = np.array([['Timestep,', 'Node']])\n",
    "    for i in range(0, (NumNodes * (NumNodes + 1))):\n",
    "        if SimulationResult[2][i] == '1':                                                                       # If the SimulationResult[2][i] (location) == 1, then that's where the salesperson goes/went\n",
    "            Path = np.concatenate((Path, np.array([[SimulationResult[j][i] for j in range(0, 2)]])), axis=0)    # Add the rows where the salesperson DOES travel to Path matrix\n",
    "            PathDict['Route'].update({k: Path[k + 1][1]})                                                       # Save the route to a dictionary\n",
    "            k += 1                                                                                              # Iterable keeps track for the dictionary, but also allows to check for constraint\n",
    "    AnalyzeResult(Path, NumNodes)                                                                               # Check if Path array satisfies other constraints as well (could integrate previous one above in function)\n",
    "\n",
    "    #############################################################################################\n",
    "    ###### Calculate the total cost of the route the salesperson made (can be in time (minutes) or in distance (km))\n",
    "    TotalRouteCost = 0\n",
    "    for trips in range(0, NumNodes):\n",
    "        TotalRouteCost = TotalRouteCost+float(CostMatrix.item(RouteMatrixElements[trips], RouteMatrixElements[trips + 1]))     # The sum of the matrix elements where the salesperson has been (determined through the indices)\n",
    "    PathDict['RouteCost'] = {'Cost':TotalRouteCost}\n",
    "\n",
    "    ##### Return the simulation result in a human understandable way =)\n",
    "    return PathDict\n",
    "\n",
    "\n",
    "############################################################################################\n",
    "##### Check whether the solution satisfies the optimization constraints \n",
    "def AnalyzeResult(Path, NumNodes):\n",
    "\n",
    "    ############################################################################################                        \n",
    "    ##### Check if the number of travels is equal to the number of nodes + 1 (for returning home)\n",
    "    if (len(Path) - 1) != NumNodes + 1:\n",
    "        raise RuntimeError('This solution is not valid -- Number of nodes visited invalid!')\n",
    "    else:\n",
    "        NumNodesPassed = NumNodes\n",
    "        print(f\"Number of nodes passed = {NumNodesPassed}. This is valid!\")\n",
    "\n",
    "    ############################################################################################                        \n",
    "    ##### Check if the nodes are different (except start/end node)\n",
    "    PastNodes = []\n",
    "    for k in range(1, len(Path) - 1):                                                                           # Start to second last node must all be different - skip header so start at 1, skip last node so - 1\n",
    "        for l in range(0, len(PastNodes)):  \n",
    "            if Path[k][1] == PastNodes[l]:\n",
    "                raise RuntimeError('This solution is not valid -- Traveled to a non-starting node more than once')\n",
    "        PastNodes.append(Path[k][1])\n",
    "    print(f\"Number of different nodes passed = {NumNodes}. This is valid!\") \n",
    "\n",
    "    ############################################################################################                        \n",
    "    ##### Check if the end node is same as the start node\n",
    "    if Path[1][1] != Path[-1][1]:\n",
    "        raise RuntimeError(f'This solution is not valid -- Start node {Path[1][1]} is not equal to end node {Path[-1][1]}')\n",
    "    print('Start and end node are the same. This is valid!')\n",
    "\n",
    "\n",
    "    print('Valid route!')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Read the results and analyze the path\n",
    "\n",
    "In reading the returned solution, the route is mapped into a human readable format.\n",
    "\n",
    "In analyzing the path, the validity of the route is checked by going over whether constraints are satisfied. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##### Call the function to interpret/convert/analyze the optimization results into a more meaningful/understandable format\n",
    "PathDict = ReadResults(route['configuration'], NodeName, CostMatrix, NumNodes)\n",
    "PathDict"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Well done! If the solver returned a valid solution, you have solved the traveling salesperson problem! If an error was raised then feel free to go back and adjust some settings and weights.\n",
    "\n",
    "## Next steps\n",
    "\n",
    "Now that you understand the problem scenario and how to define the cost function, there are a number of experiments you can perform to deepen your understanding and improve the solution defined above:\n",
    "\n",
    "- Modify the problem definition (e.g. by changing the number of nodes)\n",
    "- Rewrite the penalty functions to improve their efficiency\n",
    "- Tune the parameters (weights)\n",
    "- Try using a different solver, or a parameterized version (see [Which optimization solver should I use?](https://docs.microsoft.com/en-gb/azure/quantum/optimization-which-solver-should-you-use) for some tips)\n"
   ]
  }
 ],
 "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}