diff --git a/backend/src/optimizer.py b/backend/src/optimizer.py new file mode 100644 index 0000000..dac0df5 --- /dev/null +++ b/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) + + + + + + + +