# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3.9.1 64-bit
#     name: python391jvsc74a57bd028dde4813e0e3ae4d0d68fcd2b1a59a20aa13de3dd42e33b497643f35eca055b
# ---

# # Traveling Salesperson Problem with Azure Quantum

# # Section 1 
#
# In the [Quantum Computing Foundations](/learn/paths/quantum-computing-fundamentals?azure-portal=true) learning path, you have helped the spaceship crew optimize its asteroid mining expeditions and reparations of critical onboard emergency systems. Now you have been asked to code the navigation system for the spaceship to optimize its travel routes through the solar system. The spaceship needs to visit numerous planets, moons, and asteroids to mine and sell the space rocks. In this module, you'll use the Azure Quantum service to minimize the travel distance of the spaceship. 
#
# The design and implementation of the navigation system is a case of the [**traveling salesperson** problem](https://en.wikipedia.org/wiki/Travelling_salesman_problem?azure-portal=true). The objective is to find a path through a network of nodes such that an associated cost, such as the travel time or distance, is as small as possible. You can find applications of the traveling salesperson problem, or slight modifications of it, in a variety of fields. For example, in logistics, chemical industries, control theory, and bioinformatics. Even in your daily life, when you want to know the order in which you should visit school, the cinema, the office, and the supermarket to go the shortest way.
#
# In this module, we will cover the formulation of this **minimization problem** by modeling travel costs and penalty functions. Afterward, we will solve the problem using the Earth's Azure Quantum Optimization service. All the content will be explained in the context of a spaceship that has to travel through the solar system to mine and sell asteroids.
#
# ## Learning objectives
#
# After completing this module, you will be able:
#
# - to understand the traveling salesperson problem.
# - to evaluate the problem complexity, solvers, and the tuning process. 
# - to model travel costs of the traveling salesperson.
# - to formulate problem constraints into a penalty functions.
# - to represent the minimization problem using Azure Quantum.
# - to use the Azure Quantum Optimization service to solve optimization problems.
# - to read and analyze results returned by the Azure Quantum solvers. 
#
#
# ## Prerequisites
#
# - The latest version of the [Python SDK for Azure Quantum](/azure/quantum/optimization-install-sdk?azure-portal=true)
# - [Jupyter Notebook](https://jupyter.org/install.html?azure-portal=true)
# - An Azure Quantum workspace
# - Basic linear algebra, only needing vector-matrix multiplications.
#
# If you don't have these tools yet, we recommend that you follow the [Get started with Azure Quantum](/learn/modules/get-started-azure-quantum/?azure-portal=true) module first.
#

# # Section 2
#
# For the module, we define the traveling salesperson problem for the spaceship as follows: 
#
# You have a set of $N$ locations (planets, moons, and asteroids) that need to be visited. As the newly appointed navigation engineer, you have to code the navigation system in order to find the shortest routes through the solar system. 
#
# There are a number of conditions that you will need to fullfill in order for the navigation system to be considered 'successful' by the crew. These will be considered the problem constraints:
#
# - Conditions:
#   - `One location at a time constraint` - The spaceship can only be in one planet, moon, or asteroid at a time. No magic!
#   - `No dissapearing constraint` - The spaceship has to be in a location, it can not disappear!
#   - `No revisiting locations constraint` - The spaceship may **not** visit a location more than once, except the homebase (start/end location). 
#   - `Start and end locations constraints` - The spaceship has to start and finish in the home base, which will be Mars.
#
# - Please consider the following points for the module as well:
#   - We will ignore orbital mechanics and use real approximate mean distances between the planets, moons, and asteroids. Distances are in millions of kilometers. All locations have a terrestrial surface (no gas planet included, as that is hard to mine).
#   - We will not specify where the space rocks can, and can not be, sold and mined. However, you may expand the constraints yourself to include such criteria.
#   - The $N$ locations are (in order) : Mars, Earth, Earth's Moon, Venus, Mercury, Ganymede (moon), Titan (moon), Ceres (asteroid), Pallas (asteroid), Cybele (asteroid). Naming the planets will make it easier to read the solutions returned by the solver!
#   
#
# ## About the problem
# The goal behind reprogramming the navigation system is to find a suitable interplanetary route to minimize the spaceship's travel distances. Having a (near)-optimal route will not only lower the total distance but will also reduce the ship's energy consumption and the number of working hours for the crew. An efficient route is crucial to the overall productivity of the spaceship's mining ambitions!
#
# The main idea is to create a **cost function**, which is used to model the travel costs **and** travel constraints of the spaceship. A particular route wants to be found that minimizes the cost function as much as possible. By incorporating the constraints into the cost function through penalties, the solver can be encouraged to settle for certain solution sets. Intuitively, you can imagine this as designing the rugged optimization landscape such that the 'good' solutions (routes) have lower cost values (minima) than 'bad' solutions. Even though the problem sounds easy to solve, for a large number of locations it becomes nearly impossible to find the optimal route (global minima). The difficulty is predominantly caused by:
# - the rugged optimization landscape (non-convex). 
# - the combinatorial explosion of the solution space, meaning that you may not be able to check every possible route.
# - the variables optimized for are non-continuous, $x_{v} \in \{0,1\}.$
#
# > [!TECHNICAL NOTES]
# > The traveling salesperson problem in which an optimal route needs to be found belongs to the set of [NP-hard problems](https://en.wikipedia.org/wiki/NP-hardness?azure-portal=true). First, there is no explicit algorithm for which you can find the optimal solution, meaning that you are faced with a search of the $(N-1)!$ (factorial number of routes to pass through all the locations) possible routes. Secondly, to verify whether a candidate solution to the traveling salesperson problem is optimal, you would have to compare it to all $(N-1)!$ other candidate solutions. Such a procedure would force you to compute all routes which is an **extremely** difficult task for large $N$ (non-polynomial time $\mathcal{O}((N-1)!)$).
# > The total number of routes for the traveling salesperson problem is dependent on the problem formulation. For generality, we will consider the directed version (`directed graph`) of the traveling salesperson problem, meaning that every route through the network is unique. Therefore you can assume that there are $(N-1)!$ different routes possible, since a starting location is given. 
#
# ## Azure Quantum workspace
#
# As presented in the previous two modules, [Solve optimization problems by using quantum-inspired optimization](/learn/modules/solve-quantum-inspired-optimization-problems?azure-portal=true) and [Solve a job scheduling optimization problem by using Azure Quantum](/learn/modules/solve-job-shop-optimization-azure-quantum?azure-portal=true), the optimization problem can be submitted to Azure Quantum solvers. For this, we will use the the Python SDK and format the traveling salesperson problem for the solver with cost function `terms`. 
# In order to submit the optimization problem to the Azure Quantum solver later, you will need to have an Azure Quantum workspace. Follow these [directions](/learn/modules/get-started-azure-quantum?azure-portal=true) to set one up if you don't have one already.
#
# Create a Python file or Jupyter Notebook, and be sure to fill in the details below to connect to the Azure Quantum solvers and import relevant Python modules. In case you don't recall your `subscription_id` or `resource_group`, you can find your Workspace's information on the Azure Portal page.
# Throughtout the module we will be appending to the Python script below such that you can run and view intermediary results. 
#
# ```python
#
# from azure.quantum import Workspace
# from azure.quantum.optimization import Problem, ProblemType, Term, HardwarePlatform, Solver
# from azure.quantum.optimization import SimulatedAnnealing, ParallelTempering, Tabu, QuantumMonteCarlo
# from typing import List
#
# import numpy as np
# import math
#
# workspace = Workspace (
#     subscription_id = "",  # Add your subscription id
#     resource_group = "",   # Add your resource group
#     name = "",             # Add your workspace name
#     location = ""          # Add your workspace location (for example, "westus")
# )
#
# workspace.login()
# ```
# The first time you run this code on your device, a window might prompt in your default browser asking for your credentials.
#
# ## Definitions
#
# A list to quickly look up variables and definitions:  
#
# - Locations: The planets, moons, and asteroids. 
# - Origin: A location the starship departs from.
# - Destination: A location the starship travels to. 
# - Trip: A travel between two locations. For example, traveling from Mars to Earth.
# - Route: Defined as multiple trips. Traveling between more than two locations. For example, Mars &rarr; Earth &rarr; Celeste &rarr; Venus.
# - $N$: The total number of locations. 
# - $C$: The travel cost matrix. In this module, the costs are distances given in millions of kilometers. 
# - $i$: Variable used to index the rows of the cost matrix, which represent origin locations.
# - $j$: Variable used to index the columns of the cost matrix, which represent destination locations. 
# - $c_{i,j}$: The elements of the matrix, which represent the travel cost, that is the distance, between origin $i$ and destination $j$.
# - $x_v$: The elements of a location vector. Each represents a location before or after a trip. These are the optimization variables.
# - $w_b$: Constraint weights in the cost function. These need to be tuned to find suitable solutions.
#
# ## Problem formulation
#
# With problem background dealt with, you can start looking at modeling it. In this unit, we will first go over how to calculate the spaceship's travel costs. This will become the foundation of the cost function. In later units, you will expand the cost function by incorporating penalty functions (constraints), to find more suitable routes throughout the solar system.
#
# Let's get to work!
#
# ### Defining the travel cost matrix $C$
#
# Consider a single trip for the spaceship, from one planet to another. The first step is to give each location (planet, moon, asteroid) a unique integer in $\{0,N-1\}$ (N in total). Traveling from planet $i$ to planet $j$ will then have a travel distance
# of $c_{i,j}$, where $i$ denotes the `origin` and $j$ the `destination`.
#
# $$ \text{The origin location } (\text{node}_i) \text{ with } i \in \{0,N-1\}.$$
# $$ \text{The destination location } (\text{node}_j) \text{ with } j \in \{0,N-1\}.$$
# $$ \text{Distance from } \text{location}_i \text{ to } \text{location}_j \text{ is } c_{i,j}.$$
#
# Writing out the travel distances for every $i$-$j$ (origin-destination) combination gives the travel cost matrix $C$:
#
# Travel cost matrix:
# $$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}. $$
#
# Here, the rows are indexed by $i$ and represent the origin locations. The columns are indexed by $j$ and represent the destination locations. For example, traveling from location $0$ to location $1$ is described by:
#
# $$ C(0,1) = c_{0,1}.$$ 
#
# The travel cost matrix is your spaceflight "dictionary", the most important asset to you as the spaceship's navigation engineer. It contains the most crucial information in finding a suitable route through the solar system to sell and mine space rocks. The measure of the travel cost is arbitrary. In this module the cost we want to minimize is the distance, however you can also use time, space debris, space money, ..., or a combination of different costs. 
#
# ### Defining the location vectors
#
# With a travel cost formulation between two locations complete, a representation of the origin and destination needs to be defined to select an element of the travel cost matrix. Remember, this is still for a single trip between location $i$ and location $j$. Selecting an element of the matrix can be achieved by multiplying the matrix with a vector from the left, and from the right! The left vector and right vector specify the origin and the destination, respectively. Consider the example where the spaceship travels from planet $0$ to planet $1$: 
#
# $$ \text{Travel cost - location 0 to location 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}.$$
#
# Recall that the spaceship can only be in one location at a time, and that there is only one spaceship, therefore only one element in the origin and destination vectors can equal 1, while the rest of them are 0. In other words, the sum of elements of the origin and destination vectors must equal 1. We will work this into a constraint in a later unit.
#
# Fantastic, you now know how to extract information from the travel cost matrix $C$ and express the travel cost for a single trip. However, it is necessary to express these ideas in a mathematical format for the solver. In the previous example the trip was hard-coded to go from location $0$ to location $1$. Let's generalize for any trip:
#
# $$ x_v \in \{0,1\} \text{ for } v \in \{0,2N-1\}, $$
#
# $$ \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}. $$
#
# Each $x_v$ represents a location before or after a trip. The reason that the destination vector elements are indexed from $N+1$ to $2N$ instead of $0$ to $N-1$ is to differentiate the origin variables and destination variables for the solver. If this is not done, the solver would consider the origin and destination the same, meaning that the spaceship would never take off to its destination planet! When submitting the optimization problem, the solver will determine which $x_v$ is assigned a value 1 or 0, dependent on the respective trip's travel cost. To re-iterate an important point, the sum of the origin and destination vector elements must be 1 as the spaceship can not be in more than or less than one location at a time: 
#
# $$ \text{Sum of the origin vector elements: }\sum_{v = 0}^{N-1} x_v = 1, $$ 
# $$ \text{Sum of the destination vector elements: }\sum_{v = N}^{2N-1} x_v = 1.$$ 
#
#
# Now that the basic mathematical formulations are covered, the scope can be expanded. As the spaceship's navigation engineer, it is your job to calculate a route between multiple planets, moons, and asteroids, to help make the space rock mining endeavour more succesful. With the travel cost matrix and the origin/destination vectors defined, you are equipped with the right mathematical tools to start looking at routes through the solar system. Let's take a look at how two trips can be modeled!
#
#
# ### Defining the travel costs for a route
#
# To derive the travel costs for two trips, which can be considered a route, you will need a way to describe the **total** travel cost. As you might expect, the total cost of a route is simply the sum of the trip's travel costs that constitute the route (sum of the trips). Say you have a route $R$ in which the spacehip travels from location $1$ to location $3$ to location $2$. Then the total cost is the sum of the trips' costs:
#
# $$ \text{Cost of route: } R_{1-3-2} = c_{1,3}+c_{3,2}. $$
#
# By taking a closer look at the equation a very important relation is found. The destination for the first trip and the origin for the second trip are the same. Incorporating this relation in the cost function will help reduce the number of variables the solver has to optimize for. As a result of the simplification, the optimization problem becomes much easier and has an increased probability of finding suitable solutions! 
#
# > [!TECHNICAL NOTE]
# > For the 2-trip example this recurrence relation reduces the solver search space size (number of possible routes/solutions) from $2^N \cdot 2^N \cdot 2^N \cdot 2^N$ to $2^N \cdot 2^N \cdot 2^N$, a factor $N$ ($N$ = 10 in this module) difference. Without considering constraints, each vector element can take two values (${0,1}$) and has a length $N$, therefore each vector multiplies the solution set by $2^N$. The reduction in variables becomes even more apparent for longer routes. Visiting the 10 planets, moons, and asteroids as in this module, the relation reduces the search space size from $2^{20N}$ to $2^{11N}$! The number of variables optimized for decreases from $200$ to $110$, respectively.
#
# Recall that the travel costs can be written with vectors and matrices. Then the total travel cost of the route is: 
#
# $$ \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}.$$
#
# In the equation the destination vector of the first trip and the origin vector of the second trip contain the same variables, making use of the recurrence relation. Generalizing the 2-trip route to a $N$-trip route is achieved by including more vector-matrix multiplications in the addition. Written out as a sum:
#
# $$\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),$$
#
# which can equivalently be written as:
#  
# $$\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),$$
#
# where the $x$ variables indices are dependent on the trip number $k$, the total number of locations $N$, and the origin ($i$) and destination ($j$) locations. 
#
# Great! A function to calculate the travel costs for the spaceship's route has been found! Because you want to minimize (denoted by the 'min') the total travel cost with respect to the variables $x_v$ (written underneath the 'min'), we will write the foundation of the cost function as follows:
#
# $$\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).$$
#
# This equation will not make up the entire cost function for the solver. Penalty functions have to added to it for the constraints, otherwise the solver will return invalid solutions. These will be added to the cost function in upcoming units. 
#
# The crew is very excited with the progress you have already booked. You are ready to program the travel costs into navigation system you are building for the spaceship. Time to write some code!
#
#
# ### Progress on the navigation system's cost function
#
# - The cost function contains:
#   - The `travel costs`.
#
# $$ \text{Cost function so far: } \underset{x_0, x_1,\dots,x_{(N^2+2N)}}{min} \left(\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)\right) $$
#
# ### Coding the travel costs 
#
# For the solvers to find a suitable solution, you will need to specify how it should calculate the travel cost for a 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 above, the weighting is simply the $c_{i,j}$-th element of the cost matrix, resembling the distance between the locations. For example, weighting a trip from location $1$ to location $2$ in the second trip ($k=1$) has the following weighting:
#
# $$ c_{i,j} \cdot x_{Nk+i} \cdot x_{N(k+1)+j} =  c_{1,2} \cdot x_{N+1} \cdot x_{N+2}. $$
#
# The $x_v$ variables have an $N$ term because for the solver we need to differentiate between the variables over the trips. In other words, after each trip there are $N$ new variables that represent where the spaceship can travel to next.
#
# Below, you can find the code snippet that will be added to the Python script. We define a problem instance through the `OptProblem` function, in which we will continue to append pieces of code throughout the module's units. Later, this problem will be submitted to the Azure solvers. If you want to see how the weights are assigned to each $ x_{Nk+i} \cdot x_{N(k+1)+j}$ combination, you can uncomment the last lines in the code. Lastly, since the optimization variables $x_v$ can take values ${0,1}$, the problem type falls into the `pubo` category (polynomial unconstrained binary optimization). 
#
# ```python
#
# ##### Define variables
#
# # The number of planets/moons/asteroids.
# NumLocations = 10
#
# # Location names. Names of some of the solar system's planets/moons/asteroids.
# LocationNames = {0:'Mars', 1:'Earth', 2:"Earth's Moon", 3:'Venus', 4:'Mercury', 5:'Ganymede', 6:'Titan', 7:'Ceres', 8:'Pallas', 9:'Cybele'}
#
# # Approximate mean distances between the planets/moons/asteroids. Note that they can be very innacurate as orbital mechanics are ignored. 
# # This is a symmetric matrix since we assume distance between planets is constant for this module.
# CostMatrix = np.array([     [0,   78,       2,  120,   170,  550,  1200,  184, 600,  1.5   ],
#                             [78,   0,     0.5,   41,    92,  640,  1222,  264, 690,  0.25  ],
#                             [2,   0.5,      0,   40,    91,  639,  1221,  263, 689,  0.25  ], 
#                             [120,  41,     40,    0,    50,  670,  1320,  300, 730,  41.5  ],
#                             [170,  92,     91,   50,     0,  720,  1420,  400, 830,  141.5 ],
#                             [550,  640,   639,  670,   720,    0,   650,  363,  50,  548   ],
#                             [1200, 1222, 1221, 1320,  1420,  650,     0, 1014,  25,  625   ],  
#                             [184,  264,   263,  300,   400,  363,  1014,    0, 100,  400   ],
#                             [600,  690,   689,  730,   830,   50,    25, 100,    0,  350   ],
#                             [1.5,  0.25, 0.25, 41.5, 141.5,  548,   625, 400,  350,  0     ]
#                       ])    
#                        
# ##### If you want try running with a random cost matrix, uncomment the following:
# #maxCost = 10
# #CostMatrix = np.random.randint(maxCost, size=(NumLocations,NumLocations))
#  
# ############################################################################################
# ##### Define the optimization problem for the Quantum Inspired Solver
# def OptProblem(CostMatrix) -> Problem:
#     
#
#     #'terms' will contain the weighting terms for the solver!
#     terms = []
#
#     ############################################################################################
#     ##### Cost of traveling between locations  
#     for k in range(0,len(CostMatrix)):                          # For each trip (there are N trips to pass through all the locations and return to home base)
#         for i in range(0,len(CostMatrix)):                      # For each origin (reference location)
#             for j in range(0,len(CostMatrix)):                  # For each destination (next location w.r.t reference location)
#                 
#                 #Assign a weight to every possible trip from location i to location j - for any combination
#                 terms.append(
#                     Term(
#                         c = CostMatrix.item((i,j)),                                     # Element of the cost matrix
#                         indices = [i+(len(CostMatrix)*k), j+(len(CostMatrix)*(k+1))]    # +1 to denote dependence on next location
#                     )
#                 )
#                 ##----- Uncomment one of the below statements if you want to see how the weights are assigned! -------------------------------------------------------------------------------------------------
#                 #print(f'{i+(len(CostMatrix)*k)}, {j+(len(CostMatrix)*(k+1))}')                                                                 # Combinations between the origin and destination locations 
#                 #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)
#                 #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
#
#     return Problem(name="Spaceship navigation system", problem_type=ProblemType.pubo, terms=terms)
#
# OptimizationProblem = OptProblem(CostMatrix)
#
# ``` 
#
#
# ## Next Steps
#
# As navigation engineer of the spaceship you have completed the initial work on the route planner. But with the current code, the Azure Quantum solvers on Earth would send routes to the spaceship that would make no sense! Currently, all $x_v$ would be assigned a value 0 by the solver because that would yield a value-0 travel cost! To avoid this from happening, you will need to implement the constraints into the cost function. These constraints will penalize the routes that are incorrect. More on this in the next units! 
#
#

# # Section 3
#
# With the foundations of the navigation system in place, it is time to incorporate the constraints. The initial results made the crew enthousiastic! However, before launching the final system, you will need to implement penalty functions to ensure that the Azure Quantum solvers generate feasible routes through the solar system. 
#
# In this unit we will consider the *"one location at a time"* constraint. Recall that the spaceship can only be in one location at a time. It is impossible for the spasceship to be on both Earth and Mars simultaneously! You will need to program this into the navigation system by penalizing the solvers for routes in which such magic happens. 
#
# >[!Note]
# >From this unit onwards, the origin and destination vectors will be referred to as **location vectors** due to the recurrence relation presented in the previous unit. The location vector is used to represent the location of the spaceship after a number of trips. For example, starting from the home base will give the location vector 0. After completing the first trip, the spaceship will be at a location given by location vector 1. 
#
#
# ## "One location at a time" constraint
#
# The spaceship can only be at one planet, moon, or asteroid at a time. So far, the defined cost function only contains information about the travel costs, nothing about whether the spaceship is at multiple locations at once. To mathematically express this idea, it must be true that only one element in a location vector equals 1. For example:
#
# $$ \text{Incorrect location vector: } \begin{bmatrix} 1 \\ 1 \\ 0 \\ \vdots \\ 0 \end{bmatrix}, $$
#
# should not be allowed as the spaceship is at the first two locations simultanesouly. Only one element may be non-zero:
#
# $$ \text{Correct location vector: } \begin{bmatrix} 0 \\ 1 \\ 0 \\ \vdots \\ 0 \end{bmatrix}. $$
#
# One way of designing this constraint would be look at the sum of the location vector elements, which must equal 1: 
#
# $$ \text{Location vector 0: (home base)} \hspace{0.5cm}  x_0 + x_1 + \dots + x_{N-1} = 1, $$
# $$ \text{Location vector 1: } \hspace{0.5cm} x_{N} + x_{N+1} + \dots + x_{2N-1} = 1, $$
# $$ \text{Location vector N (home base): } \hspace{0.5cm} x_{N^2} + x_{N^2+1} + \dots + x_{N^2 + N-1} = 1. $$
#
# Enforcing this constraint over all trips would then give (N+1 because the spaceship returns to the home base on Mars):
#
# $$ \text{For all locations: } \hspace{0.5cm}  x_0 + x_1 + \dots +  x_{N^2 + N-1} = N+1. $$
#
# This equation is a valid way to model the constraint. However, there is a drawback to it as well. In the equation, only the individual locations are penalized. There is no information about being in two locatoins at once! The following values satisfy the equation:
#
# $$ \text{If } x_0=1, x_{1}=1, \dots, x_{N}=1, \text{ and the remaining } x_{N+1}=0, x_{N+3}=0, \dots, x_{N^2+N-1}=0,$$  
#
# but violate the constraint we are trying to model. With these values, the spaceship is at all planets/moons/asteroids at once before the first trip (location vector 0). Let's rethink. Consider a small example of three locations (a length-3 location vector). Then if the spaceship is in location 0, it can not be in location $1$ or $2$. Instead of using a summation to express this, another valid way would be to use products. The product of elements of a location vector must always equal zero, regardless where the spaceship is, because only one of the three $x_v$ elements can take value 1. Writing out the products of a single location vector gives:
#
# $$ x_0 \cdot x_1 = 0,$$
# $$ x_0 \cdot x_2 = 0,$$ 
# $$ x_1 \cdot x_2 = 0.$$
#
# In this format, the constraint is much more specific and stringent for the solver. These equations reflect the interrelationships between the locations more accurately than the summation. They can be implemented in a way that distinguish the variables for different location vectors, unlike the summation which adds all the variables of the location vectors together. As a result of describing the constraint more specifically, the solver will return better solutions. 
#
# >[!Note]
# >We do not want to weight combinations between locations more than once, as this would lead to inbalances (assymmetries) in the cost function. Tuning the weights of an imbalanced cost function tends to be more difficult. We therefore exclude the reverse combinations:
# >$$ x_{1}\cdot x_{0}, $$
# >$$ x_{2}\cdot x_{0}, $$
# >$$ x_{2}\cdot x_{1}. $$
#
# The next step consists on generalizing the one location at a time constraint for the spaceship, which has to pass by all locations in the solar system and return to home base (N+1 locations need to be visited). In the equation below, $l$ iterates over the number of location vectors, while $i$ and $j$ iterate over the vector elements:
#
# $$ \text{One location at a time constraint: }0=\sum_{l=0}^{N} \sum_{i=0}^{N-1} \sum_{j=0}^{N-1} x_{(i+Nl)} \cdot x_{(j+Nl)} \text{ with } \{ i,j | i<j \} $$
#
#
# ### Progress on the navigation system's cost function
#
# - The cost function contains:
#   - The `travel costs`.
#   - The `one location at a time constraint`, with constraint weight $w_1$. 
#
# $$ \text{Cost function so far: } \underset{x_0, x_1,\dots,x_{(N^2+2N)}}{min} \left(\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)  +   w_1 \left( \sum_{l=0}^{N} \sum_{i=0}^{N-1} \sum_{j=0}^{N-1} x_{(i+Nl)} \cdot x_{(j+Nl)} \text{ with } \{ i,j | i<j \} \right) \right) $$
#
# ## Coding the constraint
#
# Great! With the mathematics written out, it is time to add the constraint into the navigation system's code. As in previous modules, you'll need to assign constraint weights to penalize invalid solutions. For now, you can ignore them as they need to be tuned after coding the remaining constraints. Weight tuning will require a bit hands-on experimentation, as they will impact the satisfiability of the route returned by the solver.
#
#
# ``` python
#
#     ############################################################################################
#     ##### Constraint: One location at a time constraint - spaceship can only be in 1 location at a time.
#     for l in range(0,len(CostMatrix)+1):                # The total number of locations that are visited over the route (N+1 because returning to home base)
#         for i in range(0,len(CostMatrix)):              # For each location (iterate over the location vector)
#             for j in range(0,len(CostMatrix)):          # For each location (iterate over the location vector)
#                 if i!=j and i<j:                        # i<j because we don't want to penalize twice // i==j is forbidden (this could equal 1, that's why)
#                     terms.append(
#                         Term(
#                             c = int(0),                                    
#                             indices = [i+(len(CostMatrix)*l),j+(len(CostMatrix)*l)]                   
#                         )
#                     )
#                     ##----- Uncomment one of the below statements if you want to see how the weights are assigned! -------------------------------------------------------------------------------------------------
#                     #print(f'{i+(len(CostMatrix)*l)},{j+(len(CostMatrix)*(l))}')
#                     #print(f'Location constraint 1: x_{i+(len(CostMatrix)*l)} - x_{j+(len(CostMatrix)*(l+1))} (trip {l}) assigned weight: 0')  # In a format for the solver (as formulated in the cost function)
#
#
#     return Problem(name="Spaceship navigation system", problem_type=ProblemType.pubo, terms=terms)
#
# OptimizationProblem = OptProblem(CostMatrix)  
#
#
# ``` 
#
#
# ## Next steps
#
# The first constraint has been finished. Significant progress is being made on the spaceship's navigation system! Let's keep going and work on the implementing the second constraint, the spaceship may not disappear in the navigation system! 
#
#
#
#

# # Section 4
#
# Unfortunately, the navigation system so far is not really helpful. Due to the way we formulated the cost function, the spaceship is penalized for being 'anywhere' without regard to the mission objective of visiting each planet, moon and asteroid of the itinerary.
# Recall that the cost function includes two functions so far, the cost for traveling, and the one location at a time constraint. Currently, the Azure solvers on Earth would consistently return location vectors containing only 0 elements because that would lead to no incurred travel costs and no constraint violation, even though the spaceship is 'nowhere'. 
#
# > By looking at the code and previous sections you come to the conclusion that the minimal cost is obtained by setting all $x_v$ to 0.
#
# It is your task to make sure that the spaceship does not dissapear in the navigation system. In this unit, we will cover the *"no dissapearing"* constraint. By incorporating negative penalty terms (rewards) in the cost function, the solver will be encouraged to return routes in which the spaceship will consistently be in some location in our solar system. These negative cost terms will locally decrease cost function value(s) for suitable solution configurations. The aim is that these function valleys in the optimization landscape become low enough for the solver to settle in, depending on the constraints and weighting. 
#
#
# ## "No dissapearing" constraint
#
# By incorporating negative terms in the cost function, we can encourage the solver to return a particular set of solutions. To demonstrate this point, take the following optimization problem: 
#
# $$f(\boldsymbol{x}) := \underset{x_0, x_1, x_2}{min} x_0 + x_1 - x_2,$$ 
# $$\text{with } x_0, x_1, x_2 \in \{0,1\}.$$
#
# The minimum value for this example is achieved for the solution $x_0$ = 0, $x_1$ = 0, $x_2$ = 1, with the optimal function value equal to -1. Here the negatively weighted term encourages $x_2$ to take a value 1 instead of 0, unlike $x_0$ and $x_1$. With this idea in mind, you can prevent the spaceship from dissapearing in the navigation system. 
#
# If the spaceship has to visit $N+1$ locations (+1 because return to home base) for a mining expedition, then $N+1$ of the $x_v$ should be assigned a value 1. Written as an equation for the $N(N+1)$ variables of a route ($N+1$ vectors, each with with $N$ elements gives $N(N+1)$):
#
# $$ \sum_{v=0}^{N(N+1)-1} x_v = N+1.$$
#
# You could split this equation for each location vector separately, but if you keep the constraint linear in $x_k$ then the resulting cost function will be the same. Moving the variables to the right side of the equation for negative weighting gives:
#
# $$ 0 = (N+1) -\left( \sum_{v=0}^{N(N+1)-1} x_v \right).$$
#
# As you may be aware, there is no guarantee which particular $x_v$ the solver will assign a value 1 in this equation. However, in the previous constraint the spaceship is already forced to be in a maximum of one location at a time. Therefore, with the cost function weights tuned, it can be assumed that only one $x_v$ per location vector will be assigned a value 1. 
#
# To summarize, the `one location at a time constraint` penalizes the solver for being in more than one node at a time, while the `no dissapearing constraint` rewards the solver for being at as many locations as possible. The weights of the constraints will effectively determine how they are satisfied, a balance between the two needs to be found such that both are adhered to.
#
# >[!Note]
# > You can combine the `one location at a time constraint` and `no dissapearing constraint`. For explanatory purposes these two are split, however if you realize that -$x_v^2$ = -$x_v$, then you can expand the 'if' statement in the previous constraint for 'i==j'. 
#
# ### Progress on the navigation system's cost function
#
# - The cost function contains:
#   - The `travel costs`.
#   - The `one location at a time constraint`, with constraint weight $w_1$. 
#   - The `no dissapearing constraint`, with constraint weight $w_2$.
#
# $$ \text{Cost function so far: } \underset{x_0, x_1,\dots,x_{(N^2+2N)}}{min} \left(\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)  +   w_1 \left( \sum_{l=0}^{N} \sum_{i=0}^{N-1} \sum_{j=0}^{N-1} x_{(i+Nl)} \cdot x_{(j+Nl)} \text{ with } \{ i,j | i<j \} \right) +\left( \sum_{v=0}^{N(N+1)-1} x_v \right) \right).$$
#
#
# ## Coding the constraint
#
# Incorporating this constraint in the spaceship's navigation system's code is straightforward. We simply need to weight each variable by a negative weight. The $N+1$ can be ignored for the cost function, as it is independent of the variables optimized for and because it has no effect on the returned results. This scalar term is a constant linear offset in the cost for all solution configurations, and you can visualize it as a upward/downward shift of the entire optimization landscape.
#
# ``` python
#
#     ############################################################################################
#     ##### Constraint: No dissapearing constraint - encourage the spaceship to be 'somewhere' otherwise all x_k might be 0 (for example).
#     for v in range(0, len(CostMatrix) + len(CostMatrix) * (len(CostMatrix))):    # Select variable 
#         terms.append(
#             Term(
#                 c = -int(0),          
#                 indices = [v]   
#             )
#         )
#         ##----- Uncomment one of the below statements if you want to see how the weights are assigned! -------------------------------------------------------------------------------------------------
#         #print(v)
#         #print(f'No dissapearing constraint 2: x_{v} assigned weight: 0')                                                      # In a format for the solver (as formulated in the cost function)
#         #print(f'No dissapearing constraint 2: location_{v % NumLocations} after {np.floor(v / NumLocations)} trips assigned weight: 0')   # In a format that is easier to read for a human
#
#     return Problem(name="Spaceship navigation system", problem_type=ProblemType.pubo, terms=terms)
#
# OptimizationProblem = OptProblem(CostMatrix)  
#
#
# ``` 
#
# ## Next steps
#
# By now, the navigation system should prevent the spaceship from being in multiple locations at once, prevent it from dissapearing, and select a low-cost travel route. Great work! Yet, for the mining mission it is crucial that no planets, moons, or asteroids are not visited more than once. That would be bad for business! As you can see from the travel cost matrix, it is also the case that staying in the same planet has a zero-weighted cost (the diagonal terms of the matrix). The time has come for you to start looking at route specific constraints. 
#
#
#
#
#

# # Section 5
#
# In this unit, you will be making sure that the spaceship does not visit any planet, moon, or asteroid more than once (except the home base). Space rocks are crucial to the construction and development of the newest scientific gadgets. Therefore a key part of mining expidition is also selling these rocks to the science and astronaut communities in our solar system. On each expedition, you want to visit each location to sell and mine the space rocks. The *"no revisiting locations"* constraint will penalize routes that revisit any planets, moons, and asteroids. 
#
# > Currently, the solver may provide solutions in which the spaceship revisits a location. For example, there is no cost associated with staying in the same location, as shown by the diagonal elements of the cost matrix. It is necessary to penalize routes in which revisits occur. 
#
#
# ## "No revisiting locations" constraint
#
# The navigation system may generate routes in which a location is visited more than once. Similarly to previous modules, we will have to penalize routes in which this occurs. To start, let's consider the fourth location $x_3$ (Venus) for each trip:
#
# $$ x_3, x_{3+N}, x_{3+2N}, ... $$
#
# For a valid route, the spaceship will only have passed through this location once. That means that only one of these variables is allowed to take a value 1, the remaining have to be 0. For example, if the spaceship visits Venus after the first trip, then $x_{3+N}=1$ and $x_3=0$, x_{3+2N}=0$, and so on. As done similarly with the one location at a time constraint, the product of these variables must equal 0. The following can then be derived:
#
# $$ x_3 \cdot x_{3+N} \cdot x_{3+2N} \cdot ... = 0. $$
#
# Even though this constraint might seem correct, it is not stringent enough. If only one of these variables is 0, then the constraint is already satisfied. However, the multiplication of variables can be split into their individual products:
#
# $$ x_3 \cdot x_{3+N} =0,$$
# $$ x_3 \cdot x_{3+2N} = 0,$$
# $$ x_{3+N} \cdot x_{3+2N} = 0,$$
# $$ \vdots $$
#
# By weighting these terms in the cost function, at least one of the two variables for each respective equation has to be zero. The solver will therefore be penalized if more than one variable takes the value 1. Continuing this relation for all $N$ trips and all $N$ locations yields the penalty function:
#
# $$ 0 = \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),$$
#
# in which the first summation ($p$) assigns a reference location $x_p$, and the second summation assigns the same location but after a number of trips $x_f$ (multiple of $N$ due to the stepsize). 
#
#
# ### Progress on the navigation system's cost function
#
# - The cost function contains:
#   - The `travel costs`.
#   - The `one location at a time constraint`, with constraint weight $w_1$. 
#   - The `no dissapearing constraint`, with constraint weight $w_2$.
#   - The `no revisiting locations constraint`, with constraint weight $w_3$.
#
# $$ \text{Cost function so far: } \underset{x_0, x_1,\dots,x_{(N^2+2N)}}{min} \left(\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)  +   w_1 \left( \sum_{l=0}^{N} \sum_{i=0}^{N-1} \sum_{j=0}^{N-1} x_{(i+Nl)} \cdot x_{(j+Nl)} \text{ with } \{ i,j | i<j \} \right) +w_2\left( \sum_{k=0}^{N(N+1)-1} x_k \right) + w_3 \left(  \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)   \right)          \right) $$
#
#
# ## Coding the constraint
# Integrating this into the cost function will penalize routes in which locations are visited more than once. The derived model does not 
# penalize the last trip in which the spaceship should return to the home base. The reason is that the last trip ($N$-th trip) would always violate the constraint, as the $N$ different locations have all already been visited. Also, adding it would only make the cost function larger without adding any value to the optimization problem. In the next unit, you will look at how the spaceship can end at a specific location.
#
#
# ``` python
#
#     ############################################################################################                        
#     ##### Constraint: no revisiting locations constraint --- (in the last step we can travel without penalties (this is to make it easier to specify an end location ))
#     for p in range(0,len(CostMatrix)+len(CostMatrix)*(len(CostMatrix))):                                  # This selects a present location x: 'p' for present    
#         for f in range(p+len(CostMatrix),len(CostMatrix)*(len(CostMatrix)),len(CostMatrix)):              # This selects the same location x but after upcoming trips: 'f' for future
#             terms.append(
#                 Term(
#                     c =int(0),          
#                     indices = [p,f]   
#                 )
#             )     
#             ##----- Uncomment one of the below statements if you want to see how the weights are assigned! -------------------------------------------------------------------------------------------------
#             #print(f'x_{p},x_{f}')  # Just variable numbers 
#             #print(f'Visit once constraint: x_{p} - x_{f}  assigned weight: 0')  # In a format for the solver (as formulated in the cost function)
#             #print(f' Visit once constraint: location_{p%NumLocations} - location_{(p+f)%NumLocations} after {(f-p)/NumLocations} trips assigned weight: 0')  # In a format that is easier to read for a human
#
#
#     return Problem(name="Spaceship navigation system", problem_type=ProblemType.pubo, terms=terms)
#
# OptimizationProblem = OptProblem(CostMatrix)  
#
#
# ``` 
#
# ## Next steps
#
# Wow! The navigation system is nearly complete. Shortly, you will be able to assist the crew in their logistics and planning departments! There are still some final steps to complete before you are ready to demonstrate it. You still need a way to include some start and end locations and complete the weight tuning steps. In the next unit, the start and end locations will be added to the cost function. 
#

# # Section 6
#
# The crew of the spaceship has decided to start and end the mining expedition on their home base on Mars. You will need to code these conditions into the navigation system. Fortunately, relative to the other constraints, these are very easy to include! The only elements you will need to add to the navigation system's code are an initial location weighting, and a final location weighting. As done previously, the solver can be rewarded by finding routes which start and end on Mars. 
#
#
# ## "Start and end location" constraint 
# All that has to be done in this unit is assing negative costs for starting and ending in a specific location. Departure from the home base on Mars is represented by $x_0=1$. Returning to home base after $N$ trips is given by x_{N^2}=1. Assigning these specific variables negative weights in the cost function will promote the solver to give them that value. 
#
# This weighting procedure can also used to integrate user knowledge into the cost function. For example, visiting Venus after the $5$-th trip can be accomplished by negatively weighting $x_{3+5N}$. Alternatively, the weighting procedure can be used to promote a particular set of locations after a trip! For example, if you want to include a constraint about where the spaceship can and can't mine or sell after a trip, then you can negatively weight those set of locations. 
#
# In the second unit it was stated that we want to keep the cost function "balanced", avoiding re-weighting by reverse combinations and giving all constraint terms equal weights (currently all have 0 weights). The benefit is that it is easier to specify preferred locations here, without having to tune each variable individually, as their relative weightings to the cost function are all the same. This will also reduce the effort in tuning in upcoming units. 
#
# >[Note!]
# >The negative weights assigned for these variables here will contribute to the negative weights assigned to the respective variables in the `no dissapearing constraint`. As a best-practice, you can combine these terms to reduce the length of the cost function. 
#
# Another way of enforcing the `start and end location constraint` is to directly tell the solver to assign the respective variables a value equal to 1. The [`Problem.set_fixed_variables`](/azure/quantum/optimization-problem) is handy when you know what values certain variables **have** to be. Also, the cost function is automatically simplified for you when calling the method as the variables are filled in. The reason why the method is not implemented in this module is because you would not be able to encourage the spaceship to choose between a set of locations for a trip, a consequence of hard-coding the variable values. Even though promoting a set of locations is not part of the problem statement, it is an easy expansion of the `start and end location constraint` considered in this unit. 
#
#
# ### Progress on the navigation system's cost function
#
# - The cost function contains:
#   - The `travel costs`.
#   - The `one location at a time constraint`, with constraint weight $w_1$. 
#   - The `no dissapearing constraint`, with constraint weight $w_2$.
#   - The `no revisiting locations constraint`, with constraint weight $w_3$.
#   - The `start and end locations constraints`, with constraint weights $w_4$ and $w_5$. 
#
# $$ \text{Final cost function: } \underset{x_0, x_1,\dots,x_{(N^2+2N)}}{min} \left(\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)  +   w_1 \left( \sum_{l=0}^{N} \sum_{i=0}^{N-1} \sum_{j=0}^{N-1} x_{(i+Nl)} \cdot x_{(j+Nl)} \text{ with } \{ i,j | i<j \} \right) +w_2\left( \sum_{k=0}^{N(N+1)-1} x_k \right) + w_3 \left(  \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) \right) + w_4 ( x_0)  + w_5 ( x_{N^2}) \right)  $$
#
# ## Coding the initial conditions
#
# Here, we add negative weights to the start and end location variables. A route in which these locations are the start and finish location should have the lowest associated costs. 
#
# ``` python
#
#     #############################################################################################                        
#     ##### Begin at x_0 (Mars - home base)
#     terms.append(
#         Term(
#             c = -int(0),
#             indices = [0]   
#         )
#     )
#
#     ############################################################################################                        
#     ##### End at x_{N^2} (Mars - home base)
#     terms.append(
#         Term(
#             c = -int(0),
#             indices = [len(CostMatrix)*(len(CostMatrix))]   
#         )    
#     )
#
#
#     return Problem(name="Spaceship navigation system", problem_type=ProblemType.pubo, terms=terms)
#
# OptimizationProblem = OptProblem(CostMatrix)  
#
#
# ``` 
#
#
# ## Next steps
#
# Well done! You have finished programming the constraints for navigation system. It should not be long before the spaceship crew can use it to find their way through the solar system. In the next units we will go over constraint weight tuning and submitting the minimization problem to the Azure solvers. 
#
#

# # Section 7
#
# The minimization problem is ready to be submitted to the solvers in the Earth's Azure Quantum workspace. But before that, the navigation system now requires some fine-tuning and solver selection. In this unit you will learn how to read and analyze solutions returned by the solver, and about some differences in the available solvers. Make sure the crew buckles up, interplanetary routes will be available soon!
#
# ## Some notes about solvers
#
# There are numerous solvers (algorithms) and targets (hardware providers) available in the Azure Quantum optimization service. Choosing a which to use is not straightforward, it is recommended to review the [Microsoft QIO provider](/azure/quantum/provider-microsoft-qio?azure-portal=true) and [1QBit provider](/azure/quantum/provider-1qbit?azure-portal=true) documentation pages for further information and references. We will briefly cover some background of the solvers which might help you decide how you would like to submit the optimization problem. 
#
# Unfortunately, there is no standardized way of finding the best solver for (highly) non-convex problems. For these optimization problems you are left with two options. First, you can use `heuristic solvers`. Heuristic methods can be understood as searches of the optimization landscape in which as many minima want to be found as possible. The second option is to use linearizations of the problem, suitable for small and 'easy' non-convex problems (`constraint relaxation`, `cost function linearization`, etc.). For linearizations, a degree of accuracy is paid to reduce the problem complexity. However, for large non-convex problems linearizations usually perform poorly are as they generate **local** and **approximate** solutions. (Approximating nonlinearies with linear functions becomes increasingly difficult with the degree of nonlinearities). 
#
# For the traveling salesperson problem, it is practically required to use heuristic solvers due to the combinatioral explosion with the number of locations. The main bottleneck is how to efficiently design/weight the cost function and how to choose the right solver. Choosing an appropriate solver requires some research as each heuristic method is tailored for specific problem conditions. It can therefore be that a solver performs poorly on a certain problem instance, but well on another. Usually, finding a suitable one requires experimentation, and/or more thorough insight into the problem. To help you with that, you can consider the following notes:
#
# 1. What are the problem properties? 
#   - Is the problem convex or non-convex (rugged)? Is the optimization landscape slightly rugged(hills), or very rugged (peaks)?
#   - How large is the optimization problem (given by the number of variables $x_v$)?  
#   - Are you looking at an Ising model ($x_v \in \{-1,1\}$ ), or a PUBO/QUBO model ($x_v \in \{0,1\}$)?
#   - Can you reduce or simplify the optimization problem? Smaller cost functions will reduce computational effort for the solvers, and save time!
# 2. Finding solvers through experimentation. 
#   - Try a solver that you can use as a reference to compare to other solver performances, such as accuracy or time. For example, simulated annealing generally does well on a wide range of problems and can provide an initial metric to evaluate other solvers.
#   - Also, in early stages it can be helpful to use a solver for which parameters do not need to be tuned (`parameter-free`).
#   - Select solvers which match well with the problem, for an overview see the documentation [solvers and targets](/azure/quantum/qio-target-list?azure-portal=true).  
#   - Try other solvers and compare their performances. 
#   - Do all solvers perform poorly? Consider redesigning the cost function or a longer solver time (`timeout`). Often you can make the cost function more efficient for the solver, some tips are described in the [documentation](/azure/quantum/provider-microsoft-qio?azure-portal=true).
#   - Are there some good results? Consider these solvers as candidates for the fine-tuning process. Fine-tuning can be a lot of work, therefore it is useful to have a small subset of solvers. 
#   - Solver selection is an iterative process. If possible, it can be helpful to experiment on small-scale problem as this might be more insightful and reduce effort. 
# 3. Fine-tuning candidate solvers.
#   - Improve the performances by defining solver parameters such as, start conditions, convergence criteria, sweeps, etc. 
#   - Depending on the scale of the problem and solver type, you can consider different hardware (CPUs / FPGAs / GPUs / Annealers) and targets. 
#   - Adjust the cost function weights specifically for a type of solver.
#   - Solver tuning is also an iterative process.
#
# ## Submitting the optimization problem to Azure Quantum
#
# Submitting the navigation system's cost function to the Azure solvers will require a solver. As navigation engineer, it is important to know advantages and disadvantages of these methods, and for testing reasons, you want to give them all a go. 
#
# Below, you can find the code snippet that will be added to the python script. The solver hardware is submitted to the Microsoft QIO target and defaulted to CPU hardware. The 120 seconds for `timeout` defines when the optimization procedure should terminate (approximately). The four available solvers are `simulated annealing`, `parallel tempering`, `tabu search`, and `quantum Monte Carlo`. The first three are available parameter-free, meaning that they are tuned/inferred for you by the Azure Quantum solvers. For the `quantum Monte Carlo` algorithm this is not yet available, and therefore it has some necessary parameters predefined. 
#
# ``` python
#
# ############################################################################################                        
# ##### Define the solver --- uncomment if you wish to use a different one --- you can also change the timeout.
#
# solver = SimulatedAnnealing(workspace, timeout = 120)   
# #solver = ParallelTempering(workspace, timeout = 120)
# #solver = Tabu(workspace, timeout = 120)
# #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
#
#
# route = solver.optimize(OptimizationProblem)     # Solve the optimization problem -- wait until done.
# print(route)
#
#
# ```
#
# ## Reading and analyzing the results
#
# You should always verify the solution returned by the solver. It could be that the solver converged to an unsuitable minima, reflected by the solution configuration. By checking if there are constraint violations and the solution's cost value, you can determine the validity of a solution. By adding such a verification procedure in the code, the constraint weight tuning will be more insightful. Also, the navigation system becomes much more fault tolerant! 
#
# The verification procedure consists of the two functions defined below. The first, `ReadResults`, converts the returned solution (`Config`) into a more readable format and calculates the cost of the route. The second, `AnalyzeResult`, assesses whether the route conforms with the problem constraints. 
#
# ```python
#
# def ReadResults(Config: dict, LocationNames, CostMatrix, NumLocations):  
#
#     #############################################################################################
#     ##### Read the return result (dictionary) from the solver and sort it
#     RouteChoice = Config.items()
#     RouteChoice = [(int(k), v) for k, v in Config.items()] 
#     RouteChoice.sort(key=lambda tup: tup[0]) 
#
#     #############################################################################################
#     ##### Initialize variables to understand the routing    
#     TimeStep=[]                                      # This will contain an array of times/trips - each location is represented for each time/trip interval
#     Names = []                                       # This will contain an array of location names (string)
#     Location = []                                    # This will contain the location(s) the spaceship for each time/trip (numerical, 0 or 1)
#     RouteMatrixElements = []                         # This will contain the indices of the cost matrix representing where the spaceship has traveled (to determine total cost)
#
#     #############################################################################################
#     ##### Go through locations during each timestep/trip (all x_v) to see where the spaceship has been
#     for Index in RouteChoice:
#         TimeStep.append(math.floor(Index[0] / len(CostMatrix)))         # Time step/trip = the k-th is floor of the index divided by the number of locations
#         Names.append(LocationNames[(Index[0] % len(CostMatrix))])       # Append location names for each time step
#         Location.append(Index[1])                                       # Append location for each time step (0 or 1 if spaceship has been there or not)
#         if Index[1] == 1:                                               # Save selected location where the spaceship travels to in that trip (given by x_v = 1)
#             RouteMatrixElements.append(Index[0] % len(CostMatrix))      # Save the indices (this returns the row index)
#     SimulationResult = np.array([TimeStep, Names, Location])            # Save all the route data (also where the spaceship did not go during a turn/trip/timestep)
#  
#     #############################################################################################
#     ##### Create the route dictionary 
#     k=0                                                                                                             
#     RouteDict = {}                                                                                                                                              
#     RouteDict['Route'] = {}
#     RouteMatrix = np.array([['Timestep,', 'Location']])
#     for i in range(0, (NumLocations * (NumLocations + 1))):
#         if SimulationResult[2][i] == '1':                                                                                     # If the SimulationResult[2][i] (location) == 1, then that's where the spaceship goes/went
#             RouteMatrix = np.concatenate((RouteMatrix, np.array([[SimulationResult[j][i] for j in range(0, 2)]])), axis=0)    # append the rows where the spaceship DOES travel to to RouteMatrix
#             RouteDict['Route'].update({k: RouteMatrix[k + 1][1]})                                                             # Save the route to a dictionary
#             k += 1                                                                                                            # Iterable keeps track for the dictionary, but also allows to check for constraint
#     AnalyzeResult(RouteMatrix, NumLocations)                                                                                  # Check if RouteMatrix satisfies other constraints as well
#
#     #############################################################################################
#     ###### Calculate the total cost of the route the spaceship made (in millions of of km - distance)
#     TotalRouteCost = 0
#     for trips in range(0, NumLocations):
#         TotalRouteCost = TotalRouteCost+float(CostMatrix.item(RouteMatrixElements[trips], RouteMatrixElements[trips + 1]))     # The sum of the matrix elements where the spaceship has been (determined through the indices)
#     RouteDict['RouteCost'] = {'Cost':TotalRouteCost}
#
#     ##### Return the simulation result in a human understandable way =)
#     return RouteDict
#
#
# ############################################################################################
# ##### Check whether the solution satisfies the constraints 
# def AnalyzeResult(RouteMatrix, NumLocations):
#
#     ############################################################################################                        
#     ##### Check if the number of travels is equal to the number of locations (+ 1 for returning home)
#     if (len(RouteMatrix) - 1) != NumLocations + 1:
#         raise RuntimeError('This solution is not valid -- Number of locations visited invalid!')
#     else:
#         NumLocationsPassed = NumLocations 
#         print(f"Number of planets, moons, and asteroids passed = {NumLocationsPassed}. This is correct!")
#
#     ############################################################################################                        
#     ##### Check if the locations are different (except start/end location)
#     PastLocations = []
#     for k in range(1, len(RouteMatrix) - 1):                                                                           # Start to second last location must all be different - skip header so start at 1, skip last location so - 1
#         for l in range(0, len(PastLocations)):  
#             if RouteMatrix[k][1] == PastLocations[l]:
#                 raise RuntimeError('This route is not valid -- Traveled to a non-starting planet/moon/or asteroid more than once')
#         PastLocations.append(RouteMatrix[k][1])
#     print(f"Number of different planets, moons, asteroids passed = {NumLocations}. This is correct!") 
#
#     ############################################################################################                        
#     ##### Check if the end location is same as the start location
#     if RouteMatrix[1][1] != RouteMatrix[-1][1]:
#         raise RuntimeError(f'This route is not valid -- Start ({RouteMatrix[1][1]}) and end ({RouteMatrix[-1][1]}) location are not the same!')
#     print('Start and end location are the same. This is correct!')
#
#     ############################################################################################                        
#     ##### Check if start location and end location are correct
#     if RouteMatrix[1][1] != 'Mars' or RouteMatrix[-1][1] != 'Mars':
#         raise RuntimeError(f'This solution is not correct -- Start location {Path[1][1]} or end location {Path[-1][1]} incorrect')
#     print('Start and end location are as planned. This is correct!')
#
#     print('Valid route!')
#
# ```
#
# ## Next steps
#
# In the next unit, you will add the code pieces together and submit the problem to the Azure solvers. By tuning the constraint weights you will improve the routes for the spaceship. By the end, you should have a working navigation system ready for the crew!
#
#
#

# # Section 8
#
# The final task to complete the navigation system is to tune the constraint weights. To do this, you will need to submit the problem to Azure Quantum workspace, and evaluate the results returned by the solver. In the previous unit, verification steps were covered to check the validity of a solution configuration. For tuning, these steps will save a lot of investigatory work. Because tuning is an iterative process and dependent on many factors (solvers, constraints, cost definitions, etc.) we will only implement the **simulated annealing** solver, but feel free to uncomment other solver methods in the code and experiment!
#
#
# ## The final cost function
#
# The optimization problem, which you have programmed into the navigation system in the previous units, will be submitted to the Azure Quantum solvers to find efficient routes through the solar system! Before that, you will need to assign a value to the weights $w_1$, $w_2$, $w_3$, $w_4$, $w_5$. 
#
# $$ \text{Cost function so far: } \underset{x_0, x_1,\dots,x_{(N^2+2N)}}{min} \left(\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)  +   w_1 \left( \sum_{l=0}^{N} \sum_{i=0}^{N-1} \sum_{j=0}^{N-1} x_{(i+Nl)} \cdot x_{(j+Nl)} \text{ with } \{ i,j | i<j \} \right) +w_2\left( \sum_{k=0}^{N(N+1)-1} x_k \right) + w_3 \left(  \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)   \right)  + w_4 ( x_0)  + w_5 ( x_{N^2}) \right)  $$
#
# Remember that $w_1$ and $w_3$ should be a positive value (a cost) and that $w_2$, $w_4$, and $w_5$ should be a negative value (negative cost, reward).
#
# ### Initial weighting
#
# The travel costs of the spaceship can change when the space rocks need to be mined or sold on different planets. For the cost function it is important to take this into account, as you want the travel costs and constraints to have a constant relative weighting. Otherwise it could be the case that the solver punishes a constraint much less strictly, resulting in poor routes. Weighting the constraints as a function of the cost matrix as well, like the travel costs, will allow you to have constant relative weighting in the cost function. Then regardless of the distances between planets, moons, and asteroids considered, the navigation system should generate efficient routes. 
#
# There are a number of ways to implement a relative weighting, which boils down to a design decision. Here, the constraint weights are tuned by multiplying some factor $\alpha$ by the maximum element in the cost matrix $C$:
#
# $w_h = \alpha \cdot \max(C_{i,j}).$
#
# Other functions, for example matrix norms, are good candidates as well. You will need to re-tune the constraints if you decide to change them!
#
#
# - Some tips for the tuning process:
#   - Start by assigning the start and end locations $w_4$ and $w_5$ with a big negative weight (but keep $\alpha \leq$  10). Once this constraint is always satisfied, work on the other constraint weights.
#   - Look at the printed route, and see which constraint it violates. Then try to tune the respective constraint(s) accordingly. Iteratively try to do this and you will eventually find good routes!
#
#
# ## Submitting the code
#
# In the code below, the constraints weights have all been initialized dependent on the maximum element in the cost matrix. It is your job to experiment with around with the weights given in the constraint ($w_1$, $w_2$, $w_3$, $w_4$, $w_5$). If you want to skip the tuning process, you can use the values below:
#
# - $w_1: int(2 \cdot np.max(CostMatrix))$
# - $w_2: int(-1.65 \cdot np.max(CostMatrix))$
# - $w_3: int(2 \cdot np.max(CostMatrix))$
# - $w_4: int(-10 \cdot np.max(CostMatrix))$
# - $w_5: int(-10 \cdot np.max(CostMatrix))$
#
# The code is submitted to the Azure solvers as a `synchronous` job because you will be manually tuning and evaluating the results sequentially. Another option is submit the jobs `asynchronously`, which is helpful if you want to solve many different problem instances and compare them afterwards. Asynchronous submissions are especially helpful if you want to automate parameter tuning through grid searches, for example.
#
# >[!Note]
# >It is good practice to use `int`s instead of `float`s for the cost function weights. Doing so will give more accurate results. It is not mandatory, however.
#
# ``` python
#
# ### FILL IN THE CONSTRAINT WEIGHTS w_1, w_2, w_3, w_4, AND w_5 IN THE PROBLEM DEFINITION TO SUBMIT!
#
# from azure.quantum import Workspace
# from azure.quantum.optimization import Problem, ProblemType, Term, HardwarePlatform, Solver
# from azure.quantum.optimization import SimulatedAnnealing, ParallelTempering, Tabu, QuantumMonteCarlo
# from typing import List
#
# import numpy as np
# import math
#
# workspace = Workspace (
#     subscription_id = "",  # Add your subscription id
#     resource_group = "",   # Add your resource group
#     name = "",             # Add your workspace name
#     location = ""          # Add your workspace location (for example, "westus")
# )
#
# workspace.login()
#
# ##### Define variables
#
# # The number of planets/moons/asteroids.
# NumLocations = 10
#
# # Location names. Names of some of the solar system's planets/moons/asteroids.
# LocationNames = {0:'Mars', 1:'Earth', 2:"Earth's Moon", 3:'Venus', 4:'Mercury', 5:'Ganymede', 6:'Titan', 7:'Ceres', 8:'Pallas', 9:'Cybele'}
#
# # Approximate mean distances between the planets/moons/asteroids. Note that they can be very innacurate as orbital mechanics are ignored. 
# # This is a symmetric matrix since we assume distance between planets is constant for this module.
# CostMatrix = np.array([     [0,   78,       2,  120,   170,  550,  1200,  184, 600,  1.5   ],
#                             [78,   0,     0.5,   41,    92,  640,  1222,  264, 690,  0.25  ],
#                             [2,   0.5,      0,   40,    91,  639,  1221,  263, 689,  0.25  ], 
#                             [120,  41,     40,    0,    50,  670,  1320,  300, 730,  41.5  ],
#                             [170,  92,     91,   50,     0,  720,  1420,  400, 830,  141.5 ],
#                             [550,  640,   639,  670,   720,    0,   650,  363,  50,  548   ],
#                             [1200, 1222, 1221, 1320,  1420,  650,     0, 1014,  25,  625   ],  
#                             [184,  264,   263,  300,   400,  363,  1014,    0, 100,  400   ],
#                             [600,  690,   689,  730,   830,   50,    25, 100,    0,  350   ],
#                             [1.5,  0.25, 0.25, 41.5, 141.5,  548,   625, 400,  350,  0     ]
#                       ])    
#                        
# ##### If you want try running with a random cost matrix, uncomment the following:
# #maxCost = 10
# #CostMatrix = np.random.randint(maxCost, size=(NumLocations,NumLocations))
#  
# ############################################################################################
# ##### Define the optimization problem for the Quantum Inspired Solver
# def OptProblem(CostMatrix) -> Problem:
#     
#     #'terms' will contain the weighting terms for the solver!
#     terms = []
#
#     ############################################################################################
#     ##### Cost of traveling between locations  
#     for k in range(0,len(CostMatrix)):                          # For each trip (there are N trips to pass through all the locations and return to home base)
#         for i in range(0,len(CostMatrix)):                      # For each origin (reference location)
#             for j in range(0,len(CostMatrix)):                  # For each destination (next location w.r.t reference location)
#                 
#                 #Assign a weight to every possible trip from location i to location j - for any combination
#                 terms.append(
#                     Term(
#                         c = CostMatrix.item((i,j)),                                     # Element of the cost matrix
#                         indices = [i+(len(CostMatrix)*k), j+(len(CostMatrix)*(k+1))]    # +1 to denote dependence on next location
#                     )
#                 )
#                 ##----- Uncomment one of the below statements if you want to see how the weights are assigned! -------------------------------------------------------------------------------------------------
#                 #print(f'{i+(len(CostMatrix)*k)}, {j+(len(CostMatrix)*(k+1))}')                                                                 # Combinations between the origin and destination locations 
#                 #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)
#                 #print(f'For location_{i}, to location_{j} in trip number {k} costs: {CostMatrix.item((i,j))}')                                         # In a format that is easier to read for a human
#
#     ############################################################################################
#     ##### Constraint: One location at a time constraint - spaceship can only be in 1 location at a time.
#     for l in range(0,len(CostMatrix)+1):                # The total number of locations that are visited over the route (N+1 because returning to home base)
#         for i in range(0,len(CostMatrix)):              # For each location (iterate over the location vector)
#             for j in range(0,len(CostMatrix)):          # For each location (iterate over the location vector)
#                 if i!=j and i<j:                        # i<j because we don't want to penalize twice // i==j is forbidden (this could equal 1, that's why)
#                     terms.append(
#                         Term(
#                             c = int(w_1 * np.max(CostMatrix)),       # FILL IN - w_1 should be a positive value
#                             indices = [i+(len(CostMatrix)*l),j+(len(CostMatrix)*l)]                   
#                         )
#                     )
#                     ##----- Uncomment one of the below statements if you want to see how the weights are assigned! -------------------------------------------------------------------------------------------------
#                     #print(f'{i+(len(CostMatrix)*l)},{j+(len(CostMatrix)*(l))}')
#                     #print(f'Location constraint 1: x_{i+(len(CostMatrix)*l)} - x_{j+(len(CostMatrix)*(l+1))} (trip {l}) assigned weight: FILL IN')  # In a format for the solver (as formulated in the cost function)
#
#     ############################################################################################
#     ##### Constraint: No dissapearing constraint - encourage the spaceship to be 'somewhere' otherwise all x_k might be 0 (for example).
#     for v in range(0, len(CostMatrix) + len(CostMatrix) * (len(CostMatrix))):    # Select variable 
#         terms.append(
#             Term(
#                 c =  int(w_2 * np.max(CostMatrix)),          # FILL IN - w_2 should be a negative value
#                 indices = [v]   
#             )
#         )
#         ##----- Uncomment one of the below statements if you want to see how the weights are assigned! -------------------------------------------------------------------------------------------------
#         #print(v)
#         #print(f'No dissapearing constraint 2: x_{v} assigned weight: FILL IN')                                                      # In a format for the solver (as formulated in the cost function)
#         #print(f'No dissapearing constraint 2: location_{v % NumLocations} after {np.floor(v / NumLocations)} trips assigned weight: FILL IN')   # In a format that is easier to read for a human
#
#     ############################################################################################                        
#     ##### Constraint: no revisiting locations constraint --- (in the last step we can travel without penalties (this is to make it easier to specify an end location ))
#     for p in range(0,len(CostMatrix)+len(CostMatrix)*(len(CostMatrix))):                                  # This selects a present location x: 'p' for present    
#         for f in range(p+len(CostMatrix),len(CostMatrix)*(len(CostMatrix)),len(CostMatrix)):              # This selects the same location x but after upcoming trips: 'f' for future
#             terms.append(
#                 Term(
#                     c = int(w_3 * np.max(CostMatrix)),  # FILL IN - w_3 should be a positive value
#                     indices = [p,f]   
#                 )
#             )     
#             ##----- Uncomment one of the below statements if you want to see how the weights are assigned! -------------------------------------------------------------------------------------------------
#             #print(f'x_{p},x_{f}')  # Just variable numbers 
#             #print(f'Visit once constraint: x_{p} - x_{f}  assigned weight: FILL IN')  # In a format for the solver (as formulated in the cost function)
#             #print(f'Visit once constraint: location_{p%NumLocations} - location_{(p+f)%NumLocations} after {(f-p)/NumLocations} trips assigned weight: FILL IN')  # In a format that is easier to read for a human
#             
#     #############################################################################################                        
#     ##### Begin at x_0 (Mars - home base)
#     terms.append(
#         Term(
#             c = int(w_4 * np.max(CostMatrix)),  # FILL IN - w_4 should be a negative value,
#             indices = [0]   
#         )
#     )
#
#     ############################################################################################                        
#     ##### End at x_{N^2} (Mars - home base)
#     terms.append(
#         Term(
#             c = int(w_5 * np.max(CostMatrix)),  # FILL IN - w_5 should be a negative value,
#             indices = [len(CostMatrix)*(len(CostMatrix))]   
#         )    
#     )
#
#     return Problem(name="Spaceship navigation system", problem_type=ProblemType.pubo, terms=terms)
#     
# def ReadResults(Config: dict, LocationNames, CostMatrix, NumLocations):  
#
#     #############################################################################################
#     ##### Read the return result (dictionary) from the solver and sort it
#     RouteChoice = Config.items()
#     RouteChoice = [(int(k), v) for k, v in Config.items()] 
#     RouteChoice.sort(key=lambda tup: tup[0]) 
#
#     #############################################################################################
#     ##### Initialize variables to understand the routing    
#     TimeStep=[]                                      # This will contain an array of times/trips - each location is represented for each time/trip interval
#     Names = []                                       # This will contain an array of location names (string)
#     Location = []                                    # This will contain the location(s) the spaceship for each time/trip (numerical, 0 or 1)
#     RouteMatrixElements = []                         # This will contain the indices of the cost matrix representing where the spaceship has traveled (to determine total cost)
#
#     #############################################################################################
#     ##### Go through locations during each timestep/trip (all x_v) to see where the spaceship has been
#     for Index in RouteChoice:
#         TimeStep.append(math.floor(Index[0] / len(CostMatrix)))         # Time step/trip = the k-th is floor of the index divided by the number of locations
#         Names.append(LocationNames[(Index[0] % len(CostMatrix))])       # Append location names for each time step
#         Location.append(Index[1])                                       # Append location for each time step (0 or 1 if spaceship has been there or not)
#         if Index[1] == 1:                                               # Save selected location where the spaceship travels to in that trip (given by x_v = 1)
#             RouteMatrixElements.append(Index[0] % len(CostMatrix))      # Save the indices (this returns the row index)
#     SimulationResult = np.array([TimeStep, Names, Location])            # Save all the route data (also where the spaceship did not go during a turn/trip/timestep)
#  
#     #############################################################################################
#     ##### Create the route dictionary 
#     k=0                                                                                                             
#     RouteDict = {}                                                                                                                                              
#     RouteDict['Route'] = {}
#     RouteMatrix = np.array([['Timestep,', 'Location']])
#     for i in range(0, (NumLocations * (NumLocations + 1))):
#         if SimulationResult[2][i] == '1':                                                                                     # If the SimulationResult[2][i] (location) == 1, then that's where the spaceship goes/went
#             RouteMatrix = np.concatenate((RouteMatrix, np.array([[SimulationResult[j][i] for j in range(0, 2)]])), axis=0)    # append the rows where the spaceship DOES travel to to RouteMatrix
#             RouteDict['Route'].update({k: RouteMatrix[k + 1][1]})                                                             # Save the route to a dictionary
#             k += 1                                                                                                            # Iterable keeps track for the dictionary, but also allows to check for constraint
#     AnalyzeResult(RouteMatrix, NumLocations)                                                                                  # Check if RouteMatrix satisfies other constraints as well
#
#     #############################################################################################
#     ###### Calculate the total cost of the route the spaceship made (in millions of of km - distance)
#     TotalRouteCost = 0
#     for trips in range(0, NumLocations):
#         TotalRouteCost = TotalRouteCost+float(CostMatrix.item(RouteMatrixElements[trips], RouteMatrixElements[trips + 1]))     # The sum of the matrix elements where the spaceship has been (determined through the indices)
#     RouteDict['RouteCost'] = {'Cost':TotalRouteCost}
#
#     ##### Return the simulation result in a human understandable way =)
#     return RouteDict
#
# ############################################################################################
# ##### Check whether the solution satisfies the constraints 
# def AnalyzeResult(RouteMatrix, NumLocations):
#
#     print(RouteMatrix) # print the route to the terminal so that you can debug/re-weight the constraints
#
#     ############################################################################################                        
#     ##### Check if the number of travels is equal to the number of locations (+ 1 for returning home)
#     if (len(RouteMatrix) - 1) != NumLocations + 1:
#         raise RuntimeError('This solution is not valid -- Number of locations visited invalid!')
#     else:
#         NumLocationsPassed = NumLocations 
#         print(f"Number of planets, moons, and asteroids passed = {NumLocationsPassed}. This is correct!")
#
#     ############################################################################################                        
#     ##### Check if the locations are different (except start/end location)
#     PastLocations = []
#     for k in range(1, len(RouteMatrix) - 1):                                                                           # Start to second last location must all be different - skip header so start at 1, skip last location so - 1
#         for l in range(0, len(PastLocations)):  
#             if RouteMatrix[k][1] == PastLocations[l]:
#                 raise RuntimeError('This route is not valid -- Traveled to a non-starting planet/moon/or asteroid more than once')
#         PastLocations.append(RouteMatrix[k][1])
#     print(f"Number of different planets, moons, asteroids passed = {NumLocations}. This is correct!") 
#
#     ############################################################################################                        
#     ##### Check if the end location is same as the start location
#     if RouteMatrix[1][1] != RouteMatrix[-1][1]:
#         raise RuntimeError(f'This route is not valid -- Start ({RouteMatrix[1][1]}) and end ({RouteMatrix[-1][1]}) location are not the same!')
#     print('Start and end location are the same. This is correct!')
#
#     ############################################################################################                        
#     ##### Check if start location and end location are correct
#     if RouteMatrix[1][1] != 'Mars' or RouteMatrix[-1][1] != 'Mars':
#         raise RuntimeError(f'This solution is not correct -- Start location {Path[1][1]} or end location {Path[-1][1]} incorrect')
#     print('Start and end location are as planned. This is correct!')
#
#     print('Valid route!')
#
# ############################################################################################                        
# ##### Define the solver --- uncomment if you wish to use a different one --- you can also change the timeout.
#
# solver = SimulatedAnnealing(workspace, timeout = 120)   
# #solver = ParallelTempering(workspace, timeout = 120)
# #solver = Tabu(workspace, timeout = 120)
# #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, you will need to tune these solver parameters. 
#
# ############################################################################################                        
# ##### Submit the optimization problem to the Azure Quantum solvers --- uncomment if you wish to use a different one --- you can also change the timeout.
# OptimizationProblem = OptProblem(CostMatrix)
# route = solver.optimize(OptimizationProblem)     # Solve the optimization problem as a synchronous job -- wait until done.
#
# PathDict = ReadResults(route['configuration'], LocationNames, CostMatrix, NumLocations)
# print(PathDict)
# ```
#
# ## A working navigation system
#
# Congratulations, you have finished building the navigation system for the spaceship! With the routes provided by the navigation system there is no doubt that the mining expeditions will be succesful. To give an impression of how difficult this problem would have been to solve by hand, it would require writing out 362,880 routes. A lot of effort saved! You can now also try adding more planets, moon, and asteroids to the itinerary to see how the routes and difficulty of the problem change. 
#
#
#
#
#
#

# # Section 9 
#
#
#
# Choose the best response for each question. Then select **Check your answers**:
#
#
# ### Questions
#
# 1. 'Why is the recurrence relation between the origin and destination vectors important?'
#   - A) 'It is modeled as a constraint in the cost function.'
#   - B) 'The number of variables in the vectors is reduced.' 
#   - C) 'The number of vectors, or variables, that the solver needs to optimize for is reduced.'     
#             
# 2. 'What is the effect of assigning constraints, or terms, negative weights?'
#   - A) 'This decreases cost function values for a preferenced solution configuration. '
#   - B) 'To be certain the solver settles for a particular solution.'
#   - C) 'It creates a minima in which the solver can settle.'
#
# 3. 'Which constraints can be combined to merge cost function terms?'
#   - A) 'Some of the "one location at a time" constraint terms can be merged with some of the "no revisiting locations" constraint terms.'
#   - B) 'The "start and end locations" constraint terms can be merged with some of the "no dissapearing constraint" terms.'
#   - C) 'Some of the "no dissapearing" constraint terms can be combined with some of the "no revisiting locations" constraint terms.'
#
# 4. 'Which of the following statements is FALSE?'
#   - A) 'Each solver has its advantages and disadvantages, and their applicability depends on the considered optimization problem.'
#   - B) 'The tuning process is only about finding suitable cost function weights.'
#   - C) 'Verifying whether a solution is optimal for the traveling salesperson problem requires computing and comparison to all other possible routes.'
#
#
# ### Answers
#
# 1.  
#   - A) 'Incorrect. Even though you could form it into a constraint, using the relation to reduce the number of different vectors simplifies the problem significantly.'
#   - B) 'Incorrect. It is still necessary to represent each location before/after each trip. But the total number of vectors, or total number of variables, is greatly reduced.'
#   - C) 'Correct! By reducing the number of vectors in the optimization problem it becomes much easier to solve. The number of variables the solver needs to optimize for is much smaller.'  
#
# 2. 
#   - A) 'Correct! Negatively weighting constraint, or terms, locally decreases the cost function value(s). This increases the chance that the solver settles for a preferenced solution.'
#   - B) 'Incorrect. Because of the non-convexities, and how the solvers work underwater (random walks), it can not be guaranteed that a particlar solution is returned. This is why you have to verify the results!'
#   - C) 'Incorrect. This may be the case, but it is not a strict condition. For example, the maximum value of the cost function could decrease which probably will not immediataly make it a local minima.'
#
# 3.   
#   - A) 'Incorrect. There is no overlap between the terms in these two constraint. The first weights products of variables over the trips while the second weights products of variables for a location vector.'
#   - B) 'Correct! Both constraints weight individual terms. The terms for the first location and last location can be combined with a new weighting equal to the sum of the weights.'
#   - C) 'Incorrect. There is no overlap between the terms in these two constraints. The first weights individual variables, the second products of individual variables.'
#
# 4.  
#   - A) 'This statement is correct. Solvers have properties that may or may not (or less) suitable for a particular optimization problem.'
#   - B) 'This statement is incorrect. Tuning is also about finding appropriate solvers, solver parameters, and hardware.'
#   - C) 'This statement is correct. The NP-hardness is caused by the combinatorial explosion and because verification of optimality is non-polynomial time.'

# # Section 10
#
# Well done! You have completed the development of the spaceship's navigation system. 
#
# In this module, you learned:
# - What the traveling salesperson problem is, its applications, and its relevance in non-convex optimization.
# - How to model the travel costs, with real interplanetary mean-distances, to construct a cost function.
# - How to model problem constraints and incorporate these into the cost function.
# - Some brief theory on solvers and how to verify the results.
# - How to tune the constraint weights.
#
# ## Next steps
#
# To learn or challenge yourself more, you can look at increasing the difficulty of the problem:
# - Increase the number of planets, moons, and asteroids that need to be visited.
# - Try different solvers, and how to tune the solver parameters.
# - Incorporate more constraints on where the spaceship can be after some trip. 
# - Change the problem costs to include orbital mechanics, meaning the distances change over time (difficult). 
# - Try solving a small traveling salesperson problem with Q#, see the [Quantum Katas](https://github.com/Microsoft/QuantumKatas/) on how you could approach doing this.
#
# Alternatively, you may also consider variations of the traveling salesperson problem:
# - Vehicle routing problem (multiple spaceships)
# - Consider how many space rocks the spaceship can carry (constraint on when to mine or sell).
# - Consider having constraints on the spaceship's energy consumption.
#
# For more quantum inspired optimization samples, you can check out the [qio samples repo](https://github.com/microsoft/qio-samples/tree/main/samples). More reading on the wide-range of quantum related can be found in the [Microsoft Quantum documentation](/quantum/?azure-portal=true).
#
#
#

StarcoderdataJupyterScriptsDedupFiltered