corrected symmetry cosntraint

backend/src/optimizer.py

@@ -10,11 +10,44 @@ class landmark :
         self.attractiveness = attractiveness
         self.loc = loc
 
+
+
+
+def untangle2(resx: list) :
+    N = len(resx)                   # length of res
+    L = int(np.sqrt(N))             # number of landmarks. CAST INTO INT but should not be a problem because N = L**2 by def.
+    n_edges = resx.sum()      # number of edges
+
+    order = []
+    nonzeroind = np.nonzero(resx)[0] # the return is a little funny so I use the [0]
+
+    nonzero_tup = np.unravel_index(nonzeroind, (L,L))
+
+    indx = nonzero_tup[0].tolist()
+    indy = nonzero_tup[1].tolist()
+
+    vert = (indx[0], indy[0])
+
+    order.append(vert[0])
+    order.append(vert[1])
+    
+    while len(order) < n_edges + 1 :
+        ind = indx.index(vert[1])
+
+        vert = (indx[ind], indy[ind])
+
+        order.append(vert[1])
+
+    return order
+
+
+
+
 # Convert the result (edges from j to k like d_25 = edge between vertex 2 and vertex 5) into the list of indices corresponding to the landmarks
-def untangle(res: list) :
+def untangle(resx: list) :
-    N = len(res)                # length of res
+    N = len(resx)                # length of res
     L = int(np.sqrt(N))         # number of landmarks. CAST INTO INT but should not be a problem because N = L**2 by def.
-    n_landmarks = res.sum()     # number of visited landmarks
+    n_landmarks = resx.sum()     # number of edges
     visit_order = []
     cnt = 0
 
@@ -40,47 +73,67 @@ def untangle(res: list) :
     return visit_order
 
 # Just to print the result
-def print_res(res: list, P) :
+def print_res(res: list, landmarks: list, P) :
     X = abs(res.x)
-    order = untangle(X)
+
+    N = int(np.sqrt(len(X)))
+    for i in range(N):
+        print(X[i*N:i*N+N])
+
+    order = untangle2(X)
+
+    order_ideal = [0, 0, 0, 0, 0, 0, 1, 0]
 
     # print("Optimal value:", -res.fun)  # Minimization, so we negate to get the maximum
     # print("Optimal point:", res.x)
-    # N = int(np.sqrt(len(X)))
+    
-    # for i in range(N):
+    #for i,x in enumerate(X) : X[i] = round(x,0)
-    #     print(X[i*N:i*N+N])
+    
-    # print(order)
+    #print(order)
 
     if (X.sum()+1)**2 == len(X) : 
         print('\nAll landmarks can be visited within max_steps, the following order is most likely not the fastest')
     else :
         print('Could not visit all the landmarks, the following order could be the fastest but not sure')
     print("Order of visit :")
-    for i, elem in enumerate(landmarks) : 
+    for idx in order : 
-        if i in order : print('- ' + elem.name)
+        print('- ' + landmarks[idx].name)
 
     steps = path_length(P, abs(res.x))
     print("\nSteps walked : " + str(steps))
 
-# Constraint to use only the upper triangular indices for travel
+# Constraint to not have d14 and d41 simultaneously
-def break_sym(landmarks, A_eq, b_eq):
+def break_sym(landmarks, A_ub, b_ub):
     L = len(landmarks)
-    l = [0]*L*L
+    upper_ind = np.triu_indices(L,0,L)
-    for i in range(L) :
-        for j in range(L) :
-            if i >= j :
-                l[j+i*L] = 1
 
-    A_eq = np.vstack((A_eq,l))
+    up_ind_x = upper_ind[0]
-    b_eq.append(0)
+    up_ind_y = upper_ind[1]
 
-    return A_eq, b_eq
+    for i, _ in enumerate(up_ind_x) :
+        l = [0]*L*L
+        if up_ind_x[i] != up_ind_y[i] :
+            l[up_ind_x[i]*L + up_ind_y[i]] = 1
+            l[up_ind_y[i]*L + up_ind_x[i]] = 1
+
+            A_ub = np.vstack((A_ub,l))
+            b_ub.append(1)
+
+            """for i in range(7):
+                print(l[i*7:i*7+7])
+            print("\n")"""
+
+    return A_ub, b_ub
 
 # Constraint to respect max number of travels
 def respect_number(landmarks, A_ub, b_ub):
     h = []
     for i in range(len(landmarks)) : h.append([1]*len(landmarks))
     T = block_diag(*h)
+    """for l in T :
+        for i in range(7):
+            print(l[i*7:i*7+7])
+        print("\n")"""
     return np.vstack((A_ub, T)), b_ub + [1]*len(landmarks)
 
 # Constraint to tie the problem together and have a connected path
@@ -99,6 +152,10 @@ def respect_order(landmarks: list, A_eq, b_eq):
             A_eq = np.vstack((A_eq,l))
             b_eq.append(0)
 
+            for i in range(7):
+                print(l[i*7:i*7+7])
+            print("\n")
+
     return A_eq, b_eq
 
 # Compute manhattan distance between 2 locations
@@ -111,12 +168,19 @@ def manhattan_distance(loc1: tuple, loc2: tuple):
 def init_eq_not_stay(landmarks): 
     L = len(landmarks)
     l = [0]*L*L
+
+
     for i in range(L) :
         for j in range(L) :
             if j == i :
                 l[j + i*L] = 1
+    l[L-1] = 1      # cannot skip from start to finish
     #A_eq = np.array([np.array(xi) for xi in A_eq])                  # Must convert A_eq into an np array
     l = np.array(np.array(l))
+
+    """for i in range(7):
+        print(l[i*7:i*7+7])"""
+
     return [l], [0]
 
 # Initialize A and c. Compute the distances from all landmarks to each other and store attractiveness
@@ -135,24 +199,37 @@ def init_ub_dist(landmarks: list, max_steps: int):
     c = c*len(landmarks)
     A_ub = []
     for line in A :
+        #print(line)
         A_ub += line
     return c, A_ub, [max_steps]
 
 # Go through the landmarks and force the optimizer to use landmarks where attractiveness is set to -1
 def respect_user_mustsee(landmarks: list, A_eq: list, b_eq: list) :
     L = len(landmarks)
+    H = 0       # sort of heuristic to get an idea of the number of steps needed
+    for i in landmarks : 
+        if i.name == "départ" : elem_prev = i              # list of all matches
     for i, elem in enumerate(landmarks) :
         if elem.attractiveness == -1 :
             l = [0]*L*L
             if elem.name != "arrivée" :
                 for j in range(L) :
                     l[j +i*L] = 1
+                    
             else :                          # This ensures we go to goal
                 for k in range(L-1) :
                         l[k*L+L-1] = 1  
+
+            H += manhattan_distance(elem.loc, elem_prev.loc)
+            elem_prev = elem
+
+            """for i in range(7):
+                print(l[i*7:i*7+7])
+            print("\n")"""
+
             A_eq = np.vstack((A_eq,l))
             b_eq.append(1)
-    return A_eq, b_eq
+    return A_eq, b_eq, H
 
 # Computes the path length given path matrix (dist_table) and a result
 def path_length(P: list, resx: list) :
@@ -161,27 +238,30 @@ def path_length(P: list, resx: list) :
 # Initialize all landmarks (+ start and goal). Order matters here
 landmarks = []
 landmarks.append(landmark("départ", -1, (0, 0)))
-landmarks.append(landmark("concorde", -1, (5,5)))
+landmarks.append(landmark("tour eiffel", 99, (0,2)))                           # PUT IN JSON
-landmarks.append(landmark("tour eiffel", 99, (1,1)))                           # PUT IN JSON
+landmarks.append(landmark("arc de triomphe", 99, (0,4)))
-landmarks.append(landmark("arc de triomphe", 99, (2,3)))
+landmarks.append(landmark("louvre", 99, (0,6)))
-landmarks.append(landmark("louvre", 70, (4,2)))
+landmarks.append(landmark("montmartre", 99, (0,10)))
-landmarks.append(landmark("montmartre", 20, (0,2)))
+landmarks.append(landmark("concorde", 99, (0,8)))
 landmarks.append(landmark("arrivée", -1, (0, 0)))
 
 
 
 # CONSTRAINT TO RESPECT MAX NUMBER OF STEPS
-max_steps = 25
+max_steps = 16
 
 # SET CONSTRAINTS FOR INEQUALITY
 c, A_ub, b_ub = init_ub_dist(landmarks, max_steps)              # Add the distances from each landmark to the other
 P = A_ub                                                        # store the paths for later. Needed to compute path length
 A_ub, b_ub = respect_number(landmarks, A_ub, b_ub)              # Respect max number of visits. 
 
+# TODO : Problems with circular symmetry
+A_ub, b_ub = break_sym(landmarks, A_ub, b_ub)                  # break the symmetry. Only use the upper diagonal values
+
 # SET CONSTRAINTS FOR EQUALITY
 A_eq, b_eq = init_eq_not_stay(landmarks)                       # Force solution not to stay in same place
-A_eq, b_eq = respect_user_mustsee(landmarks, A_eq, b_eq)       # Check if there are user_defined must_see. Also takes care of start/goal
+A_eq, b_eq, H = respect_user_mustsee(landmarks, A_eq, b_eq)       # Check if there are user_defined must_see. Also takes care of start/goal
-A_eq, b_eq = break_sym(landmarks, A_eq, b_eq)                  # break the symmetry. Only use the upper diagonal values
+
 A_eq, b_eq = respect_order(landmarks, A_eq, b_eq)              # Respect order of visit (only works when max_steps is limiting factor)
 
 # Bounds for variables (x can only be 0 or 1)
@@ -192,10 +272,17 @@ res = linprog(c, A_ub=A_ub, b_ub=b_ub, A_eq=A_eq, b_eq = b_eq, bounds=x_bounds,
 
 # Raise error if no solution is found
 if not res.success :
-    raise ValueError("No solution has been found, please adapt your max steps")
+    print(f"No solution has been found within given timeframe.\nMinimum steps to visit all must_see is : {H}")
+    # Override the max_steps using the heuristic
+    for i, val in enumerate(b_ub) :
+        if val == max_steps : b_ub[i] = H
+
+    # Solve problem again :
+    res = linprog(c, A_ub=A_ub, b_ub=b_ub, A_eq=A_eq, b_eq = b_eq, bounds=x_bounds, method='highs', integrality=3)
+
 
 # Print result
-print_res(res, P)
+print_res(res, landmarks, P)