In this notebook we'll cover how to formulate an optimization problem to find a path through a 3D lattice that does not cross itself (self-avoiding). The optimization problem is solved with the Azure Quantum service which contains numerous heuristic (non-linear) optimization solvers. Finding a self-avoiding path in a 3D lattice can be considered a difficult problem, in the sense that it has suffers from exponential scaling conditioned on the number of turns and the number of dimensions (2D vs 3D, etc.).
+In this notebook, a step by step approach is taken to explain the function definitions to are necessary to create the optimization problem.
+Goal: To find a path in a 3D lattice that does not cross itself.
+To clarify some vocab and assumptions used in the notebook:
+Note: Because the number of terms to for this optimization problem grows explosively with the number of turns, you might want to try running this in the Azure Quantum notebooks. The more technical reason for this is that expansions of nonlinear terms have to be computed locally resulting in massive number of terms needing to be uploaded. Unless you are very patient, a recommendation would be to either check out the online notebook experience or rewrite this notebook's problem class to its streaming counterpart (https://docs.microsoft.com/en-us/azure/quantum/optimization-streaming-problem).
+ +import time
+from math import floor, log2
+from azure.quantum import Workspace
+from azure.quantum.optimization import Problem, ProblemType, Term, ParallelTempering, Tabu, SimulatedAnnealing
+from azure.identity import ClientSecretCredential
+from mpl_toolkits import mplot3d
+import numpy as np
+import matplotlib.pyplot as plt
+To run the sample, you'll need to have a quantum workspace. Check out this module if you don't have one yet: https://docs.microsoft.com/en-us/learn/modules/get-started-azure-quantum/. +Fill in the variables below.
+ +workspace = Workspace(
+ subscription_id = "",
+ resource_group = "",
+ name = "",
+ location = "",
+ credential = ClientSecretCredential(tenant_id="",
+ client_id="",
+ client_secret="")
+)
+Throughout the notebook you'll probably realize that printing the "terms" isn't going to provide much help, mainly because there will be too many of them and their dictionary format. A term is the way a mathematical term is expressed in terms of a dictionary for the SDK (see examples). Below a function is defined that prints the term(s) dictionaries that describe the constraint/cost function in terms of "q" optimization variables. You can call the function to check if the working is correct, or if you're unsure of what a constraint looks like. Note that for a large problem it nevertheless becomes difficult to understand the entire function output because of the number of terms.
+An example: +{'c': 1, 'ids': [0, 10, 11, 12]} ---> $1q_0q_{10}q_{11}q_{12}$
+ +def print_function(terms: list):
+
+ '''
+ Purpose:
+ Takes a list of terms and prints it out as a mathematical (cost/constraint) function.
+ Example:
+ {'c': 1, 'ids': [0, 10, 11, 12]} ---> 1q_0q_10q_11q_12
+ Inputs:
+ 1. terms: the list of terms.
+ '''
+
+ k = 0
+ final_string = ''
+ final_string = ''
+ for term in terms:
+ term = term.to_dict()
+ weight = term['c']
+ ids = term['ids']
+ string = '('
+ if weight >= 0:
+ if k == 0:
+ string = '('+str(weight)
+ else:
+ string = '+' + '('+str(weight)
+ if weight < 0:
+ if k == 0:
+ string = '(' + str(abs(weight))
+ else:
+ string = '-' + '(' + str(abs(weight))
+ for id_ in ids:
+ string = string + f'q_{id_}'
+ string = string + ')'
+ final_string = final_string + string
+ k += 1
+
+ print('[ ' + final_string + ' ]')
+Throughout the notebook we will need to cross-multiply nonlinear terms, for computing squares for example. The function below performs the expansion. An important note to keep in mind here, especially if you're looking to optimize the code, is that the function does not assemble identical terms. Since the cost functions dealt with in this notebook are still rather small, reducing the number of terms was left out of the scope.
+More clearly stated in terms of equations:
+$$ \text{Current function (see below): } \hspace{0.2cm} (a+b)^2 = a^2 + ab + ab + b^2 $$$$ \text{Optimal function: } \hspace{0.2cm} (a+b)^2 = a^2 + 2ab + b^2 $$Should be a fun programming exercise to reduce the number of terms significantly!
+ +def cross_multiply(list_a: list, list_b: list) -> list:
+
+ '''
+ Purpose: Cross multiplies two lists of terms (linear ex. [q_0 + q_1] or non-linear ex. [q_0q_1]) to return the expansion.
+ Can compute powers of groups this way (^2, ^3,...), like squaring a list of terms.
+ Calculates the expansion locally, unlike the SlcTerm class.
+ Example: (2q_0q_1+3q_2q_3)^2 => (2q_0q_1)^2 + 6q_0q_1q_2q_3 + 6q_0q_1q_2q_3 + (3q_2q_3)^2
+ Input:
+ 1. list_a: list which serves as the reference list (first 'for' loop).
+ 2. list_b: list which serves as the target list (second 'for' loop).
+ Output:
+ 1. list of term objects.
+ '''
+
+ terms = []
+ for one in list_a:
+ for uno in list_b:
+ alpha = one.to_dict()
+ beta = uno.to_dict()
+ weight = int(alpha['c'])*int(beta['c'])
+ ids = list(alpha['ids'])+list(beta['ids'])
+ terms += [Term(c = weight, indices = ids)]
+ return terms
+While running through the notebook you might want to execute and print some of the constraint functions. To do that, you'll need some variables to be defined. Only 3 turns are considered with the defined parameters below, as that will keep the function outputs small and understandable!
+ +nodes = "ABCD" # The four position/node names.
+len_seq = len(nodes) # The length of the node
+num_turns = len(nodes)-1 # The number of turns, which is one less than the number of therefore -1 turns.
+num_dim = 3 # Number of optimization variables required to describe a turn
+lambda_0 = 1 # Penalty weight for the 'distance constraint'
+lambda_1 = 40 # Penalty weight for the 'no return constraint'
+lambda_2 = 25 # Penalty weight for the invalid direction '000' constraint
+lambda_3 = 25 # Penalty weight for the invalid direction '111' constraint
+Prior to starting to define our cost function, we have to construct some basic functions and agree on some definitions. In this part, we'll define the 6 directions of the x-y-z plane in terms of the optimization variables ("q").
+A common way to encode a turn based optimization problem is through a direction-based string, in which each optimization variable substring (ex. $q_0q_1q_2$) denotes some kind of decision of which direction was takn. For example, take a car take can only move backward or forward. The decision of the driver at some time $t$ can be described as $q_t \in \{0, 1\}$, where 0 represents going backward and 1 forward, respectively. Over a multiple time steps, say 3, one can then describe the movement of the car by:
+$$ \text{Movements:} \hspace{0.25cm} q_0q_1q_2 $$The sequence of movements, "forward, backward, backward", can then be represented as:
+$$ \text{Sequence of movements:} \hspace{0.25cm} q_0q_1q_2 = 100$$Alright, if that is clear to you, then the next step is to represent the 6 directions of the x-y-z plane in such a fashion. Regrettably, a "q" can only be either 0 or 1, not 2, 3, 4, 5 to represent the six directions. But what can be done is to use multiple optimization variables to describe a direction! For example:
+$$ \text{Sequence of movements:} \hspace{0.25cm} q_0q_1q_2q_3q_4q_5q_6q_7q_8 $$$$ \text{First move:} \hspace{0.25cm} q_0q_1q_2 $$$$ \text{Second move:} \hspace{0.25cm} q_3q_4q_5 $$$$ \text{Third move:} \hspace{0.25cm} q_6q_7q_8 $$Because we need to encode 6 directions, we need at least 3 optimization variables to form unique substrings. The table explains this, we to have enough possible combinations to encode the 6 directions:
+Note, for later, you should keep in mind that 2 of the 8 combinations become irrelevant. It is important that the solver does not return those strings!
+Perfect! Now we just need to assign each direction a unique substring. Below you can find my choices, feel free to change these if you're starting from scratch.
+Now we just need to define direction variables. These are as defined in the above table. These variables should take a value 1 if a movement is made in that direction, while the other directions must all be zero. This can be achieved by the following scheme, where $k$ denotes the turn, and $\gamma = 3(k-1)$ (3 because 3 dimensions):
+To clarify a bit further through an example: +If in the first turn (k=0) the solver returns $q_0=1$, $q_1=0$, $q_2=0$, and in the second turn (k=1) returns $q_3=0$, $q_4=0$, $q_5=1$, then a "+x and +y" were taken, respecively.
+These direction variables will help understand some difficult constraints later on. Defining these direction variables additionally makes the code much more readable, since we can define everything in terms of the direction variables instead of the individual optimization variables.
+Below you can find the function for the direction variables. The formulas have been expanded (mathematically) to easily declare them in a number of terms, however it is hard-coded this way. Some inputs to the function ("sign_dir", "signpos", "lambda") are relevant for building constraints later in the notebook, but have to defined in this function. Further information can be found in the function docstring.
+ +def direction_variables(direction: str, offset: int, sign_dir: int, sign_pos: int, lambda_: int) -> list:
+
+ '''
+ Purpose:
+ Translates the direction (+x,-x,+y,-y,+z,-z) of turn 'i' as a function of three optimization variables. (Three q's because of the defined coordinate system).
+ Example:
+ Direction "+z" in the first turn (turn = 1) is translated to: q_{0+offset}q_{2+offset}-q_{0+offset}q_{1+offset}q_{2+offset}.
+ Inputs:
+ 1. direction: A direction from an x-y-z coordinate system, one of the following: ('+x','-x','+y','-y','+z','-z').
+ 2. offset: Offset gives the turn number expressed in the first "q" of that turn. Equal to gamma in the explanation.
+ Example: Turn 1 starts with "q" q_0, offset=0. Turn 2 starts with q_3, offset = 3.
+ 3. sign_dir: Changes the sign of the weights corresponding to negative directions "-x", "-y", "-z" -> necessary for finding the positions, for exmaple (+x) "-" (-x).
+ 4. sign_pos: Changes the sign of the weights corresponding to negative positions "-(x,y,z)" -> necessary for finding the distances between node "i" and node "j".
+ 5. lambda_: The weight term associated with a constraint.
+ Output:
+ 1. A list of term objects.
+ '''
+
+ terms = []
+ if direction == "+x":
+ term_0 = Term(c= 1*sign_pos*lambda_, indices=[0+offset])
+ term_1 = Term(c=-1*sign_pos*lambda_, indices=[0+offset, 1+offset])
+ term_2 = Term(c=-1*sign_pos*lambda_, indices=[0+offset, 2+offset])
+ term_3 = Term(c= 1*sign_pos*lambda_, indices=[0+offset, 1+offset, 2+offset])
+ terms = [term_0, term_1, term_2, term_3]
+ elif direction == "-x":
+ term_0 = Term(c= 1*sign_dir*sign_pos*lambda_, indices=[1+offset])
+ term_1 = Term(c=-1*sign_dir*sign_pos*lambda_, indices=[0+offset, 1+offset])
+ term_2 = Term(c=-1*sign_dir*sign_pos*lambda_, indices=[1+offset, 2+offset])
+ term_3 = Term(c= 1*sign_dir*sign_pos*lambda_, indices=[0+offset, 1+offset, 2+offset])
+ terms = [term_0, term_1, term_2, term_3]
+ elif direction == "+y":
+ term_0 = Term(c= 1*sign_pos*lambda_, indices=[2+offset])
+ term_1 = Term(c=-1*sign_pos*lambda_, indices=[0+offset, 2+offset])
+ term_2 = Term(c=-1*sign_pos*lambda_, indices=[1+offset, 2+offset])
+ term_3 = Term(c= 1*sign_pos*lambda_, indices=[0+offset, 1+offset, 2+offset])
+ terms = [term_0, term_1, term_2, term_3]
+ elif direction == "-y":
+ term_0 = Term(c= 1*sign_dir*sign_pos*lambda_, indices=[0+offset, 1+offset])
+ term_1 = Term(c=-1*sign_dir*sign_pos*lambda_, indices=[0+offset, 1+offset, 2+offset])
+ terms = [term_0, term_1]
+ elif direction == "+z":
+ term_0 = Term(c= 1*sign_pos*lambda_, indices=[0+offset, 2+offset])
+ term_1 = Term(c=-1*sign_pos*lambda_, indices=[0+offset, 1+offset, 2+offset])
+ terms = [term_0, term_1]
+ elif direction == "-z":
+ term_0 = Term(c= 1*sign_dir*sign_pos*lambda_, indices=[1+offset, 2+offset])
+ term_1 = Term(c=-1*sign_dir*sign_pos*lambda_, indices=[0+offset, 1+offset, 2+offset])
+ terms = [term_0, term_1]
+ return terms
+
+
+##### ----- Test the function and print output:
+turn = 1 # play with this value!
+dir_var = direction_variables("+x",(turn-1)*num_dim,1,1,lambda_0)
+print('dir_var term dictionary "+x": ', dir_var)
+print_function(dir_var)
+Great that we have the directions defined! Now let's use them and define a function that calculates the difference in positions, which we will need later. After each turn, a direction that appends to the path. A valid path is defined as a one that does not cross itself, meaning that we need to compare the positions over the turns. In other words, every new position appended to the path needs to be checked with all previously held positions, to verify that a position hasn't been visited twice.
+So how do we go about this?
+First we need a way to know the positions after choosing to go in a certain direction, so let's tackle that first. Consider the fact that the directions are recorded for each turn, which already gives some sort of log of the positions in the path. By summing the direction variables respective to their dimension (x/y/z) acoomplishes this, but only if a negative move can compensate a positive one. For example, after two turns the $x$-location can be expressed as:
+$$ x_2 = ( d^{1}_{+x} - d^{1}_{-x} ) + ( d^{2}_{+x} - d^{2}_{-x} )$$NOTE: In the above function ("direction_variables"), the negative sign for the negative direction variables is controlled through "sign_dir". "sign_dir" must be set to -1 to assign a negative sign to it as required in these position formulas!
+Remember that for a turn only one direction variable can be set to 1! Thus if two moves in the $-x$ direction are taken, then:
+$$ x_2 = ( d^{1}_{+x} - d^{1}_{-x} ) + ( d^{2}_{+x} - d^{2}_{-x} ) = ( 0 - 1 ) + ( 0 - 1 ) = -2 $$As with the simple example for the x-dimension, the same can be applied to the $x$ and $y$ directions. The formulas below descrive the position after some number of turns through a summation (for the maths/physics enthousiasts, these are integrals of the velocities to derive the positions :) ).
+$$ x^k = \sum_{k=1}^{k} ( d^{k}_{+x} - d^{k}_{-x} ) $$$$ y^k = \sum_{k=1}^{k} ( d^{k}_{+y} - d^{k}_{-y} ) $$$$ z^k = \sum_{k=1}^{k} ( d^{k}_{+z} - d^{k}_{-z} ) $$Great stuff! Now we have a method to find the positions along the path. To find the differences in position, for example between turn 2 and turn 1, it is only necessary to apply a subtraction:
+$$ x_2 - x_1 = \sum_{k=1}^{2} ( d^{k}_{+x} - d^{k}_{-x} ) - \sum_{k=1}^{1} ( d^{k}_{+x} - d^{k}_{-x}) = ( d^{2}_{+x} - d^{2}_{-x} ) $$$$ y_2 - y_1 = \sum_{k=1}^{2} ( d^{k}_{+y} - d^{k}_{-y} ) - \sum_{k=1}^{1} ( d^{k}_{+y} - d^{k}_{-y}) = ( d^{2}_{+y} - d^{2}_{-y} ) $$$$ z_2 - z_1 = \sum_{k=1}^{2} ( d^{k}_{+z} - d^{k}_{-z} ) - \sum_{k=1}^{1} ( d^{k}_{+z} - d^{k}_{-z}) = ( d^{2}_{+z} - d^{2}_{-z} ) $$Generalizing this to any difference between two positions:
+$$ x_j - x_i = \sum_{k=i}^{j} ( d^{k}_{+x} - d^{k}_{-x} ) $$$$ y_j - y_i = \sum_{k=i}^{j} ( d^{k}_{+y} - d^{k}_{-y} ) $$$$ z_j - z_i = \sum_{k=i}^{j} ( d^{k}_{+z} - d^{k}_{-z} ) $$Fantastic! Now on to the function code for the difference in positions. Do read the function and its docstring to understand the implementation, because it reference the previously defined function to make use of direction variables.
+ +def diff_in_pos(start_turn: int, end_turn: int, num_dim: int, lamda_: int):
+
+ '''
+ Purpose:
+ Expresses the difference in position (x,y,z) between two turns as a function of the q's encoding of the directions.
+ In other words, expresses the difference in the position after the start_turn and position after the end_turn.
+ Example:
+ Difference between turn 0 (no turns yet, initial position = (0,0,0)) and turn 2:
+ x(2) = [ move(+x, turn 1) - move(-x, turn 1) ] + [ move(+x, turn 2) - move(-x, turn 2) ]
+ Note: only one 'move' per turn gets activated as they are represented by the same q's (ex. turn 1 is represented by q_0, q_1, and q_2).
+ Inputs:
+ 1. start_turn: the initial (reference) turn.
+ 2. end_turn: the final (target) turn.
+ 3. num_dim: the number of dimensions (3, x-y-z coordinate system).
+ Outputs:
+ 1. The difference in the x direction.
+ 2. The difference in the y direction.
+ 3. The difference in the z direction.
+ '''
+
+ x_diff = y_diff = z_diff = []
+ if start_turn < end_turn and start_turn >= 0 and end_turn >= 1:
+ for turn in range(start_turn,end_turn+1):
+ x_diff += direction_variables("+x",(turn-1)*num_dim,1,1,lamda_)+direction_variables("-x",(turn-1)*num_dim,-1,1,lamda_)
+ y_diff += direction_variables("+y",(turn-1)*num_dim,1,1,lamda_)+direction_variables("-y",(turn-1)*num_dim,-1,1,lamda_)
+ z_diff += direction_variables("+z",(turn-1)*num_dim,1,1,lamda_)+direction_variables("-z",(turn-1)*num_dim,-1,1,lamda_)
+ return x_diff, y_diff, z_diff
+
+
+##### ----- Test the function and print output:
+start_turn = 1 # play with this value!
+end_turn = 2 # play with this value!
+x_diff, y_diff, z_diff = diff_in_pos(start_turn, end_turn, num_dim, lambda_0)
+print('Difference in x described as term list: \n', x_diff)
+print('\nDifference in x described as mathematical function:')
+print_function(x_diff)
+Fantastic that you've made it this far! With all the above functions defined we can start solving the problem.
+The distance constraint enforces that after each turn, the difference in positions after turns $i$ and $j$ must be larger or equal to 1 (distance between lattice points in the same dimension). +To visualize the scenario this constraint tries to penalize consider a 2D example with 4 turns, starting in $(0,0)$:
+$$ \text{Turn 1, go right: } \hspace{0.3cm} (0,0) => (0,1) $$+$$ \text{Turn 2, go up: } \hspace{0.3cm} (0,1) => (1,1) $$ +$$ \text{Turn 3, go left: } \hspace{0.3cm} (1,1) => (1,0) $$ +$$ \text{Turn 4, go down: } \hspace{0.3cm} (1,0) => (0,0) $$
+As you can see, the last step conflicts with our constraint. It is not permitted to return into a position that we've already been, namely $(0,0)$. Turn 1-3 are all valid, since they remain a distance of 1 away from all other previously held positions.
+Let's define what we want mathematically. tThe distance between the position after turn $i$ and $j$ must be larger or equal to 1. For this we'll need some basic geometry, the 3D variant of the Pythagorean theorem. From the theorem, it can be understood that hypotenuse ($c$) must always be larger or equal to one, since that is what defines the distance between two points. Mathematically speaking working this out:
+$$ C_{i,j} \geq 1 $$$$ \sqrt{(x_{j} - x_{i})^2 + (y_{j} - y_{i})^2 + (z_{j} - z_{i})^2} \geq 1 $$The square-root is a bit of a headache when writing out a cost function, it would require an approximation or some sort of factorization, not fun if you're dealing with many terms. To avoid that, we're simply going to square both sides, which in our case if fine (for maths people), since we aren't dealing with negative numbers here:
+$$ C_{i,j}^2 \geq 1 $$$$ {(x_{j} - x_{i})^2 + (y_{j} - y_{i})^2 + (z_{j} - z_{i})^2} \geq 1 $$Nice! If you're sharp, you'll start seeing where the previous functions might come in handy now. The function that computes the difference between the respective positions ("diff_in_pos") in each dimension was defined above, but also the function that computes the an expansion ("cross_multiply") was defined at the beginning of the notebook. Writing those function calls out in order :
+Alright, you might be wondering how to deal with the $>=$, since it is necessary to have an equation that equals zero in order to integrate it into the cost function. +This is where slack variables come in. The purpose of slack variables is to convert an inequality constraint into an equality constraint. Or in other words, for this equation, we need to compensate some values (explained further below) in order to convert the equation to an equality constraint. Consider the constraint in its simplest form:
+$$ C_{i,j}^2 \geq 1 $$It is clear that we would need to add some number $V_{i,j}$, on the right hand side of the equation in order to make this an equality constraint. $V_{i,j}$ can be any positive value and is conditioned on $i$ and $j$:
+$$ C_{i,j}^2 = 1+V_{i,j}$$Let's say that instead of representing a single value we want V_{i,j} to represent a range of values $[0-4]$. This is where slack variables can be introduced. By introducing additional optimization variables, which have nothing to do with the path encoding, the range can be described as following:
+$$ V_{i,j} = [0-4] = q_{s1}+q_{s2}+q_{s3}+q_{s4} $$Or more compactly:
+$$ V_{i,j} = [0-4] = q_{s1}+q_{s2}+2q_{s3} $$So any value in the range $[0-4]$ is described through these q's. You can derive 1 by assigning a single "q", which has weight 1, the value 1. Likewise, you can get the value 4 by assigning all optimization variables the value 1.
+Alright, so how do we know which range of values we need to represent? That is revealed by the turn numbers, $i$ and $j$. Over these turns, the maximum squared distance that can be achieved is when moves are only made in a single direction, the "+x" direction for example. The maximum distance is then equal to the number of moves squared $(j-i)^2$.
+$$ C_{i,j}^2 = (j-i)^2 $$Then writing out the maths gives us the necessary upper bound for the range that the optimization variables need to represent, where $S$ stands for the number of slack variables necessary:
+$$ \text{Upper bound: } \hspace{0.5cm} C_{i,j}^2 = (j-i)^2 = 1 + \sum_{s=0}^{S}q_s$$$$ \text{Upper bound: } \hspace{0.5cm} \sum_{s=0}^{S}q_s = (j-i)^2 - 1 $$For the lower bound of the range, the value is 0. The reason for that is we want to be able to represent any number equal to or larger than 1 for the inequality constraint through the addition of the slack variables. Concluding, the range the q's must represent is $[0,((j-i)^2-1)]$. More on how to calculate the number of slack variables necessary to define this range later, as you don't need ((j-i)^2-1) optimization variables!
+Piecing all the parts together we find the following equation for a $i$-$j$ combination:
+$$ {(x_{j} - x_{i})^2 + (y_{j} - y_{i})^2 + (z_{j} - z_{i})^2} = 1 + \sum_{s=0}^{S}q_s $$ +Because constants are irrelevant for the optimization landscape the value 1 can be neglected. The reason being that constants only introduce linear offsets, thus impacting the entire optimization equally. The distance constraint which needs to account for all $i$-$j$ combinations, where $j$>$i$ and $N-1$ the number of turns, is summarized as follows:
+$$ \sum_{i=1}^{N-1} \sum_{j>i}^{N-1} \lambda_0 ( {(x_{j} - x_{i})^2 + (y_{j} - y_{i})^2 + (z_{j} - z_{i})^2} - \sum_{s=0}^{S}q_s ) $$ +The function that builds this constraint is presented below. Read the docstring for how it works.
+Note, to run this function you'll need to run the next function cell to define 'generate_slack_coefficients'.
def distance_constraint(num_turns: int, num_dim: int, lambda_0: int) -> list:
+
+ '''
+ Purpose:
+ Build the distance contraint based on previosly defined functions.
+ Constraint: Distance squared between i and j must be larger or equal to 1.
+ Constraint: L_{i,j}^2 >= 1 => L_{i,j}^2 = 1+q_{slacker} (converting the inequality constraint to an equality constraint.)
+ Example/Explanation:
+ L{1,2}^2 = (x1-x2)^2 + (y1-y2)^2 + (z1-z2)^2
+ x1 = [move(+x1)-move(-x1)]
+ x2 = [move(+x1)-move(-x1)] + [move(+x2)-move(-x2)]
+ <same for other dimensions >
+ <fill into first line>
+ L{1,2}^2 = [move(+x2)-move(-x2)]^2 + [move(+y2)-move(-y2)]^2 + [move(+z2)-move(-z2)]^2
+ [move(+x2)-move(-x2)]^2 + [move(+y2)-move(-y2)]^2 + [move(+z2)-move(-z2)]^2 - q_{slacker} = 0 ---> expressed in q's == cost function
+ Inputs:
+ 1. num_turns: the number of turns.
+ 2. num_dum: the number of dimensions.
+ 3. lambda_0: the constraint weight for the distance constraint.
+ Output:
+ 1. List of term objects describing the distance constraint.
+ '''
+
+ terms = []
+ slack_indexer = 0
+ for start_turn in range(1,num_turns+1):
+ for end_turn in range(start_turn+1,num_turns+2):
+ # Calculate the differences in positions for each dimension.
+ x_diff_i_j, y_diff_i_j, z_diff_i_j = diff_in_pos(start_turn, end_turn, num_dim, lambda_0)
+ # Compute the squared distance (Pythagorean theorem) by calculating the squared expansion.
+ x_diff_i_j_2, y_diff_i_j_2, z_diff_i_j_2 = cross_multiply(x_diff_i_j,x_diff_i_j), cross_multiply(y_diff_i_j,y_diff_i_j), cross_multiply(z_diff_i_j,z_diff_i_j)
+ # Add slack variables due to inequality constraint.
+ slack_var_terms = []
+ slack_coefficients = generate_slack_coefficients(end_turn-start_turn)
+ for s in range(0,len(slack_coefficients)):
+ slack_var_terms += [Term(c=-slack_coefficients[s], indices=[num_turns*num_dim+slack_indexer+s])]
+ terms += x_diff_i_j_2 + y_diff_i_j_2 + z_diff_i_j_2 + slack_var_terms
+ slack_indexer+=len(slack_coefficients)
+ return terms
+
+##### ----- Test the function and print output:
+# The output is too large to understand, nevertheless you can run the below statements to view how visualize the term-scaling of the problem.
+number_turns = 2 #play with this value!
+dist_terms = distance_constraint(number_turns, num_dim, lambda_0)
+
+# Only print below statements if you want to get an impression of the number of terms.
+#print('Distance constraint term dictionaries: ', dist_terms) # not readable!
+#print('\nDistance constraint: ')
+#print_function(dist_terms) # not readable!
+To remove any human calculations to construct the distance constraint(s), we'll automate the computation of the number of slack variables and their weights.
+As you saw in the previous part, certain slack variables can be assigned weights to reduce the number of slack variables necessary. Recall that the range $[0,4]$ can be defined by the following two equations with slack variables (s):
+$$ V_{i,j} = [0-4] = q_{s1}+q_{s2}+q_{s3}+q_{s4} $$$$ V_{i,j} = [0-4] = q_{s1}+q_{s2}+2q_{s3} $$The bottom equation is more compact because the last variables carries a weight-2. In Lucas' paper (see 2.4 of arxiv:1302.5843) the algorithm to find such weights is presented, which is based on calculating all $2^{x}$, where $x$ iterates from zero to the log_2 of maximum value in the range (4 in this example).
+Below are two examples that may help visualize why and how the slack variables represent the ranges of values. +\ +
+\ +
+\ +
+The function is given below. In context of the distance constraint, we're still faced with a minor problem. If $(j-i)^2-1$ equals 0, then we don't need any slack variables. The reason for this is that the range becomes $[0,0]$, which can also be represented by no variables at all.
+ +def generate_slack_coefficients(turn_diff: int):
+
+ '''
+ Purpose:
+ Calculates the number of slack variables and their weights.
+ Example:
+ For the constraint: x1 + x2 + x3 + x4 <= 4 which is converted to x1 + x2 + x3 + x4 + s1 + s2 +2*s3 = 4,
+ with 's' being slack variables. This function computes the weights of these slack variables ([1,1,2] for the example).
+ Reference:
+ Lucas' paper (see 2.4 of arxiv:1302.5843)
+ Input:
+ 1. turn_diff: the differences in turns (end_turn - start_turn).
+ Output:
+ 1. y: the weights of the slack variables.
+ '''
+
+ dist_diff = (turn_diff**2)-1
+ if dist_diff == 0:
+ y = [] # no slack variables needed
+ elif dist_diff > 0:
+ M = floor(log2(dist_diff))
+ y = [-2**n for n in range(M)]
+ y.append(-(dist_diff + 1 - 2**M))
+ return y
+
+##### ----- Test the function and print output:
+turn_difference = 3 # play with this value!
+slack_weights = generate_slack_coefficients(turn_difference)
+print('Slack weights: ', slack_weights)
+The distance constraint enforces a minimum distance of 1 between the different taken positions. However, for two consecutive turns, it does not completely prevent (to a high enough degree) going to a previous position. A simplified explanation for this is that the distance constraint is described through the positional differences in $x$, $y$, and $z$. Additionally, for consecutive turns, the number of slack variables is zero ($(i-j)^2-1 = 0$ for consecutive turns} - thus there is no auxiliary benefit introduced by slack variables to compensate for being in a different position after the turn. Therefore, the solver will try to keep the positions as close to each other as possible, which leads to overlapping positions, especially in consecutive turns. To overcome this, a penalization needs to be introduced that stops returning to the previous position held in the turn before. For example, what we want to stop is the following:
+$$ \text{Turn 1: A to B.} $$$$ \text{Turn 2: B to A.} $$This can be equivalenty described as moving in opposite directions over the consecutive turns. For example, first going in the "+y" direction, and following that moving in the "-y" direction.
+Designing such a constraint is straightforward, the opposite direction variables need to be multiplied. +Preventing a revisit of a position for the "x" direction over the first and second turn can be described by the following constraint:
+$$ {\lambda}_2(d^{1}_{+x} d^{2}_{-x} + d^{1}_{-x}d^{2}_{+x}) $$ +If a move is made in the "+x" direction first, $d^{1}_{+x}$ takes value 1. If afterwards a move is made in the "-x" direction, $d^{2}_{-x}$ also takes value 1, meaning that the combined term +$d^{1}_{+x} d^{2}_{-x}$ also becomes 1, enforcing the constraint with penalty value ${\lambda}_2$.
+Generalizing the idea to any direction and for all turns, we find the "no return constraint":
+$$ {\lambda}_1 \left( \sum_{m}^{\in \{x,y,z\} }\sum_{t=1}^{N-1} d^{t}_{+m} d^{t+1}_{-m} + d^{t}_{-m}d^{t+1}_{+m} \right) $$ +In this constraint, iterations are performed over the different dimensions ($m$), and the turn ($t$). The function defintion for the no return constraint is given below. Make sure to read the docstring. In the function the 'cross_multiply' method is used to expand and calculate all the necessary terms of the constraint. The direction variables are dictionaries of polynomials, and therefore have to be expanded before submitting to the Azure QIO solvers.
+ +def no_return_constraint(num_turns: int, num_dim: int, lamda_4: int)-> list:
+
+ '''
+ Purpose:
+ Build the constraint that penalizes going back to the same position/node two turns later.
+ Example:
+ Node A ---> <move +x> ---> Node B ---> <move -x> ---> Node A => erroneous as we've been there already.
+ Two sequential moves may not be in the same dimension and in opposite directions: (+x then -x), (-x then +x), (+y then -y) etc.
+ Inputs:
+ 1. num_turns: the number of turns.
+ 2. num_dim : the number of dimensions, which is 3 for this sample.
+ 3. lambda_4 : the penalty weight for this constraint.
+ Outputs:
+ 1. List of term objects.
+ '''
+
+ terms = []
+ for i in range(0,num_turns):
+ x_out_in = cross_multiply(direction_variables("+x",i*num_dim,1,1,lamda_4), direction_variables("-x",(i+1)*num_dim,1,1,lamda_4))
+ x_in_out = cross_multiply(direction_variables("-x",i*num_dim,1,1,lamda_4), direction_variables("+x",(i+1)*num_dim,1,1,lamda_4))
+ y_right_left = cross_multiply(direction_variables("+y",i*num_dim,1,1,lamda_4), direction_variables("-y",(i+1)*num_dim,1,1,lamda_4))
+ y_left_right = cross_multiply(direction_variables("-y",i*num_dim,1,1,lamda_4), direction_variables("+y",(i+1)*num_dim,1,1,lamda_4))
+ z_up_down = cross_multiply(direction_variables("+z",i*num_dim,1,1,lamda_4), direction_variables("-z",(i+1)*num_dim,1,1,lamda_4))
+ z_down_up = cross_multiply(direction_variables("-z",i*num_dim,1,1,lamda_4), direction_variables("+z",(i+1)*num_dim,1,1,lamda_4))
+ terms += x_out_in + x_in_out + y_right_left + y_left_right + z_up_down + z_down_up
+ return terms
+
+
+##### ----- Test the function and print output:
+number_turns = 1 # play with this value!
+nrc_terms = no_return_constraint(number_turns, num_dim, lambda_1)
+print('No return constraint term dictionaries: ', nrc_terms)
+print('\n No return constraint function:')
+print_function(nrc_terms)
+Recall that there are two invalid optimization variable substrings that are not associated with any move, $000$ and $111$. A path that includes these substrings is invalid. To prevent the solver from generating these invalid directions we'll need to add two constraints. Each constraint is specific to an invalid direction, and because we're working directly with the substring we don't need the abstraction layers used in previous constraints, like directional and positional variables.
+Let's first work out the constraint for preventing substring $000$. $ q_{0+\gamma} q_{1+\gamma} q_{2+\gamma} $ may not equal $000$, where $\gamma = 3(k-1)$ and $k$ the turn number. If we were to design the constraints as $ q_{0+\gamma} q_{1+\gamma} q_{2+\gamma} $ there won't be any penalty, because all $q$'s will equal and multiply to zero. Therefore, if all these $q$'s equal zero, we need the substring to multiply to 1. This is achieved by the following, which multiplies to zero if all optimization variables take the value zero:
+$$ {\lambda}_2(1-q_{0+\gamma}) (1-q_{1+\gamma}) (1-q_{2+\gamma}) $$$$ \text{If } \hspace{0.1cm} q_{0+\gamma} = q_{1+\gamma} = q_{2+\gamma} = 0, \hspace{0.1cm} \text{then: } $$+$$ {\lambda}_2(1-q_{0+\gamma}) (1-q_{1+\gamma}) (1-q_{2+\gamma}) = {\lambda}_2 \cdot 1 $$
+Nice! By expanding the equation and neglecting constant terms a constraint is found that can be implemented:
+$$\text{Constraint for 000: }-q_{0+\gamma}−q_{1+\gamma}−q_{2+\gamma}+q_{0+\gamma}q_{1+\gamma}+q_{0+\gamma}q_{2+\gamma}+q_{1+\gamma}q_{2+\gamma}−q_{0+\gamma}q_{1+\gamma}q_{2+\gamma}$$Before going to the function definition, let's first look at penalizing the substring $111$. Luckily this is much easier as the substring multiplies to 1 if all the $q$'s have the value zero. The constraint is therefore simple to derive:
+$$ {\lambda}_3 (q_{0+\gamma} q_{1+\gamma} q_{2+\gamma}) $$$$ \text{If } \hspace{0.1cm} q_{0+\gamma} = q_{1+\gamma} = q_{2+\gamma} = 1, \hspace{0.1cm} \text{then: } $$+$$ {\lambda}_3 (q_{0+\gamma} q_{1+\gamma} q_{2+\gamma}) = {\lambda}_3 \cdot 1$$
+The constraint for the substring $111$ is:
+$$ \text{Constraint for $111$: }\hspace{0.1cm} {\lambda}_3 (q_{0+\gamma} q_{1+\gamma} q_{2+\gamma}) $$The function definitions for these two constraints are given below.
+ +def penalize_000(len_seq: int, num_dim: int, lambda_2: int) -> list:
+
+ '''
+ Purpose:
+ Build the constraint that penalizes the invalid moves associated with the 3 q's string: '000'.
+ The constraint needs to be defined over the set of all turns, as defined by the 'for' loop.
+ Example:
+ If turn 2 (q_3q_4q_5) equals '000', assign a large penalty.
+ Inputs:
+ 1. len_seq: the number nodes to consider.
+ 2. num_dim: the number of dimensions (which is 3).
+ 3. lambda_2: the penalty weight for this constraint.
+ Outputs:
+ 1. List of term objects.
+ '''
+
+ terms = []
+ for k in range(0,len_seq):
+ offset = k*num_dim
+ term_0 = Term(c=-1*lambda_2,indices=[0+offset])
+ term_1 = Term(c=-1*lambda_2,indices=[1+offset])
+ term_2 = Term(c=-1*lambda_2,indices=[2+offset])
+ term_3 = Term(c= 1*lambda_2,indices=[0+offset,1+offset])
+ term_4 = Term(c= 1*lambda_2,indices=[0+offset,2+offset])
+ term_5 = Term(c= 1*lambda_2,indices=[1+offset,2+offset])
+ term_6 = Term(c=-1*lambda_2,indices=[0+offset,1+offset,2+offset])
+ terms += [term_0, term_1, term_2, term_3, term_4, term_5, term_6]
+ return terms
+
+def penalize_111(len_seq: int, num_dim: int, lambda_3: int) -> list:
+
+ '''
+ Purpose:
+ Build the constraint that penalizes the invalid moves associated with the 3 q's string: '111'.
+ The constraint needs to be defined over the set of all turns, as defined by the 'for' loop.
+ Example:
+ If turn 2 (q_3q_4q_5) equals '111', assign a large penalty.
+ Inputs:
+ 1. len_seq: the number nodes to consider.
+ 2. num_dim: the number of dimensions (which is 3 in this sample).
+ 3. lambda_3: the penalty weight for this constraint.
+ Outputs:
+ 1. List of term objects.
+ '''
+
+ terms = []
+ for k in range(0,len_seq):
+ offset = k*num_dim
+ terms += [Term(c=1*lambda_3,indices=[0+offset,1+offset,2+offset])]
+ return terms
+
+
+##### ----- Test the function and print output:
+term_000 = penalize_000(len_seq, num_dim, lambda_2)
+term_111 = penalize_111(len_seq, num_dim, lambda_3)
+print('Term dictionary penalty constraint for 000:', term_000, '\n')
+print('Term dictionary penalty constraint for 111:', term_111, '\n')
+print('Function for penality 000: ')
+print_function(term_000)
+print('Function for penality 111: ')
+print_function(term_111)
+We've finished defining the optimization function. The next step is to start looking at how we're going to submit it to the solvers and analyze the returned results. +First let's look at how to parse, validate, and visualize the results, since that will make tuning the solvers easier!
+Below you can find the function definition that reads the solution and validates it. The solution dictionary is first read out to the optimiation variable string that describes the path. Afterward, the substrings that represent the directions are translated to linguistic terms, such that the solution can be printed in a human-readable format. Based on these two steps, the validation process checks if any constraints are violated. If constraints are violated warnings will be shown in the output with some tuning suggestions.
+ +def read_validate_solution(solution: dict, num_turns: int, num_dim: int):
+
+ '''
+ Purpose:
+ To validate the solution returned by the solver. Make it readable, and analyze if it makes sense.
+ Inputs:
+ 1. solution: The solution results dictionary which is returned by the solver (results["configuration"]).
+ 2. num_turns: The number of turns for the simulation.
+ 3. num_dim: The number of dimensions, which is 3 for this sample (3D).
+ Outputs:
+ 1. valid: A boolean variable that specifies the validity of the solution.
+ 2. pos_dit: Layered position dictionary that contains all of the nodes' locations per turn {turn: {x: x_pos, y:y_pos, z:z_pos}}.
+ 3. dir_dict: Dictionary containing the linguistic interpretation of the 6 directions.
+ 4. var_dict: Dictionary containing the spin per optimization variable.
+ 5. x_arr: Array of x positions.
+ 6. y_arr: Array of y positions.
+ 7. z_arr: Array of z positions.
+ '''
+
+ print('\n')
+ valid = True
+ move = ''
+ sol_str = ''
+ x_arr = [0]
+ y_arr = [0]
+ z_arr = [0]
+ dir_dict = {'100':'out','010':'in','001':'right','110':'left','101':'up','011':'down'}
+ pos_dict = {0:{"x":0, "y":0, "z":0}}
+ var_dict = {}
+ for key,val in solution:
+ if key<(num_turns*num_dim):
+ turn = floor(key/num_dim)+1
+ print("Turn: "+str(turn),"var: "+str(key),"spin: "+str(val))
+ var_dict |= {str(key): val}
+ if key%3<num_dim and key<(num_turns*num_dim):
+ move = move+str(val)
+ if move in dir_dict:
+ x_pos = pos_dict[turn-1]["x"]
+ y_pos = pos_dict[turn-1]["y"]
+ z_pos = pos_dict[turn-1]["z"]
+ if dir_dict[move] == "out":
+ x_pos += 1
+ elif dir_dict[move] == "in":
+ x_pos += -1
+ elif dir_dict[move] == "right":
+ y_pos += 1
+ elif dir_dict[move] == "left":
+ y_pos += -1
+ elif dir_dict[move] == "up":
+ z_pos += 1
+ elif dir_dict[move] == "down":
+ z_pos += -1
+ new_pos = {turn: {"x":x_pos, "y":y_pos, "z":z_pos}}
+ pos_dict |= new_pos
+ x_arr += [x_pos]
+ y_arr += [y_pos]
+ z_arr += [z_pos]
+ #print move
+ print(dir_dict[move], '\n')
+ sol_str = sol_str + " "+ dir_dict[move]
+ move = ''
+
+ elif move == ('111' or '000'):
+ valid = False
+ print('Illegal move')
+
+ print("solution:", sol_str, '\n')
+ print("positions: ", pos_dict)
+
+ incorrect_pos = []
+ for ref in range(0,len(pos_dict)):
+ for tar in range(ref+1,len(pos_dict)):
+ if pos_dict[ref] == pos_dict[tar]:
+ print(f"Invalid position. Position already taken. Positions {ref}-{tar}, (if difference = 2, increase lambda_4, otherwise increase lambda_1)")
+ valid = False
+ incorrect_pos += [ref,tar]
+
+ print('\n Incorrect positions',incorrect_pos)
+ return valid, pos_dict, dir_dict, var_dict, x_arr, y_arr, z_arr
+Below you can find the function that plots the self-avoiding path. Probably this is the easist way to verify whether the path is correct or not. The initial position is plotted with a red "+" and the positions with a green dot, inteconnected by the lines.
+ +def plot_path(x_arr,y_arr,z_arr):
+
+ '''
+ Purpose:
+ This function plots the path of the self-avoiding walk through the position arrays.
+ Inputs:
+ 1. x_arr: the array of x positions.
+ 2. y_arr: the array of y positions.
+ 3. z_arr: the array of z positions.
+ Outputs:
+ 1. The 3D plot of the self-avoiding path.
+ '''
+ print('\n')
+ print('x_arr: ', x_arr)
+ print('y_arr: ', y_arr)
+ print('z_arr: ', z_arr)
+
+ fig = plt.figure()
+ ax = plt.axes(projection ='3d')
+ # plotting
+ ax.plot3D(x_arr, y_arr, z_arr, 'g-o')
+ ax.set_title('3D self-avoiding path')
+ plt.plot(0,0,0,'r+') # plot initial position
+ plt.show()
+Finally, time to submit to the Azure solvers! For this notebook a small problem is chosen because the execution times grow exponentially with the number of positions (path length). A self-avoiding path of 10 positions needs to be found.
+First it is necessary to build the entire optimization problem, comprising of the constraint functions. The optimization problem is then sent to a specified Azure QIO solver, such as simulated annealing or parallel tempering (heuristic solvers), which try to find a good solution to the problem. The the nonlinear optimization constraint weights and the solvers need tuning to find a valid path, which is a difficult process. To see a valid path, you can first run the with the provided parameters, as the solver and the cost function have been tuned to find these values. Other than that, feel free to play around with the values. Note that it is easier to tune for smaller problem sizes (nodes/positions) since the simulation times are smaller, however, the parameters do not map well to larger problems. Additionally, a suggestion is to first use the parameter-free solvers to limit the tuning work to the constraint weights.
+At the end of the notebook, a few resources are given where you can find more info about the QIO solvers.
+Below I define a function 'submit' which construct the optimization function, submits to a solver, and parses the results. If you want to re-run the optimization, you might first need to close the figure!
+ +def submit():
+
+ '''
+ Purpose:
+ This function submits to the Azure solvers.
+ Adjust the variables and solver properties here.
+ This functions coordinates the entire simulation by calling the previously defined functions.
+ '''
+
+ nodes = "HHHHHHHHHHHHHHH" # The 15 nodes that will be analyzed.
+ len_seq = len(nodes) # The length of the sequence.
+ num_turns = len(nodes)-1 # The number of turns between the nodes.
+ num_dim = 3 # Number of optimization variables required to describe a turn
+ lambda_0 = 1
+ lambda_2 = 100
+ lambda_3 = 100
+ lambda_4 = 3
+
+ terms = ( distance_constraint(num_turns, num_dim, lambda_0)+penalize_000(len_seq, num_dim, lambda_2)+
+ penalize_111(len_seq, num_dim, lambda_3)+no_return_constraint(num_turns, num_dim, lambda_4) )
+
+ problem = Problem(name="Self Avoiding Walk Problem", problem_type=ProblemType.pubo, terms=terms)
+
+ # Submit
+ start = time.time()
+
+ ##### Two parameter-free solvers - these are tuned for you after submission.
+ #parameter free simulated annealing
+ #solver = SimulatedAnnealing(workspace, timeout=600)
+ #parameer free parallel tempering
+ #solver = ParallelTempering(workspace, timeout=600)
+
+ ##### Parametrized simulated annealing - needs to be tuned by you.
+ solver = SimulatedAnnealing(workspace, sweeps=50, restarts=2000, beta_start=1e-8, beta_stop=10, timeout=3600)
+
+
+ print('Submitting problem...')
+ job = solver.submit(problem)
+
+ while job.details.status != 'Succeeded' and job.details.status != 'Failed':
+ job.refresh()
+ print('waiting')
+ time.sleep(10)
+
+ # Results
+ results = job.get_results()
+ print('\n Results:', results)
+ config = results["configuration"]
+
+ time_elapsed = time.time() - start
+ print("Execution time in seconds: ", time_elapsed, '\n')
+
+ solution = [(int(k),v) for k, v in config.items()]
+ solution.sort(key=lambda tup: tup[0])
+
+ valid, posDict, dirDict, configDict, x_arr, y_arr, z_arr = read_validate_solution(solution, num_turns, num_dim)
+
+ plot_path(x_arr, y_arr, z_arr)
+
+submit()
+In this notebook we'll cover how to formulate an optimization problem to find a path through a 3D lattice that does not cross itself (self-avoiding). The optimization problem is solved with the Azure Quantum service which contains numerous heuristic (non-linear) optimization solvers. Finding a self-avoiding path in a 3D lattice can be considered a difficult problem, in the sense that it has suffers from exponential scaling conditioned on the number of turns and the number of dimensions (2D vs 3D, etc.).
-In this notebook, a step by step approach is taken to explain the function definitions to are necessary to create the optimization problem.
-Goal: To find a path in a 3D lattice that does not cross itself.
-To clarify some vocab and assumptions used in the notebook:
-Note: Because the number of terms to for this optimization problem grows explosively with the number of turns, you might want to try running this in the Azure Quantum notebooks. The more technical reason for this is that expansions of nonlinear terms have to be computed locally resulting in massive number of terms needing to be uploaded. Unless you are very patient, a recommendation would be to either check out the online notebook experience or rewrite this notebook's problem class to its streaming counterpart (https://docs.microsoft.com/en-us/azure/quantum/optimization-streaming-problem).
- -import time
-from math import floor, log2
-from azure.quantum import Workspace
-from azure.quantum.optimization import Problem, ProblemType, Term, ParallelTempering, Tabu, SimulatedAnnealing
-from azure.identity import ClientSecretCredential
-from mpl_toolkits import mplot3d
-import numpy as np
-import matplotlib.pyplot as plt
-To run the sample, you'll need to have a quantum workspace. Check out this module if you don't have one yet: https://docs.microsoft.com/en-us/learn/modules/get-started-azure-quantum/. -Fill in the variables below.
- -workspace = Workspace(
- subscription_id = "",
- resource_group = "",
- name = "",
- location = "",
- credential = ClientSecretCredential(tenant_id="",
- client_id="",
- client_secret="")
-)
-Throughout the notebook you'll probably realize that printing the "terms" isn't going to provide much help, mainly because there will be too many of them and their dictionary format. A term is the way a mathematical term is expressed in terms of a dictionary for the SDK (see examples). Below a function is defined that prints the term(s) dictionaries that describe the constraint/cost function in terms of "q" optimization variables. You can call the function to check if the working is correct, or if you're unsure of what a constraint looks like. Note that for a large problem it nevertheless becomes difficult to understand the entire function output because of the number of terms.
-An example: -{'c': 1, 'ids': [0, 10, 11, 12]} ---> $1q_0q_{10}q_{11}q_{12}$
- -def print_function(terms: list):
-
- '''
- Purpose:
- Takes a list of terms and prints it out as a mathematical (cost/constraint) function.
- Example:
- {'c': 1, 'ids': [0, 10, 11, 12]} ---> 1q_0q_10q_11q_12
- Inputs:
- 1. terms: the list of terms.
- '''
-
- k = 0
- final_string = ''
- final_string = ''
- for term in terms:
- term = term.to_dict()
- weight = term['c']
- ids = term['ids']
- string = '('
- if weight >= 0:
- if k == 0:
- string = '('+str(weight)
- else:
- string = '+' + '('+str(weight)
- if weight < 0:
- if k == 0:
- string = '(' + str(abs(weight))
- else:
- string = '-' + '(' + str(abs(weight))
- for id_ in ids:
- string = string + f'q_{id_}'
- string = string + ')'
- final_string = final_string + string
- k += 1
-
- print('[ ' + final_string + ' ]')
-Throughout the notebook we will need to cross-multiply nonlinear terms, for computing squares for example. The function below performs the expansion. An important note to keep in mind here, especially if you're looking to optimize the code, is that the function does not assemble identical terms. Since the cost functions dealt with in this notebook are still rather small, reducing the number of terms was left out of the scope.
-More clearly stated in terms of equations:
-$$ \text{Current function (see below): } \hspace{0.2cm} (a+b)^2 = a^2 + ab + ab + b^2 $$$$ \text{Optimal function: } \hspace{0.2cm} (a+b)^2 = a^2 + 2ab + b^2 $$Should be a fun programming exercise to reduce the number of terms significantly!
- -def cross_multiply(list_a: list, list_b: list) -> list:
-
- '''
- Purpose: Cross multiplies two lists of terms (linear ex. [q_0 + q_1] or non-linear ex. [q_0q_1]) to return the expansion.
- Can compute powers of groups this way (^2, ^3,...), like squaring a list of terms.
- Calculates the expansion locally, unlike the SlcTerm class.
- Example: (2q_0q_1+3q_2q_3)^2 => (2q_0q_1)^2 + 6q_0q_1q_2q_3 + 6q_0q_1q_2q_3 + (3q_2q_3)^2
- Input:
- 1. list_a: list which serves as the reference list (first 'for' loop).
- 2. list_b: list which serves as the target list (second 'for' loop).
- Output:
- 1. list of term objects.
- '''
-
- terms = []
- for one in list_a:
- for uno in list_b:
- alpha = one.to_dict()
- beta = uno.to_dict()
- weight = int(alpha['c'])*int(beta['c'])
- ids = list(alpha['ids'])+list(beta['ids'])
- terms += [Term(c = weight, indices = ids)]
- return terms
-While running through the notebook you might want to execute and print some of the constraint functions. To do that, you'll need some variables to be defined. Only 3 turns are considered with the defined parameters below, as that will keep the function outputs small and understandable!
- -nodes = "ABCD" # The four position/node names.
-len_seq = len(nodes) # The length of the node
-num_turns = len(nodes)-1 # The number of turns, which is one less than the number of therefore -1 turns.
-num_dim = 3 # Number of optimization variables required to describe a turn
-lambda_0 = 1 # Penalty weight for the 'distance constraint'
-lambda_1 = 40 # Penalty weight for the 'no return constraint'
-lambda_2 = 25 # Penalty weight for the invalid direction '000' constraint
-lambda_3 = 25 # Penalty weight for the invalid direction '111' constraint
-Prior to starting to define our cost function, we have to construct some basic functions and agree on some definitions. In this part, we'll define the 6 directions of the x-y-z plane in terms of the optimization variables ("q").
-A common way to encode a turn based optimization problem is through a direction-based string, in which each optimization variable substring (ex. $q_0q_1q_2$) denotes some kind of decision of which direction was takn. For example, take a car take can only move backward or forward. The decision of the driver at some time $t$ can be described as $q_t \in \{0, 1\}$, where 0 represents going backward and 1 forward, respectively. Over a multiple time steps, say 3, one can then describe the movement of the car by:
-$$ \text{Movements:} \hspace{0.25cm} q_0q_1q_2 $$The sequence of movements, "forward, backward, backward", can then be represented as:
-$$ \text{Sequence of movements:} \hspace{0.25cm} q_0q_1q_2 = 100$$Alright, if that is clear to you, then the next step is to represent the 6 directions of the x-y-z plane in such a fashion. Regrettably, a "q" can only be either 0 or 1, not 2, 3, 4, 5 to represent the six directions. But what can be done is to use multiple optimization variables to describe a direction! For example:
-$$ \text{Sequence of movements:} \hspace{0.25cm} q_0q_1q_2q_3q_4q_5q_6q_7q_8 $$$$ \text{First move:} \hspace{0.25cm} q_0q_1q_2 $$$$ \text{Second move:} \hspace{0.25cm} q_3q_4q_5 $$$$ \text{Third move:} \hspace{0.25cm} q_6q_7q_8 $$Because we need to encode 6 directions, we need at least 3 optimization variables to form unique substrings. The table explains this, we to have enough possible combinations to encode the 6 directions:
-Note, for later, you should keep in mind that 2 of the 8 combinations become irrelevant. It is important that the solver does not return those strings!
-Perfect! Now we just need to assign each direction a unique substring. Below you can find my choices, feel free to change these if you're starting from scratch.
-Now we just need to define direction variables. These are as defined in the above table. These variables should take a value 1 if a movement is made in that direction, while the other directions must all be zero. This can be achieved by the following scheme, where $k$ denotes the turn, and $\gamma = 3(k-1)$ (3 because 3 dimensions):
-To clarify a bit further through an example: -If in the first turn (k=0) the solver returns $q_0=1$, $q_1=0$, $q_2=0$, and in the second turn (k=1) returns $q_3=0$, $q_4=0$, $q_5=1$, then a "+x and +y" were taken, respecively.
-These direction variables will help understand some difficult constraints later on. Defining these direction variables additionally makes the code much more readable, since we can define everything in terms of the direction variables instead of the individual optimization variables.
-Below you can find the function for the direction variables. The formulas have been expanded (mathematically) to easily declare them in a number of terms, however it is hard-coded this way. Some inputs to the function ("sign_dir", "signpos", "lambda") are relevant for building constraints later in the notebook, but have to defined in this function. Further information can be found in the function docstring.
- -def direction_variables(direction: str, offset: int, sign_dir: int, sign_pos: int, lambda_: int) -> list:
-
- '''
- Purpose:
- Translates the direction (+x,-x,+y,-y,+z,-z) of turn 'i' as a function of three optimization variables. (Three q's because of the defined coordinate system).
- Example:
- Direction "+z" in the first turn (turn = 1) is translated to: q_{0+offset}q_{2+offset}-q_{0+offset}q_{1+offset}q_{2+offset}.
- Inputs:
- 1. direction: A direction from an x-y-z coordinate system, one of the following: ('+x','-x','+y','-y','+z','-z').
- 2. offset: Offset gives the turn number expressed in the first "q" of that turn. Equal to gamma in the explanation.
- Example: Turn 1 starts with "q" q_0, offset=0. Turn 2 starts with q_3, offset = 3.
- 3. sign_dir: Changes the sign of the weights corresponding to negative directions "-x", "-y", "-z" -> necessary for finding the positions, for exmaple (+x) "-" (-x).
- 4. sign_pos: Changes the sign of the weights corresponding to negative positions "-(x,y,z)" -> necessary for finding the distances between node "i" and node "j".
- 5. lambda_: The weight term associated with a constraint.
- Output:
- 1. A list of term objects.
- '''
-
- terms = []
- if direction == "+x":
- term_0 = Term(c= 1*sign_pos*lambda_, indices=[0+offset])
- term_1 = Term(c=-1*sign_pos*lambda_, indices=[0+offset, 1+offset])
- term_2 = Term(c=-1*sign_pos*lambda_, indices=[0+offset, 2+offset])
- term_3 = Term(c= 1*sign_pos*lambda_, indices=[0+offset, 1+offset, 2+offset])
- terms = [term_0, term_1, term_2, term_3]
- elif direction == "-x":
- term_0 = Term(c= 1*sign_dir*sign_pos*lambda_, indices=[1+offset])
- term_1 = Term(c=-1*sign_dir*sign_pos*lambda_, indices=[0+offset, 1+offset])
- term_2 = Term(c=-1*sign_dir*sign_pos*lambda_, indices=[1+offset, 2+offset])
- term_3 = Term(c= 1*sign_dir*sign_pos*lambda_, indices=[0+offset, 1+offset, 2+offset])
- terms = [term_0, term_1, term_2, term_3]
- elif direction == "+y":
- term_0 = Term(c= 1*sign_pos*lambda_, indices=[2+offset])
- term_1 = Term(c=-1*sign_pos*lambda_, indices=[0+offset, 2+offset])
- term_2 = Term(c=-1*sign_pos*lambda_, indices=[1+offset, 2+offset])
- term_3 = Term(c= 1*sign_pos*lambda_, indices=[0+offset, 1+offset, 2+offset])
- terms = [term_0, term_1, term_2, term_3]
- elif direction == "-y":
- term_0 = Term(c= 1*sign_dir*sign_pos*lambda_, indices=[0+offset, 1+offset])
- term_1 = Term(c=-1*sign_dir*sign_pos*lambda_, indices=[0+offset, 1+offset, 2+offset])
- terms = [term_0, term_1]
- elif direction == "+z":
- term_0 = Term(c= 1*sign_pos*lambda_, indices=[0+offset, 2+offset])
- term_1 = Term(c=-1*sign_pos*lambda_, indices=[0+offset, 1+offset, 2+offset])
- terms = [term_0, term_1]
- elif direction == "-z":
- term_0 = Term(c= 1*sign_dir*sign_pos*lambda_, indices=[1+offset, 2+offset])
- term_1 = Term(c=-1*sign_dir*sign_pos*lambda_, indices=[0+offset, 1+offset, 2+offset])
- terms = [term_0, term_1]
- return terms
-
-
-##### ----- Test the function and print output:
-turn = 1 # play with this value!
-dir_var = direction_variables("+x",(turn-1)*num_dim,1,1,lambda_0)
-print('dir_var term dictionary "+x": ', dir_var)
-print_function(dir_var)
-Great that we have the directions defined! Now let's use them and define a function that calculates the difference in positions, which we will need later. After each turn, a direction that appends to the path. A valid path is defined as a one that does not cross itself, meaning that we need to compare the positions over the turns. In other words, every new position appended to the path needs to be checked with all previously held positions, to verify that a position hasn't been visited twice.
-So how do we go about this?
-First we need a way to know the positions after choosing to go in a certain direction, so let's tackle that first. Consider the fact that the directions are recorded for each turn, which already gives some sort of log of the positions in the path. By summing the direction variables respective to their dimension (x/y/z) acoomplishes this, but only if a negative move can compensate a positive one. For example, after two turns the $x$-location can be expressed as:
-$$ x_2 = ( d^{1}_{+x} - d^{1}_{-x} ) + ( d^{2}_{+x} - d^{2}_{-x} )$$NOTE: In the above function ("direction_variables"), the negative sign for the negative direction variables is controlled through "sign_dir". "sign_dir" must be set to -1 to assign a negative sign to it as required in these position formulas!
-Remember that for a turn only one direction variable can be set to 1! Thus if two moves in the $-x$ direction are taken, then:
-$$ x_2 = ( d^{1}_{+x} - d^{1}_{-x} ) + ( d^{2}_{+x} - d^{2}_{-x} ) = ( 0 - 1 ) + ( 0 - 1 ) = -2 $$As with the simple example for the x-dimension, the same can be applied to the $x$ and $y$ directions. The formulas below descrive the position after some number of turns through a summation (for the maths/physics enthousiasts, these are integrals of the velocities to derive the positions :) ).
-$$ x^k = \sum_{k=1}^{k} ( d^{k}_{+x} - d^{k}_{-x} ) $$$$ y^k = \sum_{k=1}^{k} ( d^{k}_{+y} - d^{k}_{-y} ) $$$$ z^k = \sum_{k=1}^{k} ( d^{k}_{+z} - d^{k}_{-z} ) $$Great stuff! Now we have a method to find the positions along the path. To find the differences in position, for example between turn 2 and turn 1, it is only necessary to apply a subtraction:
-$$ x_2 - x_1 = \sum_{k=1}^{2} ( d^{k}_{+x} - d^{k}_{-x} ) - \sum_{k=1}^{1} ( d^{k}_{+x} - d^{k}_{-x}) = ( d^{2}_{+x} - d^{2}_{-x} ) $$$$ y_2 - y_1 = \sum_{k=1}^{2} ( d^{k}_{+y} - d^{k}_{-y} ) - \sum_{k=1}^{1} ( d^{k}_{+y} - d^{k}_{-y}) = ( d^{2}_{+y} - d^{2}_{-y} ) $$$$ z_2 - z_1 = \sum_{k=1}^{2} ( d^{k}_{+z} - d^{k}_{-z} ) - \sum_{k=1}^{1} ( d^{k}_{+z} - d^{k}_{-z}) = ( d^{2}_{+z} - d^{2}_{-z} ) $$Generalizing this to any difference between two positions:
-$$ x_j - x_i = \sum_{k=i}^{j} ( d^{k}_{+x} - d^{k}_{-x} ) $$$$ y_j - y_i = \sum_{k=i}^{j} ( d^{k}_{+y} - d^{k}_{-y} ) $$$$ z_j - z_i = \sum_{k=i}^{j} ( d^{k}_{+z} - d^{k}_{-z} ) $$Fantastic! Now on to the function code for the difference in positions. Do read the function and its docstring to understand the implementation, because it reference the previously defined function to make use of direction variables.
- -def diff_in_pos(start_turn: int, end_turn: int, num_dim: int, lamda_: int):
-
- '''
- Purpose:
- Expresses the difference in position (x,y,z) between two turns as a function of the q's encoding of the directions.
- In other words, expresses the difference in the position after the start_turn and position after the end_turn.
- Example:
- Difference between turn 0 (no turns yet, initial position = (0,0,0)) and turn 2:
- x(2) = [ move(+x, turn 1) - move(-x, turn 1) ] + [ move(+x, turn 2) - move(-x, turn 2) ]
- Note: only one 'move' per turn gets activated as they are represented by the same q's (ex. turn 1 is represented by q_0, q_1, and q_2).
- Inputs:
- 1. start_turn: the initial (reference) turn.
- 2. end_turn: the final (target) turn.
- 3. num_dim: the number of dimensions (3, x-y-z coordinate system).
- Outputs:
- 1. The difference in the x direction.
- 2. The difference in the y direction.
- 3. The difference in the z direction.
- '''
-
- x_diff = y_diff = z_diff = []
- if start_turn < end_turn and start_turn >= 0 and end_turn >= 1:
- for turn in range(start_turn,end_turn+1):
- x_diff += direction_variables("+x",(turn-1)*num_dim,1,1,lamda_)+direction_variables("-x",(turn-1)*num_dim,-1,1,lamda_)
- y_diff += direction_variables("+y",(turn-1)*num_dim,1,1,lamda_)+direction_variables("-y",(turn-1)*num_dim,-1,1,lamda_)
- z_diff += direction_variables("+z",(turn-1)*num_dim,1,1,lamda_)+direction_variables("-z",(turn-1)*num_dim,-1,1,lamda_)
- return x_diff, y_diff, z_diff
-
-
-##### ----- Test the function and print output:
-start_turn = 1 # play with this value!
-end_turn = 2 # play with this value!
-x_diff, y_diff, z_diff = diff_in_pos(start_turn, end_turn, num_dim, lambda_0)
-print('Difference in x described as term list: \n', x_diff)
-print('\nDifference in x described as mathematical function:')
-print_function(x_diff)
-Fantastic that you've made it this far! With all the above functions defined we can start solving the problem.
-The distance constraint enforces that after each turn, the difference in positions after turns $i$ and $j$ must be larger or equal to 1 (distance between lattice points in the same dimension). -To visualize the scenario this constraint tries to penalize consider a 2D example with 4 turns, starting in $(0,0)$:
-$$ \text{Turn 1, go right: } \hspace{0.3cm} (0,0) => (0,1) $$-$$ \text{Turn 2, go up: } \hspace{0.3cm} (0,1) => (1,1) $$ -$$ \text{Turn 3, go left: } \hspace{0.3cm} (1,1) => (1,0) $$ -$$ \text{Turn 4, go down: } \hspace{0.3cm} (1,0) => (0,0) $$
-As you can see, the last step conflicts with our constraint. It is not permitted to return into a position that we've already been, namely $(0,0)$. Turn 1-3 are all valid, since they remain a distance of 1 away from all other previously held positions.
-Let's define what we want mathematically. tThe distance between the position after turn $i$ and $j$ must be larger or equal to 1. For this we'll need some basic geometry, the 3D variant of the Pythagorean theorem. From the theorem, it can be understood that hypotenuse ($c$) must always be larger or equal to one, since that is what defines the distance between two points. Mathematically speaking working this out:
-$$ C_{i,j} \geq 1 $$$$ \sqrt{(x_{j} - x_{i})^2 + (y_{j} - y_{i})^2 + (z_{j} - z_{i})^2} \geq 1 $$The square-root is a bit of a headache when writing out a cost function, it would require an approximation or some sort of factorization, not fun if you're dealing with many terms. To avoid that, we're simply going to square both sides, which in our case if fine (for maths people), since we aren't dealing with negative numbers here:
-$$ C_{i,j}^2 \geq 1 $$$$ {(x_{j} - x_{i})^2 + (y_{j} - y_{i})^2 + (z_{j} - z_{i})^2} \geq 1 $$Nice! If you're sharp, you'll start seeing where the previous functions might come in handy now. The function that computes the difference between the respective positions ("diff_in_pos") in each dimension was defined above, but also the function that computes the an expansion ("cross_multiply") was defined at the beginning of the notebook. Writing those function calls out in order :
-Alright, you might be wondering how to deal with the $>=$, since it is necessary to have an equation that equals zero in order to integrate it into the cost function. -This is where slack variables come in. The purpose of slack variables is to convert an inequality constraint into an equality constraint. Or in other words, for this equation, we need to compensate some values (explained further below) in order to convert the equation to an equality constraint. Consider the constraint in its simplest form:
-$$ C_{i,j}^2 \geq 1 $$It is clear that we would need to add some number $V_{i,j}$, on the right hand side of the equation in order to make this an equality constraint. $V_{i,j}$ can be any positive value and is conditioned on $i$ and $j$:
-$$ C_{i,j}^2 = 1+V_{i,j}$$Let's say that instead of representing a single value we want V_{i,j} to represent a range of values $[0-4]$. This is where slack variables can be introduced. By introducing additional optimization variables, which have nothing to do with the path encoding, the range can be described as following:
-$$ V_{i,j} = [0-4] = q_{s1}+q_{s2}+q_{s3}+q_{s4} $$Or more compactly:
-$$ V_{i,j} = [0-4] = q_{s1}+q_{s2}+2q_{s3} $$So any value in the range $[0-4]$ is described through these q's. You can derive 1 by assigning a single "q", which has weight 1, the value 1. Likewise, you can get the value 4 by assigning all optimization variables the value 1.
-Alright, so how do we know which range of values we need to represent? That is revealed by the turn numbers, $i$ and $j$. Over these turns, the maximum squared distance that can be achieved is when moves are only made in a single direction, the "+x" direction for example. The maximum distance is then equal to the number of moves squared $(j-i)^2$.
-$$ C_{i,j}^2 = (j-i)^2 $$Then writing out the maths gives us the necessary upper bound for the range that the optimization variables need to represent, where $S$ stands for the number of slack variables necessary:
-$$ \text{Upper bound: } \hspace{0.5cm} C_{i,j}^2 = (j-i)^2 = 1 + \sum_{s=0}^{S}q_s$$$$ \text{Upper bound: } \hspace{0.5cm} \sum_{s=0}^{S}q_s = (j-i)^2 - 1 $$For the lower bound of the range, the value is 0. The reason for that is we want to be able to represent any number equal to or larger than 1 for the inequality constraint through the addition of the slack variables. Concluding, the range the q's must represent is $[0,((j-i)^2-1)]$. More on how to calculate the number of slack variables necessary to define this range later, as you don't need ((j-i)^2-1) optimization variables!
-Piecing all the parts together we find the following equation for a $i$-$j$ combination:
-$$ {(x_{j} - x_{i})^2 + (y_{j} - y_{i})^2 + (z_{j} - z_{i})^2} = 1 + \sum_{s=0}^{S}q_s $$ -Because constants are irrelevant for the optimization landscape the value 1 can be neglected. The reason being that constants only introduce linear offsets, thus impacting the entire optimization equally. The distance constraint which needs to account for all $i$-$j$ combinations, where $j$>$i$ and $N-1$ the number of turns, is summarized as follows:
-$$ \sum_{i=1}^{N-1} \sum_{j>i}^{N-1} \lambda_0 ( {(x_{j} - x_{i})^2 + (y_{j} - y_{i})^2 + (z_{j} - z_{i})^2} - \sum_{s=0}^{S}q_s ) $$ -The function that builds this constraint is presented below. Read the docstring for how it works.
-Note, to run this function you'll need to run the next function cell to define 'generate_slack_coefficients'.
def distance_constraint(num_turns: int, num_dim: int, lambda_0: int) -> list:
-
- '''
- Purpose:
- Build the distance contraint based on previosly defined functions.
- Constraint: Distance squared between i and j must be larger or equal to 1.
- Constraint: L_{i,j}^2 >= 1 => L_{i,j}^2 = 1+q_{slacker} (converting the inequality constraint to an equality constraint.)
- Example/Explanation:
- L{1,2}^2 = (x1-x2)^2 + (y1-y2)^2 + (z1-z2)^2
- x1 = [move(+x1)-move(-x1)]
- x2 = [move(+x1)-move(-x1)] + [move(+x2)-move(-x2)]
- <same for other dimensions >
- <fill into first line>
- L{1,2}^2 = [move(+x2)-move(-x2)]^2 + [move(+y2)-move(-y2)]^2 + [move(+z2)-move(-z2)]^2
- [move(+x2)-move(-x2)]^2 + [move(+y2)-move(-y2)]^2 + [move(+z2)-move(-z2)]^2 - q_{slacker} = 0 ---> expressed in q's == cost function
- Inputs:
- 1. num_turns: the number of turns.
- 2. num_dum: the number of dimensions.
- 3. lambda_0: the constraint weight for the distance constraint.
- Output:
- 1. List of term objects describing the distance constraint.
- '''
-
- terms = []
- slack_indexer = 0
- for start_turn in range(1,num_turns+1):
- for end_turn in range(start_turn+1,num_turns+2):
- # Calculate the differences in positions for each dimension.
- x_diff_i_j, y_diff_i_j, z_diff_i_j = diff_in_pos(start_turn, end_turn, num_dim, lambda_0)
- # Compute the squared distance (Pythagorean theorem) by calculating the squared expansion.
- x_diff_i_j_2, y_diff_i_j_2, z_diff_i_j_2 = cross_multiply(x_diff_i_j,x_diff_i_j), cross_multiply(y_diff_i_j,y_diff_i_j), cross_multiply(z_diff_i_j,z_diff_i_j)
- # Add slack variables due to inequality constraint.
- slack_var_terms = []
- slack_coefficients = generate_slack_coefficients(end_turn-start_turn)
- for s in range(0,len(slack_coefficients)):
- slack_var_terms += [Term(c=-slack_coefficients[s], indices=[num_turns*num_dim+slack_indexer+s])]
- terms += x_diff_i_j_2 + y_diff_i_j_2 + z_diff_i_j_2 + slack_var_terms
- slack_indexer+=len(slack_coefficients)
- return terms
-
-##### ----- Test the function and print output:
-# The output is too large to understand, nevertheless you can run the below statements to view how visualize the term-scaling of the problem.
-number_turns = 2 #play with this value!
-dist_terms = distance_constraint(number_turns, num_dim, lambda_0)
-
-# Only print below statements if you want to get an impression of the number of terms.
-#print('Distance constraint term dictionaries: ', dist_terms) # not readable!
-#print('\nDistance constraint: ')
-#print_function(dist_terms) # not readable!
-To remove any human calculations to construct the distance constraint(s), we'll automate the computation of the number of slack variables and their weights.
-As you saw in the previous part, certain slack variables can be assigned weights to reduce the number of slack variables necessary. Recall that the range $[0,4]$ can be defined by the following two equations with slack variables (s):
-$$ V_{i,j} = [0-4] = q_{s1}+q_{s2}+q_{s3}+q_{s4} $$$$ V_{i,j} = [0-4] = q_{s1}+q_{s2}+2q_{s3} $$The bottom equation is more compact because the last variables carries a weight-2. In Lucas' paper (see 2.4 of arxiv:1302.5843) the algorithm to find such weights is presented, which is based on calculating all $2^{x}$, where $x$ iterates from zero to the log_2 of maximum value in the range (4 in this example).
-Below are two examples that may help visualize why and how the slack variables represent the ranges of values. -\ -
-\ -
-\ -
-The function is given below. In context of the distance constraint, we're still faced with a minor problem. If $(j-i)^2-1$ equals 0, then we don't need any slack variables. The reason for this is that the range becomes $[0,0]$, which can also be represented by no variables at all.
- -def generate_slack_coefficients(turn_diff: int):
-
- '''
- Purpose:
- Calculates the number of slack variables and their weights.
- Example:
- For the constraint: x1 + x2 + x3 + x4 <= 4 which is converted to x1 + x2 + x3 + x4 + s1 + s2 +2*s3 = 4,
- with 's' being slack variables. This function computes the weights of these slack variables ([1,1,2] for the example).
- Reference:
- Lucas' paper (see 2.4 of arxiv:1302.5843)
- Input:
- 1. turn_diff: the differences in turns (end_turn - start_turn).
- Output:
- 1. y: the weights of the slack variables.
- '''
-
- dist_diff = (turn_diff**2)-1
- if dist_diff == 0:
- y = [] # no slack variables needed
- elif dist_diff > 0:
- M = floor(log2(dist_diff))
- y = [-2**n for n in range(M)]
- y.append(-(dist_diff + 1 - 2**M))
- return y
-
-##### ----- Test the function and print output:
-turn_difference = 3 # play with this value!
-slack_weights = generate_slack_coefficients(turn_difference)
-print('Slack weights: ', slack_weights)
-The distance constraint enforces a minimum distance of 1 between the different taken positions. However, for two consecutive turns, it does not completely prevent (to a high enough degree) going to a previous position. A simplified explanation for this is that the distance constraint is described through the positional differences in $x$, $y$, and $z$. Additionally, for consecutive turns, the number of slack variables is zero ($(i-j)^2-1 = 0$ for consecutive turns} - thus there is no auxiliary benefit introduced by slack variables to compensate for being in a different position after the turn. Therefore, the solver will try to keep the positions as close to each other as possible, which leads to overlapping positions, especially in consecutive turns. To overcome this, a penalization needs to be introduced that stops returning to the previous position held in the turn before. For example, what we want to stop is the following:
-$$ \text{Turn 1: A to B.} $$$$ \text{Turn 2: B to A.} $$This can be equivalenty described as moving in opposite directions over the consecutive turns. For example, first going in the "+y" direction, and following that moving in the "-y" direction.
-Designing such a constraint is straightforward, the opposite direction variables need to be multiplied. -Preventing a revisit of a position for the "x" direction over the first and second turn can be described by the following constraint:
-$$ {\lambda}_2(d^{1}_{+x} d^{2}_{-x} + d^{1}_{-x}d^{2}_{+x}) $$ -If a move is made in the "+x" direction first, $d^{1}_{+x}$ takes value 1. If afterwards a move is made in the "-x" direction, $d^{2}_{-x}$ also takes value 1, meaning that the combined term -$d^{1}_{+x} d^{2}_{-x}$ also becomes 1, enforcing the constraint with penalty value ${\lambda}_2$.
-Generalizing the idea to any direction and for all turns, we find the "no return constraint":
-$$ {\lambda}_1 \left( \sum_{m}^{\in \{x,y,z\} }\sum_{t=1}^{N-1} d^{t}_{+m} d^{t+1}_{-m} + d^{t}_{-m}d^{t+1}_{+m} \right) $$ -In this constraint, iterations are performed over the different dimensions ($m$), and the turn ($t$). The function defintion for the no return constraint is given below. Make sure to read the docstring. In the function the 'cross_multiply' method is used to expand and calculate all the necessary terms of the constraint. The direction variables are dictionaries of polynomials, and therefore have to be expanded before submitting to the Azure QIO solvers.
- -def no_return_constraint(num_turns: int, num_dim: int, lamda_4: int)-> list:
-
- '''
- Purpose:
- Build the constraint that penalizes going back to the same position/node two turns later.
- Example:
- Node A ---> <move +x> ---> Node B ---> <move -x> ---> Node A => erroneous as we've been there already.
- Two sequential moves may not be in the same dimension and in opposite directions: (+x then -x), (-x then +x), (+y then -y) etc.
- Inputs:
- 1. num_turns: the number of turns.
- 2. num_dim : the number of dimensions, which is 3 for this sample.
- 3. lambda_4 : the penalty weight for this constraint.
- Outputs:
- 1. List of term objects.
- '''
-
- terms = []
- for i in range(0,num_turns):
- x_out_in = cross_multiply(direction_variables("+x",i*num_dim,1,1,lamda_4), direction_variables("-x",(i+1)*num_dim,1,1,lamda_4))
- x_in_out = cross_multiply(direction_variables("-x",i*num_dim,1,1,lamda_4), direction_variables("+x",(i+1)*num_dim,1,1,lamda_4))
- y_right_left = cross_multiply(direction_variables("+y",i*num_dim,1,1,lamda_4), direction_variables("-y",(i+1)*num_dim,1,1,lamda_4))
- y_left_right = cross_multiply(direction_variables("-y",i*num_dim,1,1,lamda_4), direction_variables("+y",(i+1)*num_dim,1,1,lamda_4))
- z_up_down = cross_multiply(direction_variables("+z",i*num_dim,1,1,lamda_4), direction_variables("-z",(i+1)*num_dim,1,1,lamda_4))
- z_down_up = cross_multiply(direction_variables("-z",i*num_dim,1,1,lamda_4), direction_variables("+z",(i+1)*num_dim,1,1,lamda_4))
- terms += x_out_in + x_in_out + y_right_left + y_left_right + z_up_down + z_down_up
- return terms
-
-
-##### ----- Test the function and print output:
-number_turns = 1 # play with this value!
-nrc_terms = no_return_constraint(number_turns, num_dim, lambda_1)
-print('No return constraint term dictionaries: ', nrc_terms)
-print('\n No return constraint function:')
-print_function(nrc_terms)
-Recall that there are two invalid optimization variable substrings that are not associated with any move, $000$ and $111$. A path that includes these substrings is invalid. To prevent the solver from generating these invalid directions we'll need to add two constraints. Each constraint is specific to an invalid direction, and because we're working directly with the substring we don't need the abstraction layers used in previous constraints, like directional and positional variables.
-Let's first work out the constraint for preventing substring $000$. $ q_{0+\gamma} q_{1+\gamma} q_{2+\gamma} $ may not equal $000$, where $\gamma = 3(k-1)$ and $k$ the turn number. If we were to design the constraints as $ q_{0+\gamma} q_{1+\gamma} q_{2+\gamma} $ there won't be any penalty, because all $q$'s will equal and multiply to zero. Therefore, if all these $q$'s equal zero, we need the substring to multiply to 1. This is achieved by the following, which multiplies to zero if all optimization variables take the value zero:
-$$ {\lambda}_2(1-q_{0+\gamma}) (1-q_{1+\gamma}) (1-q_{2+\gamma}) $$$$ \text{If } \hspace{0.1cm} q_{0+\gamma} = q_{1+\gamma} = q_{2+\gamma} = 0, \hspace{0.1cm} \text{then: } $$-$$ {\lambda}_2(1-q_{0+\gamma}) (1-q_{1+\gamma}) (1-q_{2+\gamma}) = {\lambda}_2 \cdot 1 $$
-Nice! By expanding the equation and neglecting constant terms a constraint is found that can be implemented:
-$$\text{Constraint for 000: }-q_{0+\gamma}−q_{1+\gamma}−q_{2+\gamma}+q_{0+\gamma}q_{1+\gamma}+q_{0+\gamma}q_{2+\gamma}+q_{1+\gamma}q_{2+\gamma}−q_{0+\gamma}q_{1+\gamma}q_{2+\gamma}$$Before going to the function definition, let's first look at penalizing the substring $111$. Luckily this is much easier as the substring multiplies to 1 if all the $q$'s have the value zero. The constraint is therefore simple to derive:
-$$ {\lambda}_3 (q_{0+\gamma} q_{1+\gamma} q_{2+\gamma}) $$$$ \text{If } \hspace{0.1cm} q_{0+\gamma} = q_{1+\gamma} = q_{2+\gamma} = 1, \hspace{0.1cm} \text{then: } $$-$$ {\lambda}_3 (q_{0+\gamma} q_{1+\gamma} q_{2+\gamma}) = {\lambda}_3 \cdot 1$$
-The constraint for the substring $111$ is:
-$$ \text{Constraint for $111$: }\hspace{0.1cm} {\lambda}_3 (q_{0+\gamma} q_{1+\gamma} q_{2+\gamma}) $$The function definitions for these two constraints are given below.
- -def penalize_000(len_seq: int, num_dim: int, lambda_2: int) -> list:
-
- '''
- Purpose:
- Build the constraint that penalizes the invalid moves associated with the 3 q's string: '000'.
- The constraint needs to be defined over the set of all turns, as defined by the 'for' loop.
- Example:
- If turn 2 (q_3q_4q_5) equals '000', assign a large penalty.
- Inputs:
- 1. len_seq: the number nodes to consider.
- 2. num_dim: the number of dimensions (which is 3).
- 3. lambda_2: the penalty weight for this constraint.
- Outputs:
- 1. List of term objects.
- '''
-
- terms = []
- for k in range(0,len_seq):
- offset = k*num_dim
- term_0 = Term(c=-1*lambda_2,indices=[0+offset])
- term_1 = Term(c=-1*lambda_2,indices=[1+offset])
- term_2 = Term(c=-1*lambda_2,indices=[2+offset])
- term_3 = Term(c= 1*lambda_2,indices=[0+offset,1+offset])
- term_4 = Term(c= 1*lambda_2,indices=[0+offset,2+offset])
- term_5 = Term(c= 1*lambda_2,indices=[1+offset,2+offset])
- term_6 = Term(c=-1*lambda_2,indices=[0+offset,1+offset,2+offset])
- terms += [term_0, term_1, term_2, term_3, term_4, term_5, term_6]
- return terms
-
-def penalize_111(len_seq: int, num_dim: int, lambda_3: int) -> list:
-
- '''
- Purpose:
- Build the constraint that penalizes the invalid moves associated with the 3 q's string: '111'.
- The constraint needs to be defined over the set of all turns, as defined by the 'for' loop.
- Example:
- If turn 2 (q_3q_4q_5) equals '111', assign a large penalty.
- Inputs:
- 1. len_seq: the number nodes to consider.
- 2. num_dim: the number of dimensions (which is 3 in this sample).
- 3. lambda_3: the penalty weight for this constraint.
- Outputs:
- 1. List of term objects.
- '''
-
- terms = []
- for k in range(0,len_seq):
- offset = k*num_dim
- terms += [Term(c=1*lambda_3,indices=[0+offset,1+offset,2+offset])]
- return terms
-
-
-##### ----- Test the function and print output:
-term_000 = penalize_000(len_seq, num_dim, lambda_2)
-term_111 = penalize_111(len_seq, num_dim, lambda_3)
-print('Term dictionary penalty constraint for 000:', term_000, '\n')
-print('Term dictionary penalty constraint for 111:', term_111, '\n')
-print('Function for penality 000: ')
-print_function(term_000)
-print('Function for penality 111: ')
-print_function(term_111)
-We've finished defining the optimization function. The next step is to start looking at how we're going to submit it to the solvers and analyze the returned results. -First let's look at how to parse, validate, and visualize the results, since that will make tuning the solvers easier!
-Below you can find the function definition that reads the solution and validates it. The solution dictionary is first read out to the optimiation variable string that describes the path. Afterward, the substrings that represent the directions are translated to linguistic terms, such that the solution can be printed in a human-readable format. Based on these two steps, the validation process checks if any constraints are violated. If constraints are violated warnings will be shown in the output with some tuning suggestions.
- -def read_validate_solution(solution: dict, num_turns: int, num_dim: int):
-
- '''
- Purpose:
- To validate the solution returned by the solver. Make it readable, and analyze if it makes sense.
- Inputs:
- 1. solution: The solution results dictionary which is returned by the solver (results["configuration"]).
- 2. num_turns: The number of turns for the simulation.
- 3. num_dim: The number of dimensions, which is 3 for this sample (3D).
- Outputs:
- 1. valid: A boolean variable that specifies the validity of the solution.
- 2. pos_dit: Layered position dictionary that contains all of the nodes' locations per turn {turn: {x: x_pos, y:y_pos, z:z_pos}}.
- 3. dir_dict: Dictionary containing the linguistic interpretation of the 6 directions.
- 4. var_dict: Dictionary containing the spin per optimization variable.
- 5. x_arr: Array of x positions.
- 6. y_arr: Array of y positions.
- 7. z_arr: Array of z positions.
- '''
-
- print('\n')
- valid = True
- move = ''
- sol_str = ''
- x_arr = [0]
- y_arr = [0]
- z_arr = [0]
- dir_dict = {'100':'out','010':'in','001':'right','110':'left','101':'up','011':'down'}
- pos_dict = {0:{"x":0, "y":0, "z":0}}
- var_dict = {}
- for key,val in solution:
- if key<(num_turns*num_dim):
- turn = floor(key/num_dim)+1
- print("Turn: "+str(turn),"var: "+str(key),"spin: "+str(val))
- var_dict |= {str(key): val}
- if key%3<num_dim and key<(num_turns*num_dim):
- move = move+str(val)
- if move in dir_dict:
- x_pos = pos_dict[turn-1]["x"]
- y_pos = pos_dict[turn-1]["y"]
- z_pos = pos_dict[turn-1]["z"]
- if dir_dict[move] == "out":
- x_pos += 1
- elif dir_dict[move] == "in":
- x_pos += -1
- elif dir_dict[move] == "right":
- y_pos += 1
- elif dir_dict[move] == "left":
- y_pos += -1
- elif dir_dict[move] == "up":
- z_pos += 1
- elif dir_dict[move] == "down":
- z_pos += -1
- new_pos = {turn: {"x":x_pos, "y":y_pos, "z":z_pos}}
- pos_dict |= new_pos
- x_arr += [x_pos]
- y_arr += [y_pos]
- z_arr += [z_pos]
- #print move
- print(dir_dict[move], '\n')
- sol_str = sol_str + " "+ dir_dict[move]
- move = ''
-
- elif move == ('111' or '000'):
- valid = False
- print('Illegal move')
-
- print("solution:", sol_str, '\n')
- print("positions: ", pos_dict)
-
- incorrect_pos = []
- for ref in range(0,len(pos_dict)):
- for tar in range(ref+1,len(pos_dict)):
- if pos_dict[ref] == pos_dict[tar]:
- print(f"Invalid position. Position already taken. Positions {ref}-{tar}, (if difference = 2, increase lambda_4, otherwise increase lambda_1)")
- valid = False
- incorrect_pos += [ref,tar]
-
- print('\n Incorrect positions',incorrect_pos)
- return valid, pos_dict, dir_dict, var_dict, x_arr, y_arr, z_arr
-Below you can find the function that plots the self-avoiding path. Probably this is the easist way to verify whether the path is correct or not. The initial position is plotted with a red "+" and the positions with a green dot, inteconnected by the lines.
- -def plot_path(x_arr,y_arr,z_arr):
-
- '''
- Purpose:
- This function plots the path of the self-avoiding walk through the position arrays.
- Inputs:
- 1. x_arr: the array of x positions.
- 2. y_arr: the array of y positions.
- 3. z_arr: the array of z positions.
- Outputs:
- 1. The 3D plot of the self-avoiding path.
- '''
- print('\n')
- print('x_arr: ', x_arr)
- print('y_arr: ', y_arr)
- print('z_arr: ', z_arr)
-
- fig = plt.figure()
- ax = plt.axes(projection ='3d')
- # plotting
- ax.plot3D(x_arr, y_arr, z_arr, 'g-o')
- ax.set_title('3D self-avoiding path')
- plt.plot(0,0,0,'r+') # plot initial position
- plt.show()
-Finally, time to submit to the Azure solvers! For this notebook a small problem is chosen because the execution times grow exponentially with the number of positions (path length). A self-avoiding path of 10 positions needs to be found.
-First it is necessary to build the entire optimization problem, comprising of the constraint functions. The optimization problem is then sent to a specified Azure QIO solver, such as simulated annealing or parallel tempering (heuristic solvers), which try to find a good solution to the problem. The the nonlinear optimization constraint weights and the solvers need tuning to find a valid path, which is a difficult process. To see a valid path, you can first run the with the provided parameters, as the solver and the cost function have been tuned to find these values. Other than that, feel free to play around with the values. Note that it is easier to tune for smaller problem sizes (nodes/positions) since the simulation times are smaller, however, the parameters do not map well to larger problems. Additionally, a suggestion is to first use the parameter-free solvers to limit the tuning work to the constraint weights.
-At the end of the notebook, a few resources are given where you can find more info about the QIO solvers.
-Below I define a function 'submit' which construct the optimization function, submits to a solver, and parses the results. If you want to re-run the optimization, you might first need to close the figure!
- -def submit():
-
- '''
- Purpose:
- This function submits to the Azure solvers.
- Adjust the variables and solver properties here.
- This functions coordinates the entire simulation by calling the previously defined functions.
- '''
-
- nodes = "HHHHHHHHHHHHHHH" # The 15 nodes that will be analyzed.
- len_seq = len(nodes) # The length of the sequence.
- num_turns = len(nodes)-1 # The number of turns between the nodes.
- num_dim = 3 # Number of optimization variables required to describe a turn
- lambda_0 = 1
- lambda_2 = 100
- lambda_3 = 100
- lambda_4 = 3
-
- terms = ( distance_constraint(num_turns, num_dim, lambda_0)+penalize_000(len_seq, num_dim, lambda_2)+
- penalize_111(len_seq, num_dim, lambda_3)+no_return_constraint(num_turns, num_dim, lambda_4) )
-
- problem = Problem(name="Self Avoiding Walk Problem", problem_type=ProblemType.pubo, terms=terms)
-
- # Submit
- start = time.time()
-
- ##### Two parameter-free solvers - these are tuned for you after submission.
- #parameter free simulated annealing
- #solver = SimulatedAnnealing(workspace, timeout=600)
- #parameer free parallel tempering
- #solver = ParallelTempering(workspace, timeout=600)
-
- ##### Parametrized simulated annealing - needs to be tuned by you.
- solver = SimulatedAnnealing(workspace, sweeps=50, restarts=2000, beta_start=1e-8, beta_stop=10, timeout=3600)
-
-
- print('Submitting problem...')
- job = solver.submit(problem)
-
- while job.details.status != 'Succeeded' and job.details.status != 'Failed':
- job.refresh()
- print('waiting')
- time.sleep(10)
-
- # Results
- results = job.get_results()
- print('\n Results:', results)
- config = results["configuration"]
-
- time_elapsed = time.time() - start
- print("Execution time in seconds: ", time_elapsed, '\n')
-
- solution = [(int(k),v) for k, v in config.items()]
- solution.sort(key=lambda tup: tup[0])
-
- valid, posDict, dirDict, configDict, x_arr, y_arr, z_arr = read_validate_solution(solution, num_turns, num_dim)
-
- plot_path(x_arr, y_arr, z_arr)
-
-submit()
-In this notebook we'll cover how to formulate an optimization problem to find a path through a 3D lattice that does not cross itself (self-avoiding). The optimization problem is solved with the Azure Quantum service which contains numerous heuristic (non-linear) optimization solvers. Finding a self-avoiding path in a 3D lattice can be considered a difficult problem, in the sense that it has suffers from exponential scaling conditioned on the number of turns and the number of dimensions (2D vs 3D, etc.).
+In this notebook, a step by step approach is taken to explain the function definitions to are necessary to create the optimization problem.
+Goal: To find a path in a 3D lattice that does not cross itself.
+To clarify some vocab and assumptions used in the notebook:
+Note: Because the number of terms to for this optimization problem grows explosively with the number of turns, you might want to try running this in the Azure Quantum notebooks. The more technical reason for this is that expansions of nonlinear terms have to be computed locally resulting in massive number of terms needing to be uploaded. Unless you are very patient, a recommendation would be to either check out the online notebook experience or rewrite this notebook's problem class to its streaming counterpart (https://docs.microsoft.com/en-us/azure/quantum/optimization-streaming-problem).
+ +import time
+from math import floor, log2
+from azure.quantum import Workspace
+from azure.quantum.optimization import Problem, ProblemType, Term, ParallelTempering, Tabu, SimulatedAnnealing
+from azure.identity import ClientSecretCredential
+from mpl_toolkits import mplot3d
+import numpy as np
+import matplotlib.pyplot as plt
+To run the sample, you'll need to have a quantum workspace. Check out this module if you don't have one yet: https://docs.microsoft.com/en-us/learn/modules/get-started-azure-quantum/. +Fill in the variables below.
+ +workspace = Workspace(
+ subscription_id = "",
+ resource_group = "",
+ name = "",
+ location = "",
+ credential = ClientSecretCredential(tenant_id="",
+ client_id="",
+ client_secret="")
+)
+Throughout the notebook you'll probably realize that printing the "terms" isn't going to provide much help, mainly because there will be too many of them and their dictionary format. A term is the way a mathematical term is expressed in terms of a dictionary for the SDK (see examples). Below a function is defined that prints the term(s) dictionaries that describe the constraint/cost function in terms of "q" optimization variables. You can call the function to check if the working is correct, or if you're unsure of what a constraint looks like. Note that for a large problem it nevertheless becomes difficult to understand the entire function output because of the number of terms.
+An example: +{'c': 1, 'ids': [0, 10, 11, 12]} ---> $1q_0q_{10}q_{11}q_{12}$
+ +def print_function(terms: list):
+
+ '''
+ Purpose:
+ Takes a list of terms and prints it out as a mathematical (cost/constraint) function.
+ Example:
+ {'c': 1, 'ids': [0, 10, 11, 12]} ---> 1q_0q_10q_11q_12
+ Inputs:
+ 1. terms: the list of terms.
+ '''
+
+ k = 0
+ final_string = ''
+ final_string = ''
+ for term in terms:
+ term = term.to_dict()
+ weight = term['c']
+ ids = term['ids']
+ string = '('
+ if weight >= 0:
+ if k == 0:
+ string = '('+str(weight)
+ else:
+ string = '+' + '('+str(weight)
+ if weight < 0:
+ if k == 0:
+ string = '(' + str(abs(weight))
+ else:
+ string = '-' + '(' + str(abs(weight))
+ for id_ in ids:
+ string = string + f'q_{id_}'
+ string = string + ')'
+ final_string = final_string + string
+ k += 1
+
+ print('[ ' + final_string + ' ]')
+Throughout the notebook we will need to cross-multiply nonlinear terms, for computing squares for example. The function below performs the expansion. An important note to keep in mind here, especially if you're looking to optimize the code, is that the function does not assemble identical terms. Since the cost functions dealt with in this notebook are still rather small, reducing the number of terms was left out of the scope.
+More clearly stated in terms of equations:
+$$ \text{Current function (see below): } \hspace{0.2cm} (a+b)^2 = a^2 + ab + ab + b^2 $$$$ \text{Optimal function: } \hspace{0.2cm} (a+b)^2 = a^2 + 2ab + b^2 $$Should be a fun programming exercise to reduce the number of terms significantly!
+ +def cross_multiply(list_a: list, list_b: list) -> list:
+
+ '''
+ Purpose: Cross multiplies two lists of terms (linear ex. [q_0 + q_1] or non-linear ex. [q_0q_1]) to return the expansion.
+ Can compute powers of groups this way (^2, ^3,...), like squaring a list of terms.
+ Calculates the expansion locally, unlike the SlcTerm class.
+ Example: (2q_0q_1+3q_2q_3)^2 => (2q_0q_1)^2 + 6q_0q_1q_2q_3 + 6q_0q_1q_2q_3 + (3q_2q_3)^2
+ Input:
+ 1. list_a: list which serves as the reference list (first 'for' loop).
+ 2. list_b: list which serves as the target list (second 'for' loop).
+ Output:
+ 1. list of term objects.
+ '''
+
+ terms = []
+ for one in list_a:
+ for uno in list_b:
+ alpha = one.to_dict()
+ beta = uno.to_dict()
+ weight = int(alpha['c'])*int(beta['c'])
+ ids = list(alpha['ids'])+list(beta['ids'])
+ terms += [Term(c = weight, indices = ids)]
+ return terms
+While running through the notebook you might want to execute and print some of the constraint functions. To do that, you'll need some variables to be defined. Only 3 turns are considered with the defined parameters below, as that will keep the function outputs small and understandable!
+ +nodes = "ABCD" # The four position/node names.
+len_seq = len(nodes) # The length of the node
+num_turns = len(nodes)-1 # The number of turns, which is one less than the number of therefore -1 turns.
+num_dim = 3 # Number of optimization variables required to describe a turn
+lambda_0 = 1 # Penalty weight for the 'distance constraint'
+lambda_1 = 40 # Penalty weight for the 'no return constraint'
+lambda_2 = 25 # Penalty weight for the invalid direction '000' constraint
+lambda_3 = 25 # Penalty weight for the invalid direction '111' constraint
+Prior to starting to define our cost function, we have to construct some basic functions and agree on some definitions. In this part, we'll define the 6 directions of the x-y-z plane in terms of the optimization variables ("q").
+A common way to encode a turn based optimization problem is through a direction-based string, in which each optimization variable substring (ex. $q_0q_1q_2$) denotes some kind of decision of which direction was takn. For example, take a car take can only move backward or forward. The decision of the driver at some time $t$ can be described as $q_t \in \{0, 1\}$, where 0 represents going backward and 1 forward, respectively. Over a multiple time steps, say 3, one can then describe the movement of the car by:
+$$ \text{Movements:} \hspace{0.25cm} q_0q_1q_2 $$The sequence of movements, "forward, backward, backward", can then be represented as:
+$$ \text{Sequence of movements:} \hspace{0.25cm} q_0q_1q_2 = 100$$Alright, if that is clear to you, then the next step is to represent the 6 directions of the x-y-z plane in such a fashion. Regrettably, a "q" can only be either 0 or 1, not 2, 3, 4, 5 to represent the six directions. But what can be done is to use multiple optimization variables to describe a direction! For example:
+$$ \text{Sequence of movements:} \hspace{0.25cm} q_0q_1q_2q_3q_4q_5q_6q_7q_8 $$$$ \text{First move:} \hspace{0.25cm} q_0q_1q_2 $$$$ \text{Second move:} \hspace{0.25cm} q_3q_4q_5 $$$$ \text{Third move:} \hspace{0.25cm} q_6q_7q_8 $$Because we need to encode 6 directions, we need at least 3 optimization variables to form unique substrings. The table explains this, we to have enough possible combinations to encode the 6 directions:
+Note, for later, you should keep in mind that 2 of the 8 combinations become irrelevant. It is important that the solver does not return those strings!
+Perfect! Now we just need to assign each direction a unique substring. Below you can find my choices, feel free to change these if you're starting from scratch.
+Now we just need to define direction variables. These are as defined in the above table. These variables should take a value 1 if a movement is made in that direction, while the other directions must all be zero. This can be achieved by the following scheme, where $k$ denotes the turn, and $\gamma = 3(k-1)$ (3 because 3 dimensions):
+To clarify a bit further through an example: +If in the first turn (k=0) the solver returns $q_0=1$, $q_1=0$, $q_2=0$, and in the second turn (k=1) returns $q_3=0$, $q_4=0$, $q_5=1$, then a "+x and +y" were taken, respecively.
+These direction variables will help understand some difficult constraints later on. Defining these direction variables additionally makes the code much more readable, since we can define everything in terms of the direction variables instead of the individual optimization variables.
+Below you can find the function for the direction variables. The formulas have been expanded (mathematically) to easily declare them in a number of terms, however it is hard-coded this way. Some inputs to the function ("sign_dir", "signpos", "lambda") are relevant for building constraints later in the notebook, but have to defined in this function. Further information can be found in the function docstring.
+ +def direction_variables(direction: str, offset: int, sign_dir: int, sign_pos: int, lambda_: int) -> list:
+
+ '''
+ Purpose:
+ Translates the direction (+x,-x,+y,-y,+z,-z) of turn 'i' as a function of three optimization variables. (Three q's because of the defined coordinate system).
+ Example:
+ Direction "+z" in the first turn (turn = 1) is translated to: q_{0+offset}q_{2+offset}-q_{0+offset}q_{1+offset}q_{2+offset}.
+ Inputs:
+ 1. direction: A direction from an x-y-z coordinate system, one of the following: ('+x','-x','+y','-y','+z','-z').
+ 2. offset: Offset gives the turn number expressed in the first "q" of that turn. Equal to gamma in the explanation.
+ Example: Turn 1 starts with "q" q_0, offset=0. Turn 2 starts with q_3, offset = 3.
+ 3. sign_dir: Changes the sign of the weights corresponding to negative directions "-x", "-y", "-z" -> necessary for finding the positions, for exmaple (+x) "-" (-x).
+ 4. sign_pos: Changes the sign of the weights corresponding to negative positions "-(x,y,z)" -> necessary for finding the distances between node "i" and node "j".
+ 5. lambda_: The weight term associated with a constraint.
+ Output:
+ 1. A list of term objects.
+ '''
+
+ terms = []
+ if direction == "+x":
+ term_0 = Term(c= 1*sign_pos*lambda_, indices=[0+offset])
+ term_1 = Term(c=-1*sign_pos*lambda_, indices=[0+offset, 1+offset])
+ term_2 = Term(c=-1*sign_pos*lambda_, indices=[0+offset, 2+offset])
+ term_3 = Term(c= 1*sign_pos*lambda_, indices=[0+offset, 1+offset, 2+offset])
+ terms = [term_0, term_1, term_2, term_3]
+ elif direction == "-x":
+ term_0 = Term(c= 1*sign_dir*sign_pos*lambda_, indices=[1+offset])
+ term_1 = Term(c=-1*sign_dir*sign_pos*lambda_, indices=[0+offset, 1+offset])
+ term_2 = Term(c=-1*sign_dir*sign_pos*lambda_, indices=[1+offset, 2+offset])
+ term_3 = Term(c= 1*sign_dir*sign_pos*lambda_, indices=[0+offset, 1+offset, 2+offset])
+ terms = [term_0, term_1, term_2, term_3]
+ elif direction == "+y":
+ term_0 = Term(c= 1*sign_pos*lambda_, indices=[2+offset])
+ term_1 = Term(c=-1*sign_pos*lambda_, indices=[0+offset, 2+offset])
+ term_2 = Term(c=-1*sign_pos*lambda_, indices=[1+offset, 2+offset])
+ term_3 = Term(c= 1*sign_pos*lambda_, indices=[0+offset, 1+offset, 2+offset])
+ terms = [term_0, term_1, term_2, term_3]
+ elif direction == "-y":
+ term_0 = Term(c= 1*sign_dir*sign_pos*lambda_, indices=[0+offset, 1+offset])
+ term_1 = Term(c=-1*sign_dir*sign_pos*lambda_, indices=[0+offset, 1+offset, 2+offset])
+ terms = [term_0, term_1]
+ elif direction == "+z":
+ term_0 = Term(c= 1*sign_pos*lambda_, indices=[0+offset, 2+offset])
+ term_1 = Term(c=-1*sign_pos*lambda_, indices=[0+offset, 1+offset, 2+offset])
+ terms = [term_0, term_1]
+ elif direction == "-z":
+ term_0 = Term(c= 1*sign_dir*sign_pos*lambda_, indices=[1+offset, 2+offset])
+ term_1 = Term(c=-1*sign_dir*sign_pos*lambda_, indices=[0+offset, 1+offset, 2+offset])
+ terms = [term_0, term_1]
+ return terms
+
+
+##### ----- Test the function and print output:
+turn = 1 # play with this value!
+dir_var = direction_variables("+x",(turn-1)*num_dim,1,1,lambda_0)
+print('dir_var term dictionary "+x": ', dir_var)
+print_function(dir_var)
+Great that we have the directions defined! Now let's use them and define a function that calculates the difference in positions, which we will need later. After each turn, a direction that appends to the path. A valid path is defined as a one that does not cross itself, meaning that we need to compare the positions over the turns. In other words, every new position appended to the path needs to be checked with all previously held positions, to verify that a position hasn't been visited twice.
+So how do we go about this?
+First we need a way to know the positions after choosing to go in a certain direction, so let's tackle that first. Consider the fact that the directions are recorded for each turn, which already gives some sort of log of the positions in the path. By summing the direction variables respective to their dimension (x/y/z) acoomplishes this, but only if a negative move can compensate a positive one. For example, after two turns the $x$-location can be expressed as:
+$$ x_2 = ( d^{1}_{+x} - d^{1}_{-x} ) + ( d^{2}_{+x} - d^{2}_{-x} )$$NOTE: In the above function ("direction_variables"), the negative sign for the negative direction variables is controlled through "sign_dir". "sign_dir" must be set to -1 to assign a negative sign to it as required in these position formulas!
+Remember that for a turn only one direction variable can be set to 1! Thus if two moves in the $-x$ direction are taken, then:
+$$ x_2 = ( d^{1}_{+x} - d^{1}_{-x} ) + ( d^{2}_{+x} - d^{2}_{-x} ) = ( 0 - 1 ) + ( 0 - 1 ) = -2 $$As with the simple example for the x-dimension, the same can be applied to the $x$ and $y$ directions. The formulas below descrive the position after some number of turns through a summation (for the maths/physics enthousiasts, these are integrals of the velocities to derive the positions :) ).
+$$ x^k = \sum_{k=1}^{k} ( d^{k}_{+x} - d^{k}_{-x} ) $$$$ y^k = \sum_{k=1}^{k} ( d^{k}_{+y} - d^{k}_{-y} ) $$$$ z^k = \sum_{k=1}^{k} ( d^{k}_{+z} - d^{k}_{-z} ) $$Great stuff! Now we have a method to find the positions along the path. To find the differences in position, for example between turn 2 and turn 1, it is only necessary to apply a subtraction:
+$$ x_2 - x_1 = \sum_{k=1}^{2} ( d^{k}_{+x} - d^{k}_{-x} ) - \sum_{k=1}^{1} ( d^{k}_{+x} - d^{k}_{-x}) = ( d^{2}_{+x} - d^{2}_{-x} ) $$$$ y_2 - y_1 = \sum_{k=1}^{2} ( d^{k}_{+y} - d^{k}_{-y} ) - \sum_{k=1}^{1} ( d^{k}_{+y} - d^{k}_{-y}) = ( d^{2}_{+y} - d^{2}_{-y} ) $$$$ z_2 - z_1 = \sum_{k=1}^{2} ( d^{k}_{+z} - d^{k}_{-z} ) - \sum_{k=1}^{1} ( d^{k}_{+z} - d^{k}_{-z}) = ( d^{2}_{+z} - d^{2}_{-z} ) $$Generalizing this to any difference between two positions:
+$$ x_j - x_i = \sum_{k=i}^{j} ( d^{k}_{+x} - d^{k}_{-x} ) $$$$ y_j - y_i = \sum_{k=i}^{j} ( d^{k}_{+y} - d^{k}_{-y} ) $$$$ z_j - z_i = \sum_{k=i}^{j} ( d^{k}_{+z} - d^{k}_{-z} ) $$Fantastic! Now on to the function code for the difference in positions. Do read the function and its docstring to understand the implementation, because it reference the previously defined function to make use of direction variables.
+ +def diff_in_pos(start_turn: int, end_turn: int, num_dim: int, lamda_: int):
+
+ '''
+ Purpose:
+ Expresses the difference in position (x,y,z) between two turns as a function of the q's encoding of the directions.
+ In other words, expresses the difference in the position after the start_turn and position after the end_turn.
+ Example:
+ Difference between turn 0 (no turns yet, initial position = (0,0,0)) and turn 2:
+ x(2) = [ move(+x, turn 1) - move(-x, turn 1) ] + [ move(+x, turn 2) - move(-x, turn 2) ]
+ Note: only one 'move' per turn gets activated as they are represented by the same q's (ex. turn 1 is represented by q_0, q_1, and q_2).
+ Inputs:
+ 1. start_turn: the initial (reference) turn.
+ 2. end_turn: the final (target) turn.
+ 3. num_dim: the number of dimensions (3, x-y-z coordinate system).
+ Outputs:
+ 1. The difference in the x direction.
+ 2. The difference in the y direction.
+ 3. The difference in the z direction.
+ '''
+
+ x_diff = y_diff = z_diff = []
+ if start_turn < end_turn and start_turn >= 0 and end_turn >= 1:
+ for turn in range(start_turn,end_turn+1):
+ x_diff += direction_variables("+x",(turn-1)*num_dim,1,1,lamda_)+direction_variables("-x",(turn-1)*num_dim,-1,1,lamda_)
+ y_diff += direction_variables("+y",(turn-1)*num_dim,1,1,lamda_)+direction_variables("-y",(turn-1)*num_dim,-1,1,lamda_)
+ z_diff += direction_variables("+z",(turn-1)*num_dim,1,1,lamda_)+direction_variables("-z",(turn-1)*num_dim,-1,1,lamda_)
+ return x_diff, y_diff, z_diff
+
+
+##### ----- Test the function and print output:
+start_turn = 1 # play with this value!
+end_turn = 2 # play with this value!
+x_diff, y_diff, z_diff = diff_in_pos(start_turn, end_turn, num_dim, lambda_0)
+print('Difference in x described as term list: \n', x_diff)
+print('\nDifference in x described as mathematical function:')
+print_function(x_diff)
+Fantastic that you've made it this far! With all the above functions defined we can start solving the problem.
+The distance constraint enforces that after each turn, the difference in positions after turns $i$ and $j$ must be larger or equal to 1 (distance between lattice points in the same dimension). +To visualize the scenario this constraint tries to penalize consider a 2D example with 4 turns, starting in $(0,0)$:
+$$ \text{Turn 1, go right: } \hspace{0.3cm} (0,0) => (0,1) $$+$$ \text{Turn 2, go up: } \hspace{0.3cm} (0,1) => (1,1) $$ +$$ \text{Turn 3, go left: } \hspace{0.3cm} (1,1) => (1,0) $$ +$$ \text{Turn 4, go down: } \hspace{0.3cm} (1,0) => (0,0) $$
+As you can see, the last step conflicts with our constraint. It is not permitted to return into a position that we've already been, namely $(0,0)$. Turn 1-3 are all valid, since they remain a distance of 1 away from all other previously held positions.
+Let's define what we want mathematically. tThe distance between the position after turn $i$ and $j$ must be larger or equal to 1. For this we'll need some basic geometry, the 3D variant of the Pythagorean theorem. From the theorem, it can be understood that hypotenuse ($c$) must always be larger or equal to one, since that is what defines the distance between two points. Mathematically speaking working this out:
+$$ C_{i,j} \geq 1 $$$$ \sqrt{(x_{j} - x_{i})^2 + (y_{j} - y_{i})^2 + (z_{j} - z_{i})^2} \geq 1 $$The square-root is a bit of a headache when writing out a cost function, it would require an approximation or some sort of factorization, not fun if you're dealing with many terms. To avoid that, we're simply going to square both sides, which in our case if fine (for maths people), since we aren't dealing with negative numbers here:
+$$ C_{i,j}^2 \geq 1 $$$$ {(x_{j} - x_{i})^2 + (y_{j} - y_{i})^2 + (z_{j} - z_{i})^2} \geq 1 $$Nice! If you're sharp, you'll start seeing where the previous functions might come in handy now. The function that computes the difference between the respective positions ("diff_in_pos") in each dimension was defined above, but also the function that computes the an expansion ("cross_multiply") was defined at the beginning of the notebook. Writing those function calls out in order :
+Alright, you might be wondering how to deal with the $>=$, since it is necessary to have an equation that equals zero in order to integrate it into the cost function. +This is where slack variables come in. The purpose of slack variables is to convert an inequality constraint into an equality constraint. Or in other words, for this equation, we need to compensate some values (explained further below) in order to convert the equation to an equality constraint. Consider the constraint in its simplest form:
+$$ C_{i,j}^2 \geq 1 $$It is clear that we would need to add some number $V_{i,j}$, on the right hand side of the equation in order to make this an equality constraint. $V_{i,j}$ can be any positive value and is conditioned on $i$ and $j$:
+$$ C_{i,j}^2 = 1+V_{i,j}$$Let's say that instead of representing a single value we want V_{i,j} to represent a range of values $[0-4]$. This is where slack variables can be introduced. By introducing additional optimization variables, which have nothing to do with the path encoding, the range can be described as following:
+$$ V_{i,j} = [0-4] = q_{s1}+q_{s2}+q_{s3}+q_{s4} $$Or more compactly:
+$$ V_{i,j} = [0-4] = q_{s1}+q_{s2}+2q_{s3} $$So any value in the range $[0-4]$ is described through these q's. You can derive 1 by assigning a single "q", which has weight 1, the value 1. Likewise, you can get the value 4 by assigning all optimization variables the value 1.
+Alright, so how do we know which range of values we need to represent? That is revealed by the turn numbers, $i$ and $j$. Over these turns, the maximum squared distance that can be achieved is when moves are only made in a single direction, the "+x" direction for example. The maximum distance is then equal to the number of moves squared $(j-i)^2$.
+$$ C_{i,j}^2 = (j-i)^2 $$Then writing out the maths gives us the necessary upper bound for the range that the optimization variables need to represent, where $S$ stands for the number of slack variables necessary:
+$$ \text{Upper bound: } \hspace{0.5cm} C_{i,j}^2 = (j-i)^2 = 1 + \sum_{s=0}^{S}q_s$$$$ \text{Upper bound: } \hspace{0.5cm} \sum_{s=0}^{S}q_s = (j-i)^2 - 1 $$For the lower bound of the range, the value is 0. The reason for that is we want to be able to represent any number equal to or larger than 1 for the inequality constraint through the addition of the slack variables. Concluding, the range the q's must represent is $[0,((j-i)^2-1)]$. More on how to calculate the number of slack variables necessary to define this range later, as you don't need ((j-i)^2-1) optimization variables!
+Piecing all the parts together we find the following equation for a $i$-$j$ combination:
+$$ {(x_{j} - x_{i})^2 + (y_{j} - y_{i})^2 + (z_{j} - z_{i})^2} = 1 + \sum_{s=0}^{S}q_s $$ +Because constants are irrelevant for the optimization landscape the value 1 can be neglected. The reason being that constants only introduce linear offsets, thus impacting the entire optimization equally. The distance constraint which needs to account for all $i$-$j$ combinations, where $j$>$i$ and $N-1$ the number of turns, is summarized as follows:
+$$ \sum_{i=1}^{N-1} \sum_{j>i}^{N-1} \lambda_0 ( {(x_{j} - x_{i})^2 + (y_{j} - y_{i})^2 + (z_{j} - z_{i})^2} - \sum_{s=0}^{S}q_s ) $$ +The function that builds this constraint is presented below. Read the docstring for how it works.
+Note, to run this function you'll need to run the next function cell to define 'generate_slack_coefficients'.
def distance_constraint(num_turns: int, num_dim: int, lambda_0: int) -> list:
+
+ '''
+ Purpose:
+ Build the distance contraint based on previosly defined functions.
+ Constraint: Distance squared between i and j must be larger or equal to 1.
+ Constraint: L_{i,j}^2 >= 1 => L_{i,j}^2 = 1+q_{slacker} (converting the inequality constraint to an equality constraint.)
+ Example/Explanation:
+ L{1,2}^2 = (x1-x2)^2 + (y1-y2)^2 + (z1-z2)^2
+ x1 = [move(+x1)-move(-x1)]
+ x2 = [move(+x1)-move(-x1)] + [move(+x2)-move(-x2)]
+ <same for other dimensions >
+ <fill into first line>
+ L{1,2}^2 = [move(+x2)-move(-x2)]^2 + [move(+y2)-move(-y2)]^2 + [move(+z2)-move(-z2)]^2
+ [move(+x2)-move(-x2)]^2 + [move(+y2)-move(-y2)]^2 + [move(+z2)-move(-z2)]^2 - q_{slacker} = 0 ---> expressed in q's == cost function
+ Inputs:
+ 1. num_turns: the number of turns.
+ 2. num_dum: the number of dimensions.
+ 3. lambda_0: the constraint weight for the distance constraint.
+ Output:
+ 1. List of term objects describing the distance constraint.
+ '''
+
+ terms = []
+ slack_indexer = 0
+ for start_turn in range(1,num_turns+1):
+ for end_turn in range(start_turn+1,num_turns+2):
+ # Calculate the differences in positions for each dimension.
+ x_diff_i_j, y_diff_i_j, z_diff_i_j = diff_in_pos(start_turn, end_turn, num_dim, lambda_0)
+ # Compute the squared distance (Pythagorean theorem) by calculating the squared expansion.
+ x_diff_i_j_2, y_diff_i_j_2, z_diff_i_j_2 = cross_multiply(x_diff_i_j,x_diff_i_j), cross_multiply(y_diff_i_j,y_diff_i_j), cross_multiply(z_diff_i_j,z_diff_i_j)
+ # Add slack variables due to inequality constraint.
+ slack_var_terms = []
+ slack_coefficients = generate_slack_coefficients(end_turn-start_turn)
+ for s in range(0,len(slack_coefficients)):
+ slack_var_terms += [Term(c=-slack_coefficients[s], indices=[num_turns*num_dim+slack_indexer+s])]
+ terms += x_diff_i_j_2 + y_diff_i_j_2 + z_diff_i_j_2 + slack_var_terms
+ slack_indexer+=len(slack_coefficients)
+ return terms
+
+##### ----- Test the function and print output:
+# The output is too large to understand, nevertheless you can run the below statements to view how visualize the term-scaling of the problem.
+number_turns = 2 #play with this value!
+dist_terms = distance_constraint(number_turns, num_dim, lambda_0)
+
+# Only print below statements if you want to get an impression of the number of terms.
+#print('Distance constraint term dictionaries: ', dist_terms) # not readable!
+#print('\nDistance constraint: ')
+#print_function(dist_terms) # not readable!
+To remove any human calculations to construct the distance constraint(s), we'll automate the computation of the number of slack variables and their weights.
+As you saw in the previous part, certain slack variables can be assigned weights to reduce the number of slack variables necessary. Recall that the range $[0,4]$ can be defined by the following two equations with slack variables (s):
+$$ V_{i,j} = [0-4] = q_{s1}+q_{s2}+q_{s3}+q_{s4} $$$$ V_{i,j} = [0-4] = q_{s1}+q_{s2}+2q_{s3} $$The bottom equation is more compact because the last variables carries a weight-2. In Lucas' paper (see 2.4 of arxiv:1302.5843) the algorithm to find such weights is presented, which is based on calculating all $2^{x}$, where $x$ iterates from zero to the log_2 of maximum value in the range (4 in this example).
+Below are two examples that may help visualize why and how the slack variables represent the ranges of values. +\ +
+\ +
+\ +
+The function is given below. In context of the distance constraint, we're still faced with a minor problem. If $(j-i)^2-1$ equals 0, then we don't need any slack variables. The reason for this is that the range becomes $[0,0]$, which can also be represented by no variables at all.
+ +def generate_slack_coefficients(turn_diff: int):
+
+ '''
+ Purpose:
+ Calculates the number of slack variables and their weights.
+ Example:
+ For the constraint: x1 + x2 + x3 + x4 <= 4 which is converted to x1 + x2 + x3 + x4 + s1 + s2 +2*s3 = 4,
+ with 's' being slack variables. This function computes the weights of these slack variables ([1,1,2] for the example).
+ Reference:
+ Lucas' paper (see 2.4 of arxiv:1302.5843)
+ Input:
+ 1. turn_diff: the differences in turns (end_turn - start_turn).
+ Output:
+ 1. y: the weights of the slack variables.
+ '''
+
+ dist_diff = (turn_diff**2)-1
+ if dist_diff == 0:
+ y = [] # no slack variables needed
+ elif dist_diff > 0:
+ M = floor(log2(dist_diff))
+ y = [-2**n for n in range(M)]
+ y.append(-(dist_diff + 1 - 2**M))
+ return y
+
+##### ----- Test the function and print output:
+turn_difference = 3 # play with this value!
+slack_weights = generate_slack_coefficients(turn_difference)
+print('Slack weights: ', slack_weights)
+The distance constraint enforces a minimum distance of 1 between the different taken positions. However, for two consecutive turns, it does not completely prevent (to a high enough degree) going to a previous position. A simplified explanation for this is that the distance constraint is described through the positional differences in $x$, $y$, and $z$. Additionally, for consecutive turns, the number of slack variables is zero ($(i-j)^2-1 = 0$ for consecutive turns} - thus there is no auxiliary benefit introduced by slack variables to compensate for being in a different position after the turn. Therefore, the solver will try to keep the positions as close to each other as possible, which leads to overlapping positions, especially in consecutive turns. To overcome this, a penalization needs to be introduced that stops returning to the previous position held in the turn before. For example, what we want to stop is the following:
+$$ \text{Turn 1: A to B.} $$$$ \text{Turn 2: B to A.} $$This can be equivalenty described as moving in opposite directions over the consecutive turns. For example, first going in the "+y" direction, and following that moving in the "-y" direction.
+Designing such a constraint is straightforward, the opposite direction variables need to be multiplied. +Preventing a revisit of a position for the "x" direction over the first and second turn can be described by the following constraint:
+$$ {\lambda}_2(d^{1}_{+x} d^{2}_{-x} + d^{1}_{-x}d^{2}_{+x}) $$ +If a move is made in the "+x" direction first, $d^{1}_{+x}$ takes value 1. If afterwards a move is made in the "-x" direction, $d^{2}_{-x}$ also takes value 1, meaning that the combined term +$d^{1}_{+x} d^{2}_{-x}$ also becomes 1, enforcing the constraint with penalty value ${\lambda}_2$.
+Generalizing the idea to any direction and for all turns, we find the "no return constraint":
+$$ {\lambda}_1 \left( \sum_{m}^{\in \{x,y,z\} }\sum_{t=1}^{N-1} d^{t}_{+m} d^{t+1}_{-m} + d^{t}_{-m}d^{t+1}_{+m} \right) $$ +In this constraint, iterations are performed over the different dimensions ($m$), and the turn ($t$). The function defintion for the no return constraint is given below. Make sure to read the docstring. In the function the 'cross_multiply' method is used to expand and calculate all the necessary terms of the constraint. The direction variables are dictionaries of polynomials, and therefore have to be expanded before submitting to the Azure QIO solvers.
+ +def no_return_constraint(num_turns: int, num_dim: int, lamda_4: int)-> list:
+
+ '''
+ Purpose:
+ Build the constraint that penalizes going back to the same position/node two turns later.
+ Example:
+ Node A ---> <move +x> ---> Node B ---> <move -x> ---> Node A => erroneous as we've been there already.
+ Two sequential moves may not be in the same dimension and in opposite directions: (+x then -x), (-x then +x), (+y then -y) etc.
+ Inputs:
+ 1. num_turns: the number of turns.
+ 2. num_dim : the number of dimensions, which is 3 for this sample.
+ 3. lambda_4 : the penalty weight for this constraint.
+ Outputs:
+ 1. List of term objects.
+ '''
+
+ terms = []
+ for i in range(0,num_turns):
+ x_out_in = cross_multiply(direction_variables("+x",i*num_dim,1,1,lamda_4), direction_variables("-x",(i+1)*num_dim,1,1,lamda_4))
+ x_in_out = cross_multiply(direction_variables("-x",i*num_dim,1,1,lamda_4), direction_variables("+x",(i+1)*num_dim,1,1,lamda_4))
+ y_right_left = cross_multiply(direction_variables("+y",i*num_dim,1,1,lamda_4), direction_variables("-y",(i+1)*num_dim,1,1,lamda_4))
+ y_left_right = cross_multiply(direction_variables("-y",i*num_dim,1,1,lamda_4), direction_variables("+y",(i+1)*num_dim,1,1,lamda_4))
+ z_up_down = cross_multiply(direction_variables("+z",i*num_dim,1,1,lamda_4), direction_variables("-z",(i+1)*num_dim,1,1,lamda_4))
+ z_down_up = cross_multiply(direction_variables("-z",i*num_dim,1,1,lamda_4), direction_variables("+z",(i+1)*num_dim,1,1,lamda_4))
+ terms += x_out_in + x_in_out + y_right_left + y_left_right + z_up_down + z_down_up
+ return terms
+
+
+##### ----- Test the function and print output:
+number_turns = 1 # play with this value!
+nrc_terms = no_return_constraint(number_turns, num_dim, lambda_1)
+print('No return constraint term dictionaries: ', nrc_terms)
+print('\n No return constraint function:')
+print_function(nrc_terms)
+Recall that there are two invalid optimization variable substrings that are not associated with any move, $000$ and $111$. A path that includes these substrings is invalid. To prevent the solver from generating these invalid directions we'll need to add two constraints. Each constraint is specific to an invalid direction, and because we're working directly with the substring we don't need the abstraction layers used in previous constraints, like directional and positional variables.
+Let's first work out the constraint for preventing substring $000$. $ q_{0+\gamma} q_{1+\gamma} q_{2+\gamma} $ may not equal $000$, where $\gamma = 3(k-1)$ and $k$ the turn number. If we were to design the constraints as $ q_{0+\gamma} q_{1+\gamma} q_{2+\gamma} $ there won't be any penalty, because all $q$'s will equal and multiply to zero. Therefore, if all these $q$'s equal zero, we need the substring to multiply to 1. This is achieved by the following, which multiplies to zero if all optimization variables take the value zero:
+$$ {\lambda}_2(1-q_{0+\gamma}) (1-q_{1+\gamma}) (1-q_{2+\gamma}) $$$$ \text{If } \hspace{0.1cm} q_{0+\gamma} = q_{1+\gamma} = q_{2+\gamma} = 0, \hspace{0.1cm} \text{then: } $$+$$ {\lambda}_2(1-q_{0+\gamma}) (1-q_{1+\gamma}) (1-q_{2+\gamma}) = {\lambda}_2 \cdot 1 $$
+Nice! By expanding the equation and neglecting constant terms a constraint is found that can be implemented:
+$$\text{Constraint for 000: }-q_{0+\gamma}−q_{1+\gamma}−q_{2+\gamma}+q_{0+\gamma}q_{1+\gamma}+q_{0+\gamma}q_{2+\gamma}+q_{1+\gamma}q_{2+\gamma}−q_{0+\gamma}q_{1+\gamma}q_{2+\gamma}$$Before going to the function definition, let's first look at penalizing the substring $111$. Luckily this is much easier as the substring multiplies to 1 if all the $q$'s have the value zero. The constraint is therefore simple to derive:
+$$ {\lambda}_3 (q_{0+\gamma} q_{1+\gamma} q_{2+\gamma}) $$$$ \text{If } \hspace{0.1cm} q_{0+\gamma} = q_{1+\gamma} = q_{2+\gamma} = 1, \hspace{0.1cm} \text{then: } $$+$$ {\lambda}_3 (q_{0+\gamma} q_{1+\gamma} q_{2+\gamma}) = {\lambda}_3 \cdot 1$$
+The constraint for the substring $111$ is:
+$$ \text{Constraint for $111$: }\hspace{0.1cm} {\lambda}_3 (q_{0+\gamma} q_{1+\gamma} q_{2+\gamma}) $$The function definitions for these two constraints are given below.
+ +def penalize_000(len_seq: int, num_dim: int, lambda_2: int) -> list:
+
+ '''
+ Purpose:
+ Build the constraint that penalizes the invalid moves associated with the 3 q's string: '000'.
+ The constraint needs to be defined over the set of all turns, as defined by the 'for' loop.
+ Example:
+ If turn 2 (q_3q_4q_5) equals '000', assign a large penalty.
+ Inputs:
+ 1. len_seq: the number nodes to consider.
+ 2. num_dim: the number of dimensions (which is 3).
+ 3. lambda_2: the penalty weight for this constraint.
+ Outputs:
+ 1. List of term objects.
+ '''
+
+ terms = []
+ for k in range(0,len_seq):
+ offset = k*num_dim
+ term_0 = Term(c=-1*lambda_2,indices=[0+offset])
+ term_1 = Term(c=-1*lambda_2,indices=[1+offset])
+ term_2 = Term(c=-1*lambda_2,indices=[2+offset])
+ term_3 = Term(c= 1*lambda_2,indices=[0+offset,1+offset])
+ term_4 = Term(c= 1*lambda_2,indices=[0+offset,2+offset])
+ term_5 = Term(c= 1*lambda_2,indices=[1+offset,2+offset])
+ term_6 = Term(c=-1*lambda_2,indices=[0+offset,1+offset,2+offset])
+ terms += [term_0, term_1, term_2, term_3, term_4, term_5, term_6]
+ return terms
+
+def penalize_111(len_seq: int, num_dim: int, lambda_3: int) -> list:
+
+ '''
+ Purpose:
+ Build the constraint that penalizes the invalid moves associated with the 3 q's string: '111'.
+ The constraint needs to be defined over the set of all turns, as defined by the 'for' loop.
+ Example:
+ If turn 2 (q_3q_4q_5) equals '111', assign a large penalty.
+ Inputs:
+ 1. len_seq: the number nodes to consider.
+ 2. num_dim: the number of dimensions (which is 3 in this sample).
+ 3. lambda_3: the penalty weight for this constraint.
+ Outputs:
+ 1. List of term objects.
+ '''
+
+ terms = []
+ for k in range(0,len_seq):
+ offset = k*num_dim
+ terms += [Term(c=1*lambda_3,indices=[0+offset,1+offset,2+offset])]
+ return terms
+
+
+##### ----- Test the function and print output:
+term_000 = penalize_000(len_seq, num_dim, lambda_2)
+term_111 = penalize_111(len_seq, num_dim, lambda_3)
+print('Term dictionary penalty constraint for 000:', term_000, '\n')
+print('Term dictionary penalty constraint for 111:', term_111, '\n')
+print('Function for penality 000: ')
+print_function(term_000)
+print('Function for penality 111: ')
+print_function(term_111)
+We've finished defining the optimization function. The next step is to start looking at how we're going to submit it to the solvers and analyze the returned results. +First let's look at how to parse, validate, and visualize the results, since that will make tuning the solvers easier!
+Below you can find the function definition that reads the solution and validates it. The solution dictionary is first read out to the optimiation variable string that describes the path. Afterward, the substrings that represent the directions are translated to linguistic terms, such that the solution can be printed in a human-readable format. Based on these two steps, the validation process checks if any constraints are violated. If constraints are violated warnings will be shown in the output with some tuning suggestions.
+ +def read_validate_solution(solution: dict, num_turns: int, num_dim: int):
+
+ '''
+ Purpose:
+ To validate the solution returned by the solver. Make it readable, and analyze if it makes sense.
+ Inputs:
+ 1. solution: The solution results dictionary which is returned by the solver (results["configuration"]).
+ 2. num_turns: The number of turns for the simulation.
+ 3. num_dim: The number of dimensions, which is 3 for this sample (3D).
+ Outputs:
+ 1. valid: A boolean variable that specifies the validity of the solution.
+ 2. pos_dit: Layered position dictionary that contains all of the nodes' locations per turn {turn: {x: x_pos, y:y_pos, z:z_pos}}.
+ 3. dir_dict: Dictionary containing the linguistic interpretation of the 6 directions.
+ 4. var_dict: Dictionary containing the spin per optimization variable.
+ 5. x_arr: Array of x positions.
+ 6. y_arr: Array of y positions.
+ 7. z_arr: Array of z positions.
+ '''
+
+ print('\n')
+ valid = True
+ move = ''
+ sol_str = ''
+ x_arr = [0]
+ y_arr = [0]
+ z_arr = [0]
+ dir_dict = {'100':'out','010':'in','001':'right','110':'left','101':'up','011':'down'}
+ pos_dict = {0:{"x":0, "y":0, "z":0}}
+ var_dict = {}
+ for key,val in solution:
+ if key<(num_turns*num_dim):
+ turn = floor(key/num_dim)+1
+ print("Turn: "+str(turn),"var: "+str(key),"spin: "+str(val))
+ var_dict |= {str(key): val}
+ if key%3<num_dim and key<(num_turns*num_dim):
+ move = move+str(val)
+ if move in dir_dict:
+ x_pos = pos_dict[turn-1]["x"]
+ y_pos = pos_dict[turn-1]["y"]
+ z_pos = pos_dict[turn-1]["z"]
+ if dir_dict[move] == "out":
+ x_pos += 1
+ elif dir_dict[move] == "in":
+ x_pos += -1
+ elif dir_dict[move] == "right":
+ y_pos += 1
+ elif dir_dict[move] == "left":
+ y_pos += -1
+ elif dir_dict[move] == "up":
+ z_pos += 1
+ elif dir_dict[move] == "down":
+ z_pos += -1
+ new_pos = {turn: {"x":x_pos, "y":y_pos, "z":z_pos}}
+ pos_dict |= new_pos
+ x_arr += [x_pos]
+ y_arr += [y_pos]
+ z_arr += [z_pos]
+ #print move
+ print(dir_dict[move], '\n')
+ sol_str = sol_str + " "+ dir_dict[move]
+ move = ''
+
+ elif move == ('111' or '000'):
+ valid = False
+ print('Illegal move')
+
+ print("solution:", sol_str, '\n')
+ print("positions: ", pos_dict)
+
+ incorrect_pos = []
+ for ref in range(0,len(pos_dict)):
+ for tar in range(ref+1,len(pos_dict)):
+ if pos_dict[ref] == pos_dict[tar]:
+ print(f"Invalid position. Position already taken. Positions {ref}-{tar}, (if difference = 2, increase lambda_4, otherwise increase lambda_1)")
+ valid = False
+ incorrect_pos += [ref,tar]
+
+ print('\n Incorrect positions',incorrect_pos)
+ return valid, pos_dict, dir_dict, var_dict, x_arr, y_arr, z_arr
+Below you can find the function that plots the self-avoiding path. Probably this is the easist way to verify whether the path is correct or not. The initial position is plotted with a red "+" and the positions with a green dot, inteconnected by the lines.
+ +def plot_path(x_arr,y_arr,z_arr):
+
+ '''
+ Purpose:
+ This function plots the path of the self-avoiding walk through the position arrays.
+ Inputs:
+ 1. x_arr: the array of x positions.
+ 2. y_arr: the array of y positions.
+ 3. z_arr: the array of z positions.
+ Outputs:
+ 1. The 3D plot of the self-avoiding path.
+ '''
+ print('\n')
+ print('x_arr: ', x_arr)
+ print('y_arr: ', y_arr)
+ print('z_arr: ', z_arr)
+
+ fig = plt.figure()
+ ax = plt.axes(projection ='3d')
+ # plotting
+ ax.plot3D(x_arr, y_arr, z_arr, 'g-o')
+ ax.set_title('3D self-avoiding path')
+ plt.plot(0,0,0,'r+') # plot initial position
+ plt.show()
+Finally, time to submit to the Azure solvers! For this notebook a small problem is chosen because the execution times grow exponentially with the number of positions (path length). A self-avoiding path of 10 positions needs to be found.
+First it is necessary to build the entire optimization problem, comprising of the constraint functions. The optimization problem is then sent to a specified Azure QIO solver, such as simulated annealing or parallel tempering (heuristic solvers), which try to find a good solution to the problem. The the nonlinear optimization constraint weights and the solvers need tuning to find a valid path, which is a difficult process. To see a valid path, you can first run the with the provided parameters, as the solver and the cost function have been tuned to find these values. Other than that, feel free to play around with the values. Note that it is easier to tune for smaller problem sizes (nodes/positions) since the simulation times are smaller, however, the parameters do not map well to larger problems. Additionally, a suggestion is to first use the parameter-free solvers to limit the tuning work to the constraint weights.
+At the end of the notebook, a few resources are given where you can find more info about the QIO solvers.
+Below I define a function 'submit' which construct the optimization function, submits to a solver, and parses the results. If you want to re-run the optimization, you might first need to close the figure!
+ +def submit():
+
+ '''
+ Purpose:
+ This function submits to the Azure solvers.
+ Adjust the variables and solver properties here.
+ This functions coordinates the entire simulation by calling the previously defined functions.
+ '''
+
+ nodes = "HHHHHHHHHHHHHHH" # The 15 nodes that will be analyzed.
+ len_seq = len(nodes) # The length of the sequence.
+ num_turns = len(nodes)-1 # The number of turns between the nodes.
+ num_dim = 3 # Number of optimization variables required to describe a turn
+ lambda_0 = 1
+ lambda_2 = 100
+ lambda_3 = 100
+ lambda_4 = 3
+
+ terms = ( distance_constraint(num_turns, num_dim, lambda_0)+penalize_000(len_seq, num_dim, lambda_2)+
+ penalize_111(len_seq, num_dim, lambda_3)+no_return_constraint(num_turns, num_dim, lambda_4) )
+
+ problem = Problem(name="Self Avoiding Walk Problem", problem_type=ProblemType.pubo, terms=terms)
+
+ # Submit
+ start = time.time()
+
+ ##### Two parameter-free solvers - these are tuned for you after submission.
+ #parameter free simulated annealing
+ #solver = SimulatedAnnealing(workspace, timeout=600)
+ #parameer free parallel tempering
+ #solver = ParallelTempering(workspace, timeout=600)
+
+ ##### Parametrized simulated annealing - needs to be tuned by you.
+ solver = SimulatedAnnealing(workspace, sweeps=50, restarts=2000, beta_start=1e-8, beta_stop=10, timeout=3600)
+
+
+ print('Submitting problem...')
+ job = solver.submit(problem)
+
+ while job.details.status != 'Succeeded' and job.details.status != 'Failed':
+ job.refresh()
+ print('waiting')
+ time.sleep(10)
+
+ # Results
+ results = job.get_results()
+ print('\n Results:', results)
+ config = results["configuration"]
+
+ time_elapsed = time.time() - start
+ print("Execution time in seconds: ", time_elapsed, '\n')
+
+ solution = [(int(k),v) for k, v in config.items()]
+ solution.sort(key=lambda tup: tup[0])
+
+ valid, posDict, dirDict, configDict, x_arr, y_arr, z_arr = read_validate_solution(solution, num_turns, num_dim)
+
+ plot_path(x_arr, y_arr, z_arr)
+
+submit()
+