massive numpy optimization and more tests

backend/src/main.py

@@ -109,9 +109,12 @@ def new_trip(preferences: Preferences,
     start_time = time.time()
 
     # Second stage optimization
-    refined_tour = refiner.refine_optimization(landmarks, base_tour,
+    try :
+        refined_tour = refiner.refine_optimization(landmarks, base_tour,
                                                preferences.max_time_minute,
                                                preferences.detour_tolerance_minute)
+    except Exception as exc :
+        raise HTTPException(status_code=500, detail=f"An unexpected error occurred: {str(exc)}") from exc
 
     t_second_stage = time.time() - start_time
     logger.debug(f'Generating landmarks      : {round(t_generate_landmarks,3)} seconds')

backend/src/tests/test_main.py

@@ -35,7 +35,6 @@ def test_turckheim(client, request):    # pylint: disable=redefined-outer-name
             }
         )
     result = response.json()
-    print(result)
     landmarks = load_trip_landmarks(client, result['first_landmark_uuid'])
 
 
@@ -45,13 +44,16 @@ def test_turckheim(client, request):    # pylint: disable=redefined-outer-name
     # Add details to report
     log_trip_details(request, landmarks, result['total_time'], duration_minutes)
 
+    # for elem in landmarks :
+    #     print(elem)
+
     # checks :
     assert response.status_code == 200  # check for successful planning
     assert isinstance(landmarks, list)  # check that the return type is a list
     assert duration_minutes*0.8 < int(result['total_time']) < duration_minutes*1.2
     assert len(landmarks) > 2           # check that there is something to visit
     assert comp_time < 30, f"Computation time exceeded 30 seconds: {comp_time:.2f} seconds"
-    assert 2==3
+    # assert 2==3
 
 def test_bellecour(client, request) :   # pylint: disable=redefined-outer-name
     """
@@ -84,16 +86,15 @@ def test_bellecour(client, request) :   # pylint: disable=redefined-outer-name
     # Add details to report
     log_trip_details(request, landmarks, result['total_time'], duration_minutes)
 
-    for elem in landmarks :
-        print(elem)
-        print(elem.osm_id)
+    # for elem in landmarks :
+    #     print(elem)
 
     # checks :
     assert response.status_code == 200  # check for successful planning
-    assert duration_minutes*0.8 < int(result['total_time']) < duration_minutes*1.2
     assert comp_time < 30, f"Computation time exceeded 30 seconds: {comp_time:.2f} seconds"
+    assert duration_minutes*0.8 < int(result['total_time']) < duration_minutes*1.2
+    # assert 2 == 3
 
-'''
 
 def test_Paris(client, request) :   # pylint: disable=redefined-outer-name
     """
@@ -126,14 +127,13 @@ def test_Paris(client, request) :   # pylint: disable=redefined-outer-name
     # Add details to report
     log_trip_details(request, landmarks, result['total_time'], duration_minutes)
 
-    for elem in landmarks :
-        print(elem)
-        print(elem.osm_id)
+    # for elem in landmarks :
+    #     print(elem)
 
     # checks :
     assert response.status_code == 200  # check for successful planning
-    assert duration_minutes*0.8 < int(result['total_time']) < duration_minutes*1.2
     assert comp_time < 30, f"Computation time exceeded 30 seconds: {comp_time:.2f} seconds"
+    assert duration_minutes*0.8 < int(result['total_time']) < duration_minutes*1.2
 
 
 def test_New_York(client, request) :   # pylint: disable=redefined-outer-name
@@ -167,14 +167,14 @@ def test_New_York(client, request) :   # pylint: disable=redefined-outer-name
     # Add details to report
     log_trip_details(request, landmarks, result['total_time'], duration_minutes)
 
-    for elem in landmarks :
-        print(elem)
-        print(elem.osm_id)
+    # for elem in landmarks :
+    #     print(elem)
 
     # checks :
     assert response.status_code == 200  # check for successful planning
-    assert duration_minutes*0.8 < int(result['total_time']) < duration_minutes*1.2
     assert comp_time < 30, f"Computation time exceeded 30 seconds: {comp_time:.2f} seconds"
+    assert duration_minutes*0.8 < int(result['total_time']) < duration_minutes*1.2
+
 
 def test_shopping(client, request) :   # pylint: disable=redefined-outer-name
     """
@@ -207,11 +207,15 @@ def test_shopping(client, request) :   # pylint: disable=redefined-outer-name
     # Add details to report
     log_trip_details(request, landmarks, result['total_time'], duration_minutes)
 
+    # for elem in landmarks :
+    #     print(elem)
+
     # checks :
     assert response.status_code == 200  # check for successful planning
-    assert duration_minutes*0.8 < int(result['total_time']) < duration_minutes*1.2
     assert comp_time < 30, f"Computation time exceeded 30 seconds: {comp_time:.2f} seconds"
-'''
+    assert duration_minutes*0.8 < int(result['total_time']) < duration_minutes*1.2
+
+
 # def test_new_trip_single_prefs(client):
 #     response = client.post(
 #         "/trip/new",

backend/src/tests/test_utils.py

@@ -42,7 +42,7 @@ def fetch_landmark(client, landmark_uuid: str):
 
     try:
         json_data = response.json()
-        logger.info(f'API Response: {json_data}')
+        # logger.info(f'API Response: {json_data}')
     except ValueError as e:
         logger.error(f'Failed to parse response as JSON: {response.text}')
         raise HTTPException(status_code=500, detail="Invalid response format from API") from e

backend/src/utils/cluster_manager.py

@@ -12,6 +12,10 @@ from ..utils.get_time_separation import get_distance
 from ..constants import OSM_CACHE_DIR
 
 
+# silence the overpass logger
+logging.getLogger('OSMPythonTools').setLevel(level=logging.CRITICAL)
+
+
 class Cluster(BaseModel):
     """"
     A class representing an interesting area for shopping or sightseeing.
@@ -102,7 +106,6 @@ class ClusterManager:
                 points.append(coords)
 
             self.all_points = np.array(points)
-            self.valid = True
 
             # Apply DBSCAN to find clusters. Choose different settings for different cities.
             if self.cluster_type == 'shopping' and len(self.all_points) > 200 :
@@ -114,12 +117,17 @@ class ClusterManager:
 
             labels = dbscan.fit_predict(self.all_points)
 
-            # Separate clustered points and noise points
-            self.cluster_points = self.all_points[labels != -1]
-            self.cluster_labels = labels[labels != -1]
+            # Check that there are at least 2 different clusters
+            if len(set(labels)) > 2 :
+                self.logger.debug(f"Found {len(set(labels))} different clusters.")
+                # Separate clustered points and noise points
+                self.cluster_points = self.all_points[labels != -1]
+                self.cluster_labels = labels[labels != -1]
+                self.filter_clusters()      # ValueError here sometimes. I dont know why. # Filter the clusters to keep only the largest ones.
+                self.valid = True
 
-            # filter the clusters to keep only the largest ones
-            self.filter_clusters()        
+            else : 
+                self.valid = False      
 
 
     def generate_clusters(self) -> list[Landmark]:

backend/src/utils/get_time_separation.py

@@ -23,23 +23,23 @@ def get_time(p1: tuple[float, float], p2: tuple[float, float]) -> int:
     """
 
 
-    if p1 == p2:
-        return 0
-    else:
-        # Compute the distance in km along the surface of the Earth
-        # (assume spherical Earth)
-        # this is the haversine formula, stolen from stackoverflow
-        # in order to not use any external libraries
-        lat1, lon1 = radians(p1[0]), radians(p1[1])
-        lat2, lon2 = radians(p2[0]), radians(p2[1])
+    # if p1 == p2:
+    #     return 0
+    # else:
+    # Compute the distance in km along the surface of the Earth
+    # (assume spherical Earth)
+    # this is the haversine formula, stolen from stackoverflow
+    # in order to not use any external libraries
+    lat1, lon1 = radians(p1[0]), radians(p1[1])
+    lat2, lon2 = radians(p2[0]), radians(p2[1])
 
-        dlon = lon2 - lon1
-        dlat = lat2 - lat1
+    dlon = lon2 - lon1
+    dlat = lat2 - lat1
 
-        a = sin(dlat / 2)**2 + cos(lat1) * cos(lat2) * sin(dlon / 2)**2
-        c = 2 * atan2(sqrt(a), sqrt(1 - a))
+    a = sin(dlat / 2)**2 + cos(lat1) * cos(lat2) * sin(dlon / 2)**2
+    c = 2 * atan2(sqrt(a), sqrt(1 - a))
 
-        distance = EARTH_RADIUS_KM * c
+    distance = EARTH_RADIUS_KM * c
 
     # Consider the detour factor for average an average city
     walk_distance = distance * DETOUR_FACTOR
@@ -47,7 +47,7 @@ def get_time(p1: tuple[float, float], p2: tuple[float, float]) -> int:
     # Time to walk this distance (in minutes)
     walk_time = walk_distance / AVERAGE_WALKING_SPEED * 60
 
-    return round(walk_time)
+    return min(round(walk_time), 32765)
 
 
 def get_distance(p1: tuple[float, float], p2: tuple[float, float]) -> int:

backend/src/utils/landmarks_manager.py

@@ -53,6 +53,8 @@ class LandmarkManager:
         self.overpass = Overpass()
         CachingStrategy.use(JSON, cacheDir=OSM_CACHE_DIR)
 
+        self.logger.info('LandmakManager successfully initialized.')
+
 
     def generate_landmarks_list(self, center_coordinates: tuple[float, float], preferences: Preferences) -> tuple[list[Landmark], list[Landmark]]:
         """
@@ -71,7 +73,7 @@ class LandmarkManager:
         - A list of all existing landmarks.
         - A list of the most important landmarks based on the user's preferences.
         """
-
+        self.logger.debug('Starting to fetch landmarks...')
         max_walk_dist = (preferences.max_time_minute/2)/60*self.walking_speed*1000/self.detour_factor
         reachable_bbox_side = min(max_walk_dist, self.max_bbox_side)
 
@@ -83,25 +85,32 @@ class LandmarkManager:
 
         # list for sightseeing
         if preferences.sightseeing.score != 0:
+            self.logger.debug('Fetching sightseeing landmarks...')
             score_function = lambda score: score * 10 * preferences.sightseeing.score / 5
             current_landmarks = self.fetch_landmarks(bbox, self.amenity_selectors['sightseeing'], preferences.sightseeing.type, score_function)
             all_landmarks.update(current_landmarks)
+            self.logger.debug('Fetching sightseeing clusters...')
 
             # special pipeline for historic neighborhoods
             neighborhood_manager = ClusterManager(bbox, 'sightseeing')
             historic_clusters = neighborhood_manager.generate_clusters()
             all_landmarks.update(historic_clusters)
+            self.logger.debug('Sightseeing clusters done')
 
         # list for nature
         if preferences.nature.score != 0:
+            self.logger.debug('Fetching nature landmarks...')
             score_function = lambda score: score * 10 * self.nature_coeff * preferences.nature.score / 5
             current_landmarks = self.fetch_landmarks(bbox, self.amenity_selectors['nature'], preferences.nature.type, score_function)
             all_landmarks.update(current_landmarks)
 
+
         # list for shopping
         if preferences.shopping.score != 0:
+            self.logger.debug('Fetching shopping landmarks...')
             score_function = lambda score: score * 10 * preferences.shopping.score / 5
             current_landmarks = self.fetch_landmarks(bbox, self.amenity_selectors['shopping'], preferences.shopping.type, score_function)
+            self.logger.debug('Fetching shopping clusters...')
 
             # set time for all shopping activites :
             for landmark in current_landmarks : landmark.duration = 30
@@ -111,18 +120,19 @@ class LandmarkManager:
             shopping_manager = ClusterManager(bbox, 'shopping')
             shopping_clusters = shopping_manager.generate_clusters()
             all_landmarks.update(shopping_clusters)
+            self.logger.debug('Shopping clusters done')
 
 
 
         landmarks_constrained = take_most_important(all_landmarks, self.N_important)
-        self.logger.info(f'Generated {len(all_landmarks)} landmarks around {center_coordinates}, and constrained to {len(landmarks_constrained)} most important ones.')
+        self.logger.info(f'All landmarks generated : {len(all_landmarks)} landmarks around {center_coordinates}, and constrained to {len(landmarks_constrained)} most important ones.')
 
         return all_landmarks, landmarks_constrained
 
 
-
+    """
     def count_elements_close_to(self, coordinates: tuple[float, float]) -> int:
-        """
+        
         Count the number of OpenStreetMap elements (nodes, ways, relations) within a specified radius of the given location.
 
         This function constructs a bounding box around the specified coordinates based on the radius. It then queries
@@ -134,7 +144,7 @@ class LandmarkManager:
         Returns:
             int: The number of elements (nodes, ways, relations) within the specified radius. Returns 0 if no elements
                 are found or if an error occurs during the query.
-        """
+        
         
         lat = coordinates[0]
         lon = coordinates[1]
@@ -162,6 +172,7 @@ class LandmarkManager:
             return N_elem
         except:
             return 0
+    """
 
 
     # def create_bbox(self, coordinates: tuple[float, float], reachable_bbox_side: int) -> tuple[float, float, float, float]:
@@ -211,7 +222,7 @@ class LandmarkManager:
         # caution, when applying a list of selectors, overpass will search for elements that match ALL selectors simultaneously
         # we need to split the selectors into separate queries and merge the results
         for sel in dict_to_selector_list(amenity_selector):
-            self.logger.debug(f"Current selector: {sel}")
+            # self.logger.debug(f"Current selector: {sel}")
 
             element_types = ['way', 'relation']
 
@@ -230,7 +241,7 @@ class LandmarkManager:
                 includeCenter = True,
                 out = 'center'
                 )
-            self.logger.debug(f"Query: {query}")
+            # self.logger.debug(f"Query: {query}")
 
             try:
                 result = self.overpass.query(query)

backend/src/utils/optimizer.py

@@ -35,11 +35,7 @@ class Optimizer:
         """
         Initialize the objective function coefficients and inequality constraints.
         -> Adds 1 row of constraints
-
-                1 row
-            + L-1 rows
-
-        -> Pre-allocates A_ub for the rest of the computations
+        -> Pre-allocates A_ub for the rest of the computations with 2*L rows
 
         This function computes the distances between all landmarks and stores
         their attractiveness to maximize sightseeing. The goal is to maximize 
@@ -59,36 +55,42 @@ class Optimizer:
         c = np.zeros(L, dtype=np.int16)
 
         # Coefficients of inequality constraints (left-hand side)
-        A_first = np.zeros((L, L), dtype=np.int16)
+        A_ub = np.zeros((2*L, L*L), dtype=np.int16)
+        b_ub = np.zeros(2*L, dtype=np.int16)
+
+        # Fill in first row
+        b_ub[0] = round(max_time*self.overshoot)
 
         for i, spot1 in enumerate(landmarks) :
             c[i] = -spot1.attractiveness
             for j in range(i+1, L) :
                 if i !=j :
                     t = get_time(spot1.location, landmarks[j].location) + spot1.duration
-                    A_first[i,j] = t
-                    A_first[j,i] = t
+                    A_ub[0, i*L + j] = t
+                    A_ub[0, j*L + i] = t
+
+        # Expand 'c' to L*L for every decision variable
+        c = np.tile(c, L)     
 
         # Now sort and modify A_ub for each row
         if L > 22 :
             for i in range(L):
-                # Get indices of the 20 smallest values in row i
-                closest_indices = np.argpartition(A_first[i, :], 20)[:20]
+                # Get indices of the 4 smallest values in row i
+                row_values = A_ub[0, i*L:i*L+L]
+                closest_indices = np.argpartition(row_values, 22)[:22]
 
                 # Create a mask for non-closest landmarks
                 mask = np.ones(L, dtype=bool)
                 mask[closest_indices] = False
 
-                # Set non-closest landmarks to 32700
-                A_first[i, mask] = 32765
-        
-        # Replicate the objective function 'c' for each decision variable (L times)
-        c = np.tile(c, L)  # This correctly expands 'c' to L*L
+                # Set non-closest landmarks to 32765
+                row_values[mask] = 32765
+                A_ub[0, i*L:i*L+L] = row_values
 
-        return c, A_first.flatten(), [max_time*self.overshoot]
+        return c, A_ub, b_ub
 
 
-    def respect_number(self, L, max_landmarks: int):
+    def respect_number(self, A, b, L, max_landmarks: int):
         """
         Generate constraints to ensure each landmark is visited only once and cap the total number of visited landmarks.
         -> Adds L-1 rows of constraints
@@ -100,24 +102,25 @@ class Optimizer:
             tuple[np.ndarray, list[int]]: Inequality constraint coefficients and the right-hand side of the inequality constraints.
         """
         # First constraint: each landmark is visited exactly once
-        A = np.zeros((L-1, L*L), dtype=np.int8)
-        b = []
+        # A = np.zeros((L-1, L*L), dtype=np.int8)
+        # b = []
+        # Fill-in row 2 until row L-2
         for i in range(1, L-1):
-            A[i-1, L*i:L*(i+1)] = np.ones(L, dtype=np.int8)
-            b.append(1)
+            A[i, L*i:L*(i+1)] = np.ones(L, dtype=np.int16)
+            b[i] = 1
 
+        # Fill-in row L-1
         # Second constraint: cap the total number of visits
-        A[-1, :] = np.ones(L*L, dtype=np.int8)
-        b.append(max_landmarks+2)
-        return A, b
+        A[L-1, :] = np.ones(L*L, dtype=np.int16)
+        b[L-1] = max_landmarks+2
 
 
-    def break_sym(self, L):
+    def break_sym(self, A, b, L):
         """
         Generate constraints to prevent simultaneous travel between two landmarks
         in both directions. Constraint to not have d14 and d41 simultaneously. 
         Does not prevent cyclic paths with more elements
-        -> Adds a variable number of rows of constraints
+        -> Adds L rows of constraints (some of which might be zero)
 
         Args:
             L (int): Number of landmarks.
@@ -126,25 +129,27 @@ class Optimizer:
             tuple[np.ndarray, list[int]]:   Inequality constraint coefficients and 
                                             the right-hand side of the inequality constraints.
         """
-        b = []
+        # b = []
         upper_ind = np.triu_indices(L,0,L)
         up_ind_x = upper_ind[0]
         up_ind_y = upper_ind[1]
 
-        A = np.zeros((len(up_ind_x[1:]),L*L), dtype=np.int8)
-        for i, _ in enumerate(up_ind_x[1:]) :
+        # A = np.zeros((len(up_ind_x[1:]),L*L), dtype=np.int8)
+        # Fill-in rows L to 2*L-1
+        for i in range(L) :
             if up_ind_x[i] != up_ind_y[i] :
-                A[i, up_ind_x[i]*L + up_ind_y[i]] = 1
-                A[i, up_ind_y[i]*L + up_ind_x[i]] = 1
-                b.append(1)
+                A[L+i, up_ind_x[i]*L + up_ind_y[i]] = 1
+                A[L+i, up_ind_y[i]*L + up_ind_x[i]] = 1
+                b[L+i] = 1
 
-        return A[~np.all(A == 0, axis=1)], b
+        # return A[~np.all(A == 0, axis=1)], b
 
 
-    def init_eq_not_stay(self, L: int): 
+    def init_eq_not_stay(self, landmarks: list): 
         """
         Generate constraints to prevent staying in the same position (e.g., removing d11, d22, d33, etc.).
         -> Adds 1 row of constraints
+        -> Pre-allocates A_eq for the rest of the computations with (L+2 + dynamic incr) rows
 
         Args:
             L (int): Number of landmarks.
@@ -152,15 +157,78 @@ class Optimizer:
         Returns:
             tuple[list[np.ndarray], list[int]]: Equality constraint coefficients and the right-hand side of the equality constraints.
         """
+        L = len(landmarks)
+        incr = 0
+        for i, elem in enumerate(landmarks) :
+            if (elem.must_do or elem.must_avoid) and i not in [0, L-1]:
+                incr += 1
+
+        A_eq = np.zeros((L+2+incr, L*L), dtype=np.int8)
+        b_eq = np.zeros(L+2+incr, dtype=np.int8)
         l = np.zeros((L, L), dtype=np.int8)
         
         # Set diagonal elements to 1 (to prevent staying in the same position)
         np.fill_diagonal(l, 1)
 
-        return l.flatten(), [0]
+        # Fill-in first row
+        A_eq[0,:] = l.flatten()
+        b_eq[0] = 0
+
+        return A_eq, b_eq
 
 
-    def respect_user_must_do(self, landmarks: list[Landmark]) :
+    # Constraint to ensure start at start and finish at goal
+    def respect_start_finish(self, A, b, L: int):
+        """
+        Generate constraints to ensure that the optimization starts at the designated
+        start landmark and finishes at the goal landmark.
+        -> Adds 3 rows of constraints
+
+        Args:
+            L (int): Number of landmarks.
+
+        Returns:
+            tuple[np.ndarray, list[int]]: Inequality constraint coefficients and the right-hand side of the inequality constraints.
+        """
+        # Fill-in row 1.
+        A[1, :L] = np.ones(L, dtype=np.int8)        # sets departures only for start (horizontal ones)
+        for k in range(L-1) :
+            if k != 0 :
+                # Fill-in row 2
+                A[2, k*L+L-1] = 1                   # sets arrivals only for finish (vertical ones)
+            # Fill-in row 3
+            A[3, k*L] = 1
+
+        A[3, L*(L-1):] = np.ones(L, dtype=np.int8)  # prevents arrivals at start and departures from goal
+        b[1:4] = [1, 1, 0]
+
+        # return A, b
+
+
+    def respect_order(self, A, b, L: int): 
+        """
+        Generate constraints to tie the optimization problem together and prevent 
+        stacked ones, although this does not fully prevent circles.
+        -> Adds L-2 rows of constraints
+
+        Args:
+            L (int): Number of landmarks.
+
+        Returns:
+            tuple[np.ndarray, list[int]]: Inequality constraint coefficients and the right-hand side of the inequality constraints.
+        """
+        ones = np.ones(L, dtype=np.int8)
+
+        # Fill-in rows 4 to L+2
+        for i in range(1, L-1) :           # Prevent stacked ones
+            for j in range(L) :
+                A[i-1+4, i + j*L] = -1
+            A[i-1+4, i*L:(i+1)*L] = ones
+
+        b[4:L+2] = np.zeros(L-2, dtype=np.int8)
+
+
+    def respect_user_must(self, A, b, landmarks: list[Landmark]) :
         """
         Generate constraints to ensure that landmarks marked as 'must_do' are included in the optimization.
         -> Adds a variable number of rows of constraints BUT CAN BE PRE COMPUTED
@@ -173,88 +241,20 @@ class Optimizer:
             tuple[np.ndarray, list[int]]: Inequality constraint coefficients and the right-hand side of the inequality constraints.
         """
         L = len(landmarks)
-        A = np.zeros((L, L*L), dtype=np.int8)
-        b = []
-
-        for i, elem in enumerate(landmarks) :
-            if elem.must_do is True and elem.name not in ['finish', 'start']:
-                A[i, i*L:i*L+L] = np.ones(L, dtype=np.int8)
-                b.append(1)
-
-        return A[~np.all(A == 0, axis=1)], b
-    
-
-    def respect_user_must_avoid(self, landmarks: list[Landmark]) :
-        """
-        Generate constraints to ensure that landmarks marked as 'must_avoid' are skipped
-        in the optimization.
-        -> Adds a variable number of rows of constraints BUT CAN BE PRE COMPUTED
-
-        Args:
-            landmarks (list[Landmark]): List of landmarks, where some are marked as 'must_avoid'.
-
-        Returns:
-            tuple[np.ndarray, list[int]]: Inequality constraint coefficients and the right-hand side of the inequality constraints.
-        """
-        L = len(landmarks)
-        A = np.zeros((L, L*L), dtype=np.int8)
-        b = []
+        ones = np.ones(L, dtype=np.int8)
+        incr = 0
 
         for i, elem in enumerate(landmarks) :
             if elem.must_do is True and i not in [0, L-1]:
-                A[i, i*L:i*L+L] = np.ones(L, dtype=np.int8)
-                b.append(0)
-
-        return A[~np.all(A == 0, axis=1)], b
-
-
-    # Constraint to ensure start at start and finish at goal
-    def respect_start_finish(self, L: int):
-        """
-        Generate constraints to ensure that the optimization starts at the designated
-        start landmark and finishes at the goal landmark.
-        -> Adds 3 rows of constraints
-        Args:
-            L (int): Number of landmarks.
-
-        Returns:
-            tuple[np.ndarray, list[int]]: Inequality constraint coefficients and the right-hand side of the inequality constraints.
-        """
-
-        A = np.zeros((3, L*L), dtype=np.int8)
-
-        A[0, :L] = np.ones(L, dtype=np.int8)        # sets departures only for start (horizontal ones)
-        for k in range(L-1) :
-            A[2, k*L] = 1
-            if k != 0 :
-                A[1, k*L+L-1] = 1                   # sets arrivals only for finish (vertical ones)
-        A[2, L*(L-1):] = np.ones(L, dtype=np.int8)  # prevents arrivals at start and departures from goal
-        b = [1, 1, 0]
-
-        return A, b
-
-
-    def respect_order(self, L: int): 
-        """
-        Generate constraints to tie the optimization problem together and prevent 
-        stacked ones, although this does not fully prevent circles.
-        -> Adds L-2 rows of constraints
-
-        Args:
-            L (int): Number of landmarks.
-
-        Returns:
-            tuple[np.ndarray, list[int]]: Inequality constraint coefficients and the right-hand side of the inequality constraints.
-        """
-
-        A = np.zeros((L-2, L*L), dtype=np.int8)
-        b = [0]*(L-2)
-        for i in range(1, L-1) :           # Prevent stacked ones
-            for j in range(L) :
-                A[i-1, i + j*L] = -1
-            A[i-1, i*L:(i+1)*L] = np.ones(L, dtype=np.int8)
-
-        return A, b
+                # First part of the dynamic infill
+                A[L+2+incr, i*L:i*L+L] = ones
+                b[L+2+incr] = 1
+                incr += 1
+            if elem.must_avoid is True and i not in [0, L-1]:
+                # Second part of the dynamic infill
+                A[L+2+incr, i*L:i*L+L] = ones
+                b[L+2+incr] = 0
+                incr += 1
 
 
     # Prevent the use of a particular solution
@@ -282,13 +282,14 @@ class Optimizer:
         vertices_visited = ind_a
         vertices_visited.remove(0)
 
-        ones = [1]*L
-        h = [0]*N
+        ones = np.ones(L, dtype=np.int8)
+        h = np.zeros(L*L, dtype=np.int8)
+
         for i in range(L) :
             if i in vertices_visited :
                 h[i*L:i*L+L] = ones
 
-        return h, [len(vertices_visited)-1]
+        return h, np.array([len(vertices_visited)-1])
 
 
     # Prevents the creation of the same circle (both directions)
@@ -303,22 +304,21 @@ class Optimizer:
         Returns:
             tuple[np.ndarray, list[int]]: A tuple containing a new row for constraint matrix and new value for upper bound vector.
         """
+        l = np.zeros((2, L*L), dtype=np.int8)
 
-        l1 = [0]*L*L
-        l2 = [0]*L*L
         for i, node in enumerate(circle_vertices[:-1]) :
             next = circle_vertices[i+1]
 
-            l1[node*L + next] = 1
-            l2[next*L + node] = 1
+            l[0, node*L + next] = 1
+            l[1, next*L + node] = 1
 
         s = circle_vertices[0]
         g = circle_vertices[-1]
 
-        l1[g*L + s] = 1
-        l2[s*L + g] = 1
+        l[0, g*L + s] = 1
+        l[1, s*L + g] = 1
 
-        return np.vstack((l1, l2)), [0, 0]
+        return l, np.zeros(2, dtype=np.int8)
 
 
     def is_connected(self, resx) :
@@ -331,10 +331,90 @@ class Optimizer:
         Returns:
             tuple[list[int], Optional[list[list[int]]]]: A tuple containing the visit order and a list of any detected circles.
         """
+        resx = np.round(resx).astype(np.int8)  # round all elements and cast them to int
+        
+        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.
+
+        nonzeroind = np.nonzero(resx)[0] # the return is a little funny so I use the [0]
+        nonzero_tup = np.unravel_index(nonzeroind, (L,L))
+
+        ind_a = nonzero_tup[0]  # removed .tolist()
+        ind_b = nonzero_tup[1]
+
+        # Extract all journeys
+        all_journeys_nodes = []
+        visited_nodes = set()
+
+        for node in ind_a:
+            if node not in visited_nodes:
+                journey_nodes = self.get_journey(node, ind_a, ind_b)
+                all_journeys_nodes.append(journey_nodes)
+                visited_nodes.update(journey_nodes)
+
+        for l in all_journeys_nodes :
+            if 0 in l :
+                all_journeys_nodes.remove(l)
+                break
+
+        if not all_journeys_nodes :
+            return None
+
+        return all_journeys_nodes
+
+
+    def get_journey(self, start, ind_a, ind_b):
+        """
+        Trace the journey starting from a given node and follow the connections between landmarks.
+        This method constructs a graph from two lists of landmark connections, `ind_a` and `ind_b`,
+        where each element in `ind_a` is connected to the corresponding element in `ind_b`.
+        It then performs a depth-first search (DFS) starting from the `start` node to determine
+        the path (journey) by following the connections.
+
+        Args:
+            start (int): The starting node of the journey.
+            ind_a (list[int]): List of "from" nodes, representing the starting points of each connection.
+            ind_b (list[int]): List of "to" nodes, representing the endpoints of each connection.
+        
+        Returns:
+            list[int]: A list of nodes representing the order of the journey, starting from the `start` node.
+
+        Example:
+            If `ind_a = [0, 1, 2]` and `ind_b = [1, 2, 3]`, starting from node 0, the journey would be `[0, 1, 2, 3]`.
+        """
+        graph = defaultdict(list)
+        for a, b in zip(ind_a, ind_b):
+            graph[a].append(b)
+
+        journey_nodes = []
+        visited = set()
+        stack = deque([start])
+
+        while stack:
+            node = stack.pop()
+            if node not in visited:
+                visited.add(node)
+                journey_nodes.append(node)
+                for neighbor in graph[node]:
+                    if neighbor not in visited:
+                        stack.append(neighbor)
+
+        return journey_nodes
+        
+
+    def get_order(self, resx):
+        """
+        Determine the order of visits given the result of the optimization.
+
+        Args:
+            resx (list): List of edge weights.
+
+        Returns:
+            list[int]: A list containing the visit order.
+        """
 
         # first round the results to have only 0-1 values
-        for i, elem in enumerate(resx):
-            resx[i] = round(elem)
+        resx = np.round(resx).astype(np.uint8)  # round all elements and cast them to int
         
         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.
@@ -345,48 +425,32 @@ class Optimizer:
         ind_a = nonzero_tup[0].tolist()
         ind_b = nonzero_tup[1].tolist()
 
-        # Step 1: Create a graph representation
-        graph = defaultdict(list)
-        for a, b in zip(ind_a, ind_b):
-            graph[a].append(b)
+        order = [0]
+        current = 0
+        used_indices = set()  # Track visited index pairs
+        
+        while True:
+            # Find index of the current node in ind_a
+            try:
+                i = ind_a.index(current)
+            except ValueError:
+                break  # No more links, stop the search
+            
+            if i in used_indices:
+                break  # Prevent infinite loops
 
-        # Step 2: Function to perform BFS/DFS to extract journeys
-        def get_journey(start):
-            journey_nodes = []
-            visited = set()
-            stack = deque([start])
+            used_indices.add(i)  # Mark this index as visited
+            next_node = ind_b[i]  # Get the corresponding node in ind_b
+            order.append(next_node)  # Add it to the path
 
-            while stack:
-                node = stack.pop()
-                if node not in visited:
-                    visited.add(node)
-                    journey_nodes.append(node)
-                    for neighbor in graph[node]:
-                        if neighbor not in visited:
-                            stack.append(neighbor)
+            # Switch roles, now look for next_node in ind_a
+            try:
+                current = next_node
+            except ValueError:
+                break  # No further connections, end the path
 
-            return journey_nodes
+        return order
 
-        # Step 3: Extract all journeys
-        all_journeys_nodes = []
-        visited_nodes = set()
-
-        for node in ind_a:
-            if node not in visited_nodes:
-                journey_nodes = get_journey(node)
-                all_journeys_nodes.append(journey_nodes)
-                visited_nodes.update(journey_nodes)
-
-        for l in all_journeys_nodes :
-            if 0 in l :
-                order = l
-                all_journeys_nodes.remove(l)
-                break
-
-        if len(all_journeys_nodes) == 0 :
-            return order, None
-
-        return order, all_journeys_nodes
 
 
     def link_list(self, order: list[int], landmarks: list[Landmark])->list[Landmark] :
@@ -449,33 +513,34 @@ class Optimizer:
         L = len(landmarks)
 
         # SET CONSTRAINTS FOR INEQUALITY
-        c, A_ub, b_ub = self.init_ub_time(landmarks, max_time)          # Add the distances from each landmark to the other
-        
-        A, b = self.respect_number(L, max_landmarks)                    # Respect max number of visits (no more possible stops than landmarks). 
-        A_ub = np.vstack((A_ub, A))
-        b_ub += b
-
-        A, b = self.break_sym(L)                                         # break the 'zig-zag' symmetry
-        A_ub = np.vstack((A_ub, A))
-        b_ub += b
-
+        c, A_ub, b_ub = self.init_ub_time(landmarks, max_time)  # Adds the distances from each landmark to the other.
+        self.respect_number(A_ub, b_ub, L, max_landmarks)       # Respects max number of visits (no more possible stops than landmarks).
+        self.break_sym(A_ub, b_ub, L)                           # Breaks the 'zig-zag' symmetry. Avoids d12 and d21 but not larger cirlces.
 
         # SET CONSTRAINTS FOR EQUALITY
-        A_eq, b_eq = self.init_eq_not_stay(L)                            # Force solution not to stay in same place
-        A, b = self.respect_user_must_do(landmarks)                      # Check if there are user_defined must_see. Also takes care of start/goal
-        if len(b) > 0 :
-            A_eq = np.vstack((A_eq, A), dtype=np.int8)
-            b_eq += b
-        A, b = self.respect_user_must_avoid(landmarks)                      # Check if there are user_defined must_see. Also takes care of start/goal
-        if len(b) > 0 :
-            A_eq = np.vstack((A_eq, A), dtype=np.int8)
-            b_eq += b
-        A, b = self.respect_start_finish(L)                  # Force start and finish positions
-        A_eq = np.vstack((A_eq, A), dtype=np.int8)
-        b_eq += b
-        A, b = self.respect_order(L)                         # Respect order of visit (only works when max_time is limiting factor)
-        A_eq = np.vstack((A_eq, A), dtype=np.int8)
-        b_eq += b
+        A_eq, b_eq = self.init_eq_not_stay(landmarks)           # Force solution not to stay in same place
+        self.respect_start_finish(A_eq, b_eq, L)                # Force start and finish positions
+        self.respect_order(A_eq, b_eq, L)                       # Respect order of visit (only works when max_time is limiting factor)
+        self.respect_user_must(A_eq, b_eq, landmarks)           # Force to do/avoid landmarks set by user.
+
+        self.logger.debug(f"Optimizing with {A_ub.shape[0]} + {A_eq.shape[0]} = {A_ub.shape[0] + A_eq.shape[0]} constraints.")
+
+
+
+        # A, b = self.respect_user_must_do(landmarks)                      # Check if there are user_defined must_see. Also takes care of start/goal
+        # if len(b) > 0 :
+        #     A_eq = np.vstack((A_eq, A), dtype=np.int8)
+        #     b_eq += b
+        # A, b = self.respect_user_must_avoid(landmarks)                      # Check if there are user_defined must_see. Also takes care of start/goal
+        # if len(b) > 0 :
+        #     A_eq = np.vstack((A_eq, A), dtype=np.int8)
+        #     b_eq += b
+        # A, b = self.respect_start_finish(L)                  # Force start and finish positions
+        # A_eq = np.vstack((A_eq, A), dtype=np.int8)
+        # b_eq += b
+        # A, b = self.respect_order(L)                         # Respect order of visit (only works when max_time is limiting factor)
+        # A_eq = np.vstack((A_eq, A), dtype=np.int8)
+        # b_eq += b
         # until here opti
         
         # SET BOUNDS FOR DECISION VARIABLE (x can only be 0 or 1)
@@ -484,39 +549,47 @@ class Optimizer:
         # 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)
 
+        self.logger.debug("First results are out. Looking out for circles and correcting.")
+
         # Raise error if no solution is found. FIXME: for now this throws the internal server error
         if not res.success :
+            self.logger.error("The problem is overconstrained, no solution on first try.")
             raise ArithmeticError("No solution could be found, the problem is overconstrained. Try with a longer trip (>30 minutes).")
 
         # If there is a solution, we're good to go, just check for connectiveness
-        order, circles = self.is_connected(res.x)
+        circles = self.is_connected(res.x)
         #nodes, edges = is_connected(res.x)
         i = 0
         timeout = 80
         while circles is not None and i < timeout:
+            i += 1
+            # print(f"Iteration {i} of fixing circles")
             A, b = self.prevent_config(res.x)
             A_ub = np.vstack((A_ub, A))
-            b_ub += b
-            #A_ub, b_ub = prevent_circle(order, len(landmarks), A_ub, b_ub)
+            b_ub = np.concatenate((b_ub, b))
+
             for circle in circles :
                 A, b = self.prevent_circle(circle, L)
                 A_eq = np.vstack((A_eq, A))
-                b_eq += b
+                b_eq = np.concatenate((b_eq, b))
+
             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)
+
             if not res.success :
+                self.logger.error(f'Unexpected error after {timeout} iterations of fixing circles.')
                 raise ArithmeticError("Solving failed because of overconstrained problem")
-            order, circles = self.is_connected(res.x)
+            circles = self.is_connected(res.x)
             #nodes, edges = is_connected(res.x)
             if circles is None :
                 break
-            # print(i)
-            i += 1
-        
+
         if i == timeout :
+            self.logger.error(f'Timeout: No solution found after {timeout} iterations.')
             raise TimeoutError(f"Optimization took too long. No solution found after {timeout} iterations.")
 
-        #sort the landmarks in the order of the solution
+        # Sort the landmarks in the order of the solution
+        order = self.get_order(res.x)
         tour =  [landmarks[i] for i in order] 
-        
+
         self.logger.debug(f"Re-optimized {i} times, score: {int(-res.fun)}")
         return tour