From e600f40c1a991844730d8cf7352185ef02c360a0 Mon Sep 17 00:00:00 2001 From: Kilian Scheidecker Date: Mon, 20 May 2024 01:37:35 +0200 Subject: [PATCH] corrected symmetry cosntraint --- backend/src/optimizer.py | 153 ++++++++++++++++++++++++++++++--------- 1 file changed, 120 insertions(+), 33 deletions(-) diff --git a/backend/src/optimizer.py b/backend/src/optimizer.py index dac0df5..82d136a 100644 --- a/backend/src/optimizer.py +++ b/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) : - N = len(res) # length of res +def untangle(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_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): - # print(X[i*N:i*N+N]) - # print(order) + + #for i,x in enumerate(X) : X[i] = round(x,0) + + #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) + for idx in order : + 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 -def break_sym(landmarks, A_eq, b_eq): +# Constraint to not have d14 and d41 simultaneously +def break_sym(landmarks, A_ub, b_ub): 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 + upper_ind = np.triu_indices(L,0,L) - A_eq = np.vstack((A_eq,l)) - b_eq.append(0) + up_ind_x = upper_ind[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 + 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, (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("tour eiffel", 99, (0,2))) # PUT IN JSON +landmarks.append(landmark("arc de triomphe", 99, (0,4))) +landmarks.append(landmark("louvre", 99, (0,6))) +landmarks.append(landmark("montmartre", 99, (0,10))) +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 = break_sym(landmarks, A_eq, b_eq) # break the symmetry. Only use the upper diagonal values +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 = 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)