backend/src/optimizer.py

@@ -0,0 +1,206 @@
+from scipy.optimize import linprog
+import numpy as np
+from scipy.linalg import block_diag
+
+
+# Defines the landmark class (aka some place there is to visit)
+class landmark :
+    def __init__(self, name: str, attractiveness: int, loc: tuple):
+        self.name = name
+        self.attractiveness = attractiveness
+        self.loc = loc
+
+# 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) :
+    N = len(res)                # 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
+    visit_order = []
+    cnt = 0
+
+    if n_landmarks % 2 == 1 :                                     # if odd number of visited checkpoints
+        for i in range(L) :
+            for j in range(L) :
+                if res[i*L + j] == 1 :              # if index is 1
+                    cnt += 1                        # increment counter
+                    if cnt % 2 == 1 :               # if counter odd
+                        visit_order.append(i)
+                        visit_order.append(j)
+    else :                                   # if even number of ones
+        for i in range(L) :
+            for j in range(L) :
+                if res[i*L + j] == 1 :              # if index is one
+                    cnt += 1                        # increment counter
+                    if j % (L-1) == 0 :             # if last node
+                        visit_order.append(j)       # append only the last index
+                        return visit_order          # return
+                    if cnt % 2 == 1 : 
+                        visit_order.append(i)
+                        visit_order.append(j)
+    return visit_order
+
+# Just to print the result
+def print_res(res: list, P) :
+    X = abs(res.x)
+    order = untangle(X)
+
+    # 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):
+    #     print(X[i*N:i*N+N])
+    # 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) : 
+        if i in order : print('- ' + elem.name)
+
+    steps = path_length(P, abs(res.x))
+    print("\nSteps walked : " + str(steps))
+
+# Constraint to use only the upper triangular indices for travel
+def break_sym(landmarks, A_eq, b_eq):
+    L = len(landmarks)
+    l = [0]*L*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))
+    b_eq.append(0)
+
+    return A_eq, b_eq
+
+# 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)
+    return np.vstack((A_ub, T)), b_ub + [1]*len(landmarks)
+
+# Constraint to tie the problem together and have a connected path
+def respect_order(landmarks: list, A_eq, b_eq): 
+    N = len(landmarks)
+    for i in range(N-1) :     # Prevent stacked ones
+        if i == 0 :
+            continue
+        else : 
+            l = [0]*N
+            l[i] = -1
+            l = l*N
+            for j in range(N) :
+                l[i*N + j] = 1
+
+            A_eq = np.vstack((A_eq,l))
+            b_eq.append(0)
+
+    return A_eq, b_eq
+
+# Compute manhattan distance between 2 locations
+def manhattan_distance(loc1: tuple, loc2: tuple):
+    x1, y1 = loc1
+    x2, y2 = loc2
+    return abs(x1 - x2) + abs(y1 - y2)
+
+# Constraint to not stay in position
+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
+    #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))
+    return [l], [0]
+
+# Initialize A and c. Compute the distances from all landmarks to each other and store attractiveness
+# We want to maximize the sightseeing :  max(c) st. A*x < b   and   A_eq*x = b_eq
+def init_ub_dist(landmarks: list, max_steps: int):
+    # Objective function coefficients. a*x1 + b*x2 + c*x3 + ...
+    c = []
+    # Coefficients of inequality constraints (left-hand side)
+    A = []
+    for i, spot1 in enumerate(landmarks) :
+        dist_table = [0]*len(landmarks)
+        c.append(-spot1.attractiveness)
+        for j, spot2 in enumerate(landmarks) :
+            dist_table[j] = manhattan_distance(spot1.loc, spot2.loc)
+        A.append(dist_table)
+    c = c*len(landmarks)
+    A_ub = []
+    for line in A :
+        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)
+    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
+            A_eq = np.vstack((A_eq,l))
+            b_eq.append(1)
+    return A_eq, b_eq
+
+# Computes the path length given path matrix (dist_table) and a result
+def path_length(P: list, resx: list) :
+    return np.dot(P, resx)
+
+# 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, (1,1)))                           # PUT IN JSON
+landmarks.append(landmark("arc de triomphe", 99, (2,3)))
+landmarks.append(landmark("louvre", 70, (4,2)))
+landmarks.append(landmark("montmartre", 20, (0,2)))
+landmarks.append(landmark("arrivée", -1, (0, 0)))
+
+
+
+# CONSTRAINT TO RESPECT MAX NUMBER OF STEPS
+max_steps = 25
+
+# 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. 
+
+# 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 = 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)
+x_bounds = [(0, 1)] * len(c)
+
+# Solve linear programming problem
+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)
+
+# Raise error if no solution is found
+if not res.success :
+    raise ValueError("No solution has been found, please adapt your max steps")
+
+# Print result
+print_res(res, P)
+
+
+
+
+
+
+
+