scgraph.bmssp
1from heapq import heappush, heappop 2from math import ceil, log 3 4inf = float("inf") 5 6 7class BmsspDataStructure: 8 """ 9 Data structure for inserting, updating and pulling the M smallest key-value pairs 10 together with a lower bound on the remaining values (or B if empty), as required by Alg. 3. 11 """ 12 13 def __init__(self, subset_size: int, upper_bound: int | float): 14 # subset_size: how many items to return per pull (must match Alg. 3 for level l -> Given as M) 15 self.subset_size = max(1, subset_size) 16 self.upper_bound = upper_bound 17 self.best = {} 18 self.heap = [] 19 20 def insert_key_value(self, key: int, value: int | float): 21 """ 22 Insert/refresh a key-value pair; keeps only the best value per key. 23 """ 24 if value < self.best.get(key, inf): 25 self.best[key] = value 26 heappush(self.heap, (value, key)) 27 28 def __pop_current__(self): 29 """ 30 Pop the current minimum key that matches self.best. 31 Returns None if heap is exhausted of current items. 32 """ 33 while self.heap: 34 value, key = heappop(self.heap) 35 if self.best.get(key, inf) == value: 36 self.best.pop(key, None) # Remove stale key 37 return key 38 return None 39 40 def is_empty(self) -> bool: 41 """ 42 Check for empty data structure. 43 """ 44 return len(self.best) == 0 45 46 def pull(self): 47 """ 48 Return (remaining_best, subset) where subset is up to self.subset_size keys with *globally* smallest values. 49 Remove the returned keys from the structure (matching Alg. 3 semantics). 50 remaining_best is the smallest value still present after removal, or self.upper_bound if empty. 51 """ 52 subset = set() 53 count = 0 54 55 # Take up to M distinct current keys 56 while count < self.subset_size: 57 key = self.__pop_current__() 58 if key is None: 59 break 60 subset.add(key) 61 count += 1 62 63 # Compute lower bound for remaining 64 remaining_best = ( 65 min(self.best.values()) if self.best else self.upper_bound 66 ) 67 return remaining_best, subset 68 69 70class BmsspSolver: 71 def __init__(self, graph: list[dict[int, int | float]], origin_id: int): 72 """ 73 Function: 74 75 - Initialize the BMSSP solver with a graph represented as an adjacency list. 76 77 Required Arguments: 78 79 - graph: 80 - Type: list[dict[int, int | float]] 81 - Description: The graph is represented as an adjacency list, where each node points to a dictionary of its neighbors and their edge weights. 82 - See: https://connor-makowski.github.io/scgraph/scgraph/graph.html#Graph.validate_graph 83 - origin_id: 84 - Type: int 85 - What: The ID of the starting node for the BMSSP algorithm. 86 87 Optional Arguments: 88 89 - None 90 """ 91 self.graph = graph 92 self.original_graph_length = len(graph) 93 self.graph_length = len(self.graph) 94 self.distance_matrix = [inf] * self.graph_length 95 self.predecessor = [-1] * self.graph_length 96 self.distance_matrix[origin_id] = 0 97 98 if self.graph_length <= 2: 99 raise ValueError("Your provided graph must have more than 2 nodes") 100 101 # Practical choices (k and t) based on n 102 self.pivot_relaxation_steps = max( 103 2, int(log(self.graph_length) ** (1 / 3.0)) 104 ) # k 105 # Modification: Change int to ceil 106 self.target_tree_depth = max( 107 2, ceil(log(self.graph_length) ** (2 / 3.0)) 108 ) # t 109 110 # Compute max_tree_depth based on k and t 111 self.max_tree_depth = int( 112 ceil( 113 log(max(2, self.graph_length)) / max(1, self.target_tree_depth) 114 ) 115 ) 116 117 # Run the solver algorithm 118 upper_bound, frontier = self.recursive_bmssp( 119 self.max_tree_depth, inf, {origin_id} 120 ) 121 122 def find_pivots( 123 self, upper_bound: int | float, frontier: set[int] 124 ) -> tuple[set[int], set[int]]: 125 """ 126 Function: 127 128 - Finds pivot sets pivots and temp_frontier according to Algorithm 1. 129 130 Required Arguments: 131 132 - upper_bound: 133 - Type: int | float 134 - What: The upper bound threshold (B) 135 - frontier: 136 - Type: set[int] 137 - What: Set of vertices (S) 138 139 Optional Arguments: 140 141 - None 142 143 Returns: 144 145 - pivots: 146 - Type: Set[int] 147 - What: Set of pivot vertices 148 - frontier: 149 - Type: Set[int] 150 - What: Return a new frontier set of vertices within the upper_bound 151 """ 152 temp_frontier = set(frontier) 153 prev_frontier = set(frontier) 154 155 # Multi-step limited relaxation from current frontier 156 for _ in range(self.pivot_relaxation_steps): 157 curr_frontier = set() 158 for prev_frontier_idx in prev_frontier: 159 prev_frontier_distance = self.distance_matrix[prev_frontier_idx] 160 for connection_idx, connection_distance in self.graph[ 161 prev_frontier_idx 162 ].items(): 163 new_distance = prev_frontier_distance + connection_distance 164 # Important: Allow equality on relaxations 165 if new_distance <= self.distance_matrix[connection_idx]: 166 # Addition: Add predecessor tracking 167 if new_distance < self.distance_matrix[connection_idx]: 168 self.predecessor[connection_idx] = prev_frontier_idx 169 self.distance_matrix[connection_idx] = new_distance 170 if new_distance < upper_bound: 171 curr_frontier.add(connection_idx) 172 temp_frontier.update(curr_frontier) 173 # If the search balloons, take the current frontier as pivots 174 if len(temp_frontier) > self.pivot_relaxation_steps * len(frontier): 175 pivots = set(frontier) 176 return pivots, temp_frontier 177 prev_frontier = curr_frontier 178 179 # Build tight-edge forest F on temp_frontier: edges (u -> v) with db[u] + w == db[v] 180 forest_adj = {i: set() for i in temp_frontier} 181 indegree = {i: 0 for i in temp_frontier} 182 for frontier_idx in temp_frontier: 183 frontier_distance = self.distance_matrix[frontier_idx] 184 for connection_idx, connection_distance in self.graph[ 185 frontier_idx 186 ].items(): 187 if ( 188 connection_idx in temp_frontier 189 and abs( 190 (frontier_distance + connection_distance) 191 - self.distance_matrix[connection_idx] 192 ) 193 < 1e-12 194 ): 195 # direction is frontier_idx -> connection_idx (parent to child) 196 forest_adj[frontier_idx].add(connection_idx) 197 indegree[connection_idx] += 1 198 199 # Non-sticky DFS that counts size of the reachable tree 200 def dfs_count(root): 201 seen = set() 202 stack = [root] 203 cnt = 0 204 while stack: 205 x = stack.pop() 206 if x in seen: 207 continue 208 seen.add(x) 209 cnt += 1 210 stack.extend(forest_adj[x]) 211 return cnt 212 213 pivots = set() 214 for frontier_idx in frontier: 215 if indegree.get(frontier_idx, 0) == 0: 216 size = dfs_count(frontier_idx) 217 if size >= self.pivot_relaxation_steps: 218 pivots.add(frontier_idx) 219 220 return pivots, temp_frontier 221 222 def base_case( 223 self, upper_bound: int | float, frontier: set[int] 224 ) -> tuple[int | float, set[int]]: 225 """ 226 Function: 227 228 - Implements Algorithm 2: Base Case of BMSSP 229 230 Required Arguments: 231 - upper_bound: 232 - Type: int | float 233 - frontier: 234 - Type: set 235 - What: Set with a single vertex x (complete) 236 237 Returns: 238 - new_upper_bound: 239 - Type: int | float 240 - What: The new upper bound for the search 241 - new_frontier: 242 - Type: set[int] 243 - What: Set of vertices v such that distance_matrix[v] < new_upper_bound 244 """ 245 assert len(frontier) == 1, "Frontier must be a singleton set" 246 first_frontier = next(iter(frontier)) 247 248 new_frontier = set() 249 heap = [] 250 heappush(heap, (self.distance_matrix[first_frontier], first_frontier)) 251 # Addition: Add visited check to prevent reprocessing and dropping into infinite loops 252 visited = set() 253 # grow until we exceed pivot_relaxation_steps (practical limit), as in Algorithm 2 254 while heap and len(new_frontier) < self.pivot_relaxation_steps + 1: 255 frontier_distance, frontier_idx = heappop(heap) 256 if frontier_idx in visited: 257 continue 258 visited.add(frontier_idx) 259 new_frontier.add(frontier_idx) 260 for connection_idx, connection_distance in self.graph[ 261 frontier_idx 262 ].items(): 263 new_distance = frontier_distance + connection_distance 264 if ( 265 new_distance <= self.distance_matrix[connection_idx] 266 and new_distance < upper_bound 267 ): 268 # Addition: Add predecessor tracking 269 if new_distance < self.distance_matrix[connection_idx]: 270 self.predecessor[connection_idx] = frontier_idx 271 self.distance_matrix[connection_idx] = new_distance 272 heappush(heap, (new_distance, connection_idx)) 273 274 # If we exceeded k, trim by new boundary B' = max db over visited and return new_frontier = {db < B'} 275 if len(new_frontier) > self.pivot_relaxation_steps: 276 new_upper_bound = max(self.distance_matrix[i] for i in new_frontier) 277 new_frontier = { 278 i 279 for i in new_frontier 280 if self.distance_matrix[i] < new_upper_bound 281 } 282 return new_upper_bound, new_frontier 283 else: 284 # Success for this base case: return current upper_bound unchanged and the completed set 285 return upper_bound, new_frontier 286 287 def recursive_bmssp( 288 self, recursion_depth: int, upper_bound: int | float, frontier: set[int] 289 ) -> tuple[int | float, set[int]]: 290 """ 291 Function: 292 293 - Implements Algorithm 3: Bounded Multi-Source Shortest Path (BMSSP) 294 295 Required Arguments: 296 297 - recursion_depth: 298 - Type: int 299 - What: The depth of the recursion 300 - upper_bound: 301 - Type: float 302 - What: The upper bound for the search 303 - frontier: 304 - Type: set[int] 305 - What: The set of vertices to explore 306 307 Returns: 308 309 - new_upper_bound: 310 - Type: int | float 311 - What: The new upper bound for the search 312 - new_frontier: 313 - Type: set[int] 314 - What: Set of vertices v such that distance_matrix[v] < new_upper_bound 315 """ 316 # Base case 317 if recursion_depth == 0: 318 return self.base_case(upper_bound, frontier) 319 320 # Step 4: Find pivots and temporary frontier 321 pivots, temp_frontier = self.find_pivots(upper_bound, frontier) 322 323 # Step 5–6: initialize data_struct with pivots 324 # subset_size = 2^((l-1) * t) 325 subset_size = 2 ** ((recursion_depth - 1) * self.target_tree_depth) 326 data_struct = BmsspDataStructure( 327 subset_size=subset_size, upper_bound=upper_bound 328 ) 329 for pivot in pivots: 330 data_struct.insert_key_value(pivot, self.distance_matrix[pivot]) 331 332 # Track new_frontier and B' according to Algorithm 3 333 new_frontier = set() 334 min_pivot_distance = min( 335 (self.distance_matrix[p] for p in pivots), default=upper_bound 336 ) 337 last_min_pivot_distance = min_pivot_distance 338 339 # Work budget that scales with level: k^(2*l*t) 340 # k = self.pivot_relaxation_steps 341 # t = self.target_tree_depth 342 work_budget = self.pivot_relaxation_steps ** ( 343 2 * recursion_depth * self.target_tree_depth 344 ) 345 346 # Main loop 347 while len(new_frontier) < work_budget and not data_struct.is_empty(): 348 # Step 10: Pull from data_struct: get data_struct_frontier_i and upper_bound_i 349 data_struct_frontier_bound_i, data_struct_frontier_i = ( 350 data_struct.pull() 351 ) 352 if not data_struct_frontier_i: 353 # data_struct is empty -> success at this level 354 break 355 356 # Step 11: Recurse on (l-1, data_struct_frontier_bound_i, data_struct_frontier_i) 357 last_min_pivot_distance_i, new_frontier_i = self.recursive_bmssp( 358 recursion_depth - 1, 359 data_struct_frontier_bound_i, 360 data_struct_frontier_i, 361 ) 362 363 # Track results 364 new_frontier.update(new_frontier_i) 365 last_min_pivot_distance = last_min_pivot_distance_i # If we later stop due to budget, we must return the last_min_pivot_distance 366 367 # Step 13: Initialize intermediate_frontier to batch-prepend 368 intermediate_frontier = set() 369 370 # Step 14–20: relax edges from new_frontier_i and enqueue into D or intermediate_frontier per their interval 371 for new_frontier_idx in new_frontier_i: 372 new_frontier_distance = self.distance_matrix[new_frontier_idx] 373 for connection_idx, connection_distance in self.graph[ 374 new_frontier_idx 375 ].items(): 376 # Addition: Avoid self-loops 377 if connection_idx == new_frontier_idx: 378 continue 379 new_distance = new_frontier_distance + connection_distance 380 if new_distance <= self.distance_matrix[connection_idx]: 381 # Addition: Add predecessor tracking 382 if new_distance < self.distance_matrix[connection_idx]: 383 self.predecessor[connection_idx] = new_frontier_idx 384 self.distance_matrix[connection_idx] = new_distance 385 # Insert based on which interval the new distance falls into 386 if ( 387 data_struct_frontier_bound_i 388 <= new_distance 389 < upper_bound 390 ): 391 data_struct.insert_key_value( 392 connection_idx, new_distance 393 ) 394 elif ( 395 last_min_pivot_distance_i 396 <= new_distance 397 < data_struct_frontier_bound_i 398 ): 399 intermediate_frontier.add( 400 (connection_idx, new_distance) 401 ) 402 403 # Step 21: Batch prepend intermediate_frontier plus filtered data_struct_frontier_i in last_min_pivot_distance_i, data_struct_frontier_bound_i) 404 data_struct_frontier_i_filtered = { 405 (x, self.distance_matrix[x]) 406 for x in data_struct_frontier_i 407 if last_min_pivot_distance_i 408 <= self.distance_matrix[x] 409 < data_struct_frontier_bound_i 410 } 411 for frontier_idx, frontier_distance in ( 412 intermediate_frontier | data_struct_frontier_i_filtered 413 ): 414 data_struct.insert_key_value(frontier_idx, frontier_distance) 415 416 # Step 22: Final return 417 return min(last_min_pivot_distance, upper_bound), new_frontier | { 418 v 419 for v in temp_frontier 420 if self.distance_matrix[v] < last_min_pivot_distance 421 } 422 423 424# Add a special test case directly to the bmssp algorithm to test it if this file is directly called 425if __name__ == "__main__": 426 graph = [{1: 1, 2: 1}, {2: 1, 3: 3}, {3: 1, 4: 2}, {4: 1}, {}] 427 solver = BmsspSolver(graph, 0) 428 if solver.distance_matrix != [0, 1, 1, 2, 3]: 429 print("BMSSP Test: Failed") 430 else: 431 print("BMSSP Test: Passed")
inf =
inf
class
BmsspDataStructure:
8class BmsspDataStructure: 9 """ 10 Data structure for inserting, updating and pulling the M smallest key-value pairs 11 together with a lower bound on the remaining values (or B if empty), as required by Alg. 3. 12 """ 13 14 def __init__(self, subset_size: int, upper_bound: int | float): 15 # subset_size: how many items to return per pull (must match Alg. 3 for level l -> Given as M) 16 self.subset_size = max(1, subset_size) 17 self.upper_bound = upper_bound 18 self.best = {} 19 self.heap = [] 20 21 def insert_key_value(self, key: int, value: int | float): 22 """ 23 Insert/refresh a key-value pair; keeps only the best value per key. 24 """ 25 if value < self.best.get(key, inf): 26 self.best[key] = value 27 heappush(self.heap, (value, key)) 28 29 def __pop_current__(self): 30 """ 31 Pop the current minimum key that matches self.best. 32 Returns None if heap is exhausted of current items. 33 """ 34 while self.heap: 35 value, key = heappop(self.heap) 36 if self.best.get(key, inf) == value: 37 self.best.pop(key, None) # Remove stale key 38 return key 39 return None 40 41 def is_empty(self) -> bool: 42 """ 43 Check for empty data structure. 44 """ 45 return len(self.best) == 0 46 47 def pull(self): 48 """ 49 Return (remaining_best, subset) where subset is up to self.subset_size keys with *globally* smallest values. 50 Remove the returned keys from the structure (matching Alg. 3 semantics). 51 remaining_best is the smallest value still present after removal, or self.upper_bound if empty. 52 """ 53 subset = set() 54 count = 0 55 56 # Take up to M distinct current keys 57 while count < self.subset_size: 58 key = self.__pop_current__() 59 if key is None: 60 break 61 subset.add(key) 62 count += 1 63 64 # Compute lower bound for remaining 65 remaining_best = ( 66 min(self.best.values()) if self.best else self.upper_bound 67 ) 68 return remaining_best, subset
Data structure for inserting, updating and pulling the M smallest key-value pairs together with a lower bound on the remaining values (or B if empty), as required by Alg. 3.
def
insert_key_value(self, key: int, value: int | float):
21 def insert_key_value(self, key: int, value: int | float): 22 """ 23 Insert/refresh a key-value pair; keeps only the best value per key. 24 """ 25 if value < self.best.get(key, inf): 26 self.best[key] = value 27 heappush(self.heap, (value, key))
Insert/refresh a key-value pair; keeps only the best value per key.
def
is_empty(self) -> bool:
41 def is_empty(self) -> bool: 42 """ 43 Check for empty data structure. 44 """ 45 return len(self.best) == 0
Check for empty data structure.
def
pull(self):
47 def pull(self): 48 """ 49 Return (remaining_best, subset) where subset is up to self.subset_size keys with *globally* smallest values. 50 Remove the returned keys from the structure (matching Alg. 3 semantics). 51 remaining_best is the smallest value still present after removal, or self.upper_bound if empty. 52 """ 53 subset = set() 54 count = 0 55 56 # Take up to M distinct current keys 57 while count < self.subset_size: 58 key = self.__pop_current__() 59 if key is None: 60 break 61 subset.add(key) 62 count += 1 63 64 # Compute lower bound for remaining 65 remaining_best = ( 66 min(self.best.values()) if self.best else self.upper_bound 67 ) 68 return remaining_best, subset
Return (remaining_best, subset) where subset is up to self.subset_size keys with globally smallest values. Remove the returned keys from the structure (matching Alg. 3 semantics). remaining_best is the smallest value still present after removal, or self.upper_bound if empty.
class
BmsspSolver:
71class BmsspSolver: 72 def __init__(self, graph: list[dict[int, int | float]], origin_id: int): 73 """ 74 Function: 75 76 - Initialize the BMSSP solver with a graph represented as an adjacency list. 77 78 Required Arguments: 79 80 - graph: 81 - Type: list[dict[int, int | float]] 82 - Description: The graph is represented as an adjacency list, where each node points to a dictionary of its neighbors and their edge weights. 83 - See: https://connor-makowski.github.io/scgraph/scgraph/graph.html#Graph.validate_graph 84 - origin_id: 85 - Type: int 86 - What: The ID of the starting node for the BMSSP algorithm. 87 88 Optional Arguments: 89 90 - None 91 """ 92 self.graph = graph 93 self.original_graph_length = len(graph) 94 self.graph_length = len(self.graph) 95 self.distance_matrix = [inf] * self.graph_length 96 self.predecessor = [-1] * self.graph_length 97 self.distance_matrix[origin_id] = 0 98 99 if self.graph_length <= 2: 100 raise ValueError("Your provided graph must have more than 2 nodes") 101 102 # Practical choices (k and t) based on n 103 self.pivot_relaxation_steps = max( 104 2, int(log(self.graph_length) ** (1 / 3.0)) 105 ) # k 106 # Modification: Change int to ceil 107 self.target_tree_depth = max( 108 2, ceil(log(self.graph_length) ** (2 / 3.0)) 109 ) # t 110 111 # Compute max_tree_depth based on k and t 112 self.max_tree_depth = int( 113 ceil( 114 log(max(2, self.graph_length)) / max(1, self.target_tree_depth) 115 ) 116 ) 117 118 # Run the solver algorithm 119 upper_bound, frontier = self.recursive_bmssp( 120 self.max_tree_depth, inf, {origin_id} 121 ) 122 123 def find_pivots( 124 self, upper_bound: int | float, frontier: set[int] 125 ) -> tuple[set[int], set[int]]: 126 """ 127 Function: 128 129 - Finds pivot sets pivots and temp_frontier according to Algorithm 1. 130 131 Required Arguments: 132 133 - upper_bound: 134 - Type: int | float 135 - What: The upper bound threshold (B) 136 - frontier: 137 - Type: set[int] 138 - What: Set of vertices (S) 139 140 Optional Arguments: 141 142 - None 143 144 Returns: 145 146 - pivots: 147 - Type: Set[int] 148 - What: Set of pivot vertices 149 - frontier: 150 - Type: Set[int] 151 - What: Return a new frontier set of vertices within the upper_bound 152 """ 153 temp_frontier = set(frontier) 154 prev_frontier = set(frontier) 155 156 # Multi-step limited relaxation from current frontier 157 for _ in range(self.pivot_relaxation_steps): 158 curr_frontier = set() 159 for prev_frontier_idx in prev_frontier: 160 prev_frontier_distance = self.distance_matrix[prev_frontier_idx] 161 for connection_idx, connection_distance in self.graph[ 162 prev_frontier_idx 163 ].items(): 164 new_distance = prev_frontier_distance + connection_distance 165 # Important: Allow equality on relaxations 166 if new_distance <= self.distance_matrix[connection_idx]: 167 # Addition: Add predecessor tracking 168 if new_distance < self.distance_matrix[connection_idx]: 169 self.predecessor[connection_idx] = prev_frontier_idx 170 self.distance_matrix[connection_idx] = new_distance 171 if new_distance < upper_bound: 172 curr_frontier.add(connection_idx) 173 temp_frontier.update(curr_frontier) 174 # If the search balloons, take the current frontier as pivots 175 if len(temp_frontier) > self.pivot_relaxation_steps * len(frontier): 176 pivots = set(frontier) 177 return pivots, temp_frontier 178 prev_frontier = curr_frontier 179 180 # Build tight-edge forest F on temp_frontier: edges (u -> v) with db[u] + w == db[v] 181 forest_adj = {i: set() for i in temp_frontier} 182 indegree = {i: 0 for i in temp_frontier} 183 for frontier_idx in temp_frontier: 184 frontier_distance = self.distance_matrix[frontier_idx] 185 for connection_idx, connection_distance in self.graph[ 186 frontier_idx 187 ].items(): 188 if ( 189 connection_idx in temp_frontier 190 and abs( 191 (frontier_distance + connection_distance) 192 - self.distance_matrix[connection_idx] 193 ) 194 < 1e-12 195 ): 196 # direction is frontier_idx -> connection_idx (parent to child) 197 forest_adj[frontier_idx].add(connection_idx) 198 indegree[connection_idx] += 1 199 200 # Non-sticky DFS that counts size of the reachable tree 201 def dfs_count(root): 202 seen = set() 203 stack = [root] 204 cnt = 0 205 while stack: 206 x = stack.pop() 207 if x in seen: 208 continue 209 seen.add(x) 210 cnt += 1 211 stack.extend(forest_adj[x]) 212 return cnt 213 214 pivots = set() 215 for frontier_idx in frontier: 216 if indegree.get(frontier_idx, 0) == 0: 217 size = dfs_count(frontier_idx) 218 if size >= self.pivot_relaxation_steps: 219 pivots.add(frontier_idx) 220 221 return pivots, temp_frontier 222 223 def base_case( 224 self, upper_bound: int | float, frontier: set[int] 225 ) -> tuple[int | float, set[int]]: 226 """ 227 Function: 228 229 - Implements Algorithm 2: Base Case of BMSSP 230 231 Required Arguments: 232 - upper_bound: 233 - Type: int | float 234 - frontier: 235 - Type: set 236 - What: Set with a single vertex x (complete) 237 238 Returns: 239 - new_upper_bound: 240 - Type: int | float 241 - What: The new upper bound for the search 242 - new_frontier: 243 - Type: set[int] 244 - What: Set of vertices v such that distance_matrix[v] < new_upper_bound 245 """ 246 assert len(frontier) == 1, "Frontier must be a singleton set" 247 first_frontier = next(iter(frontier)) 248 249 new_frontier = set() 250 heap = [] 251 heappush(heap, (self.distance_matrix[first_frontier], first_frontier)) 252 # Addition: Add visited check to prevent reprocessing and dropping into infinite loops 253 visited = set() 254 # grow until we exceed pivot_relaxation_steps (practical limit), as in Algorithm 2 255 while heap and len(new_frontier) < self.pivot_relaxation_steps + 1: 256 frontier_distance, frontier_idx = heappop(heap) 257 if frontier_idx in visited: 258 continue 259 visited.add(frontier_idx) 260 new_frontier.add(frontier_idx) 261 for connection_idx, connection_distance in self.graph[ 262 frontier_idx 263 ].items(): 264 new_distance = frontier_distance + connection_distance 265 if ( 266 new_distance <= self.distance_matrix[connection_idx] 267 and new_distance < upper_bound 268 ): 269 # Addition: Add predecessor tracking 270 if new_distance < self.distance_matrix[connection_idx]: 271 self.predecessor[connection_idx] = frontier_idx 272 self.distance_matrix[connection_idx] = new_distance 273 heappush(heap, (new_distance, connection_idx)) 274 275 # If we exceeded k, trim by new boundary B' = max db over visited and return new_frontier = {db < B'} 276 if len(new_frontier) > self.pivot_relaxation_steps: 277 new_upper_bound = max(self.distance_matrix[i] for i in new_frontier) 278 new_frontier = { 279 i 280 for i in new_frontier 281 if self.distance_matrix[i] < new_upper_bound 282 } 283 return new_upper_bound, new_frontier 284 else: 285 # Success for this base case: return current upper_bound unchanged and the completed set 286 return upper_bound, new_frontier 287 288 def recursive_bmssp( 289 self, recursion_depth: int, upper_bound: int | float, frontier: set[int] 290 ) -> tuple[int | float, set[int]]: 291 """ 292 Function: 293 294 - Implements Algorithm 3: Bounded Multi-Source Shortest Path (BMSSP) 295 296 Required Arguments: 297 298 - recursion_depth: 299 - Type: int 300 - What: The depth of the recursion 301 - upper_bound: 302 - Type: float 303 - What: The upper bound for the search 304 - frontier: 305 - Type: set[int] 306 - What: The set of vertices to explore 307 308 Returns: 309 310 - new_upper_bound: 311 - Type: int | float 312 - What: The new upper bound for the search 313 - new_frontier: 314 - Type: set[int] 315 - What: Set of vertices v such that distance_matrix[v] < new_upper_bound 316 """ 317 # Base case 318 if recursion_depth == 0: 319 return self.base_case(upper_bound, frontier) 320 321 # Step 4: Find pivots and temporary frontier 322 pivots, temp_frontier = self.find_pivots(upper_bound, frontier) 323 324 # Step 5–6: initialize data_struct with pivots 325 # subset_size = 2^((l-1) * t) 326 subset_size = 2 ** ((recursion_depth - 1) * self.target_tree_depth) 327 data_struct = BmsspDataStructure( 328 subset_size=subset_size, upper_bound=upper_bound 329 ) 330 for pivot in pivots: 331 data_struct.insert_key_value(pivot, self.distance_matrix[pivot]) 332 333 # Track new_frontier and B' according to Algorithm 3 334 new_frontier = set() 335 min_pivot_distance = min( 336 (self.distance_matrix[p] for p in pivots), default=upper_bound 337 ) 338 last_min_pivot_distance = min_pivot_distance 339 340 # Work budget that scales with level: k^(2*l*t) 341 # k = self.pivot_relaxation_steps 342 # t = self.target_tree_depth 343 work_budget = self.pivot_relaxation_steps ** ( 344 2 * recursion_depth * self.target_tree_depth 345 ) 346 347 # Main loop 348 while len(new_frontier) < work_budget and not data_struct.is_empty(): 349 # Step 10: Pull from data_struct: get data_struct_frontier_i and upper_bound_i 350 data_struct_frontier_bound_i, data_struct_frontier_i = ( 351 data_struct.pull() 352 ) 353 if not data_struct_frontier_i: 354 # data_struct is empty -> success at this level 355 break 356 357 # Step 11: Recurse on (l-1, data_struct_frontier_bound_i, data_struct_frontier_i) 358 last_min_pivot_distance_i, new_frontier_i = self.recursive_bmssp( 359 recursion_depth - 1, 360 data_struct_frontier_bound_i, 361 data_struct_frontier_i, 362 ) 363 364 # Track results 365 new_frontier.update(new_frontier_i) 366 last_min_pivot_distance = last_min_pivot_distance_i # If we later stop due to budget, we must return the last_min_pivot_distance 367 368 # Step 13: Initialize intermediate_frontier to batch-prepend 369 intermediate_frontier = set() 370 371 # Step 14–20: relax edges from new_frontier_i and enqueue into D or intermediate_frontier per their interval 372 for new_frontier_idx in new_frontier_i: 373 new_frontier_distance = self.distance_matrix[new_frontier_idx] 374 for connection_idx, connection_distance in self.graph[ 375 new_frontier_idx 376 ].items(): 377 # Addition: Avoid self-loops 378 if connection_idx == new_frontier_idx: 379 continue 380 new_distance = new_frontier_distance + connection_distance 381 if new_distance <= self.distance_matrix[connection_idx]: 382 # Addition: Add predecessor tracking 383 if new_distance < self.distance_matrix[connection_idx]: 384 self.predecessor[connection_idx] = new_frontier_idx 385 self.distance_matrix[connection_idx] = new_distance 386 # Insert based on which interval the new distance falls into 387 if ( 388 data_struct_frontier_bound_i 389 <= new_distance 390 < upper_bound 391 ): 392 data_struct.insert_key_value( 393 connection_idx, new_distance 394 ) 395 elif ( 396 last_min_pivot_distance_i 397 <= new_distance 398 < data_struct_frontier_bound_i 399 ): 400 intermediate_frontier.add( 401 (connection_idx, new_distance) 402 ) 403 404 # Step 21: Batch prepend intermediate_frontier plus filtered data_struct_frontier_i in last_min_pivot_distance_i, data_struct_frontier_bound_i) 405 data_struct_frontier_i_filtered = { 406 (x, self.distance_matrix[x]) 407 for x in data_struct_frontier_i 408 if last_min_pivot_distance_i 409 <= self.distance_matrix[x] 410 < data_struct_frontier_bound_i 411 } 412 for frontier_idx, frontier_distance in ( 413 intermediate_frontier | data_struct_frontier_i_filtered 414 ): 415 data_struct.insert_key_value(frontier_idx, frontier_distance) 416 417 # Step 22: Final return 418 return min(last_min_pivot_distance, upper_bound), new_frontier | { 419 v 420 for v in temp_frontier 421 if self.distance_matrix[v] < last_min_pivot_distance 422 }
BmsspSolver(graph: list[dict[int, int | float]], origin_id: int)
72 def __init__(self, graph: list[dict[int, int | float]], origin_id: int): 73 """ 74 Function: 75 76 - Initialize the BMSSP solver with a graph represented as an adjacency list. 77 78 Required Arguments: 79 80 - graph: 81 - Type: list[dict[int, int | float]] 82 - Description: The graph is represented as an adjacency list, where each node points to a dictionary of its neighbors and their edge weights. 83 - See: https://connor-makowski.github.io/scgraph/scgraph/graph.html#Graph.validate_graph 84 - origin_id: 85 - Type: int 86 - What: The ID of the starting node for the BMSSP algorithm. 87 88 Optional Arguments: 89 90 - None 91 """ 92 self.graph = graph 93 self.original_graph_length = len(graph) 94 self.graph_length = len(self.graph) 95 self.distance_matrix = [inf] * self.graph_length 96 self.predecessor = [-1] * self.graph_length 97 self.distance_matrix[origin_id] = 0 98 99 if self.graph_length <= 2: 100 raise ValueError("Your provided graph must have more than 2 nodes") 101 102 # Practical choices (k and t) based on n 103 self.pivot_relaxation_steps = max( 104 2, int(log(self.graph_length) ** (1 / 3.0)) 105 ) # k 106 # Modification: Change int to ceil 107 self.target_tree_depth = max( 108 2, ceil(log(self.graph_length) ** (2 / 3.0)) 109 ) # t 110 111 # Compute max_tree_depth based on k and t 112 self.max_tree_depth = int( 113 ceil( 114 log(max(2, self.graph_length)) / max(1, self.target_tree_depth) 115 ) 116 ) 117 118 # Run the solver algorithm 119 upper_bound, frontier = self.recursive_bmssp( 120 self.max_tree_depth, inf, {origin_id} 121 )
Function:
- Initialize the BMSSP solver with a graph represented as an adjacency list.
Required Arguments:
- graph:
- Type: list[dict[int, int | float]]
- Description: The graph is represented as an adjacency list, where each node points to a dictionary of its neighbors and their edge weights.
- See: https://connor-makowski.github.io/scgraph/scgraph/graph.html#Graph.validate_graph
- origin_id:
- Type: int
- What: The ID of the starting node for the BMSSP algorithm.
Optional Arguments:
- None
def
find_pivots( self, upper_bound: int | float, frontier: set[int]) -> tuple[set[int], set[int]]:
123 def find_pivots( 124 self, upper_bound: int | float, frontier: set[int] 125 ) -> tuple[set[int], set[int]]: 126 """ 127 Function: 128 129 - Finds pivot sets pivots and temp_frontier according to Algorithm 1. 130 131 Required Arguments: 132 133 - upper_bound: 134 - Type: int | float 135 - What: The upper bound threshold (B) 136 - frontier: 137 - Type: set[int] 138 - What: Set of vertices (S) 139 140 Optional Arguments: 141 142 - None 143 144 Returns: 145 146 - pivots: 147 - Type: Set[int] 148 - What: Set of pivot vertices 149 - frontier: 150 - Type: Set[int] 151 - What: Return a new frontier set of vertices within the upper_bound 152 """ 153 temp_frontier = set(frontier) 154 prev_frontier = set(frontier) 155 156 # Multi-step limited relaxation from current frontier 157 for _ in range(self.pivot_relaxation_steps): 158 curr_frontier = set() 159 for prev_frontier_idx in prev_frontier: 160 prev_frontier_distance = self.distance_matrix[prev_frontier_idx] 161 for connection_idx, connection_distance in self.graph[ 162 prev_frontier_idx 163 ].items(): 164 new_distance = prev_frontier_distance + connection_distance 165 # Important: Allow equality on relaxations 166 if new_distance <= self.distance_matrix[connection_idx]: 167 # Addition: Add predecessor tracking 168 if new_distance < self.distance_matrix[connection_idx]: 169 self.predecessor[connection_idx] = prev_frontier_idx 170 self.distance_matrix[connection_idx] = new_distance 171 if new_distance < upper_bound: 172 curr_frontier.add(connection_idx) 173 temp_frontier.update(curr_frontier) 174 # If the search balloons, take the current frontier as pivots 175 if len(temp_frontier) > self.pivot_relaxation_steps * len(frontier): 176 pivots = set(frontier) 177 return pivots, temp_frontier 178 prev_frontier = curr_frontier 179 180 # Build tight-edge forest F on temp_frontier: edges (u -> v) with db[u] + w == db[v] 181 forest_adj = {i: set() for i in temp_frontier} 182 indegree = {i: 0 for i in temp_frontier} 183 for frontier_idx in temp_frontier: 184 frontier_distance = self.distance_matrix[frontier_idx] 185 for connection_idx, connection_distance in self.graph[ 186 frontier_idx 187 ].items(): 188 if ( 189 connection_idx in temp_frontier 190 and abs( 191 (frontier_distance + connection_distance) 192 - self.distance_matrix[connection_idx] 193 ) 194 < 1e-12 195 ): 196 # direction is frontier_idx -> connection_idx (parent to child) 197 forest_adj[frontier_idx].add(connection_idx) 198 indegree[connection_idx] += 1 199 200 # Non-sticky DFS that counts size of the reachable tree 201 def dfs_count(root): 202 seen = set() 203 stack = [root] 204 cnt = 0 205 while stack: 206 x = stack.pop() 207 if x in seen: 208 continue 209 seen.add(x) 210 cnt += 1 211 stack.extend(forest_adj[x]) 212 return cnt 213 214 pivots = set() 215 for frontier_idx in frontier: 216 if indegree.get(frontier_idx, 0) == 0: 217 size = dfs_count(frontier_idx) 218 if size >= self.pivot_relaxation_steps: 219 pivots.add(frontier_idx) 220 221 return pivots, temp_frontier
Function:
- Finds pivot sets pivots and temp_frontier according to Algorithm 1.
Required Arguments:
- upper_bound:
- Type: int | float
- What: The upper bound threshold (B)
- frontier:
- Type: set[int]
- What: Set of vertices (S)
Optional Arguments:
- None
Returns:
- pivots:
- Type: Set[int]
- What: Set of pivot vertices
- frontier:
- Type: Set[int]
- What: Return a new frontier set of vertices within the upper_bound
def
base_case( self, upper_bound: int | float, frontier: set[int]) -> tuple[int | float, set[int]]:
223 def base_case( 224 self, upper_bound: int | float, frontier: set[int] 225 ) -> tuple[int | float, set[int]]: 226 """ 227 Function: 228 229 - Implements Algorithm 2: Base Case of BMSSP 230 231 Required Arguments: 232 - upper_bound: 233 - Type: int | float 234 - frontier: 235 - Type: set 236 - What: Set with a single vertex x (complete) 237 238 Returns: 239 - new_upper_bound: 240 - Type: int | float 241 - What: The new upper bound for the search 242 - new_frontier: 243 - Type: set[int] 244 - What: Set of vertices v such that distance_matrix[v] < new_upper_bound 245 """ 246 assert len(frontier) == 1, "Frontier must be a singleton set" 247 first_frontier = next(iter(frontier)) 248 249 new_frontier = set() 250 heap = [] 251 heappush(heap, (self.distance_matrix[first_frontier], first_frontier)) 252 # Addition: Add visited check to prevent reprocessing and dropping into infinite loops 253 visited = set() 254 # grow until we exceed pivot_relaxation_steps (practical limit), as in Algorithm 2 255 while heap and len(new_frontier) < self.pivot_relaxation_steps + 1: 256 frontier_distance, frontier_idx = heappop(heap) 257 if frontier_idx in visited: 258 continue 259 visited.add(frontier_idx) 260 new_frontier.add(frontier_idx) 261 for connection_idx, connection_distance in self.graph[ 262 frontier_idx 263 ].items(): 264 new_distance = frontier_distance + connection_distance 265 if ( 266 new_distance <= self.distance_matrix[connection_idx] 267 and new_distance < upper_bound 268 ): 269 # Addition: Add predecessor tracking 270 if new_distance < self.distance_matrix[connection_idx]: 271 self.predecessor[connection_idx] = frontier_idx 272 self.distance_matrix[connection_idx] = new_distance 273 heappush(heap, (new_distance, connection_idx)) 274 275 # If we exceeded k, trim by new boundary B' = max db over visited and return new_frontier = {db < B'} 276 if len(new_frontier) > self.pivot_relaxation_steps: 277 new_upper_bound = max(self.distance_matrix[i] for i in new_frontier) 278 new_frontier = { 279 i 280 for i in new_frontier 281 if self.distance_matrix[i] < new_upper_bound 282 } 283 return new_upper_bound, new_frontier 284 else: 285 # Success for this base case: return current upper_bound unchanged and the completed set 286 return upper_bound, new_frontier
Function:
- Implements Algorithm 2: Base Case of BMSSP
Required Arguments:
- upper_bound:
- Type: int | float
- frontier:
- Type: set
- What: Set with a single vertex x (complete)
Returns:
- new_upper_bound:
- Type: int | float
- What: The new upper bound for the search
- new_frontier:
- Type: set[int]
- What: Set of vertices v such that distance_matrix[v] < new_upper_bound
def
recursive_bmssp( self, recursion_depth: int, upper_bound: int | float, frontier: set[int]) -> tuple[int | float, set[int]]:
288 def recursive_bmssp( 289 self, recursion_depth: int, upper_bound: int | float, frontier: set[int] 290 ) -> tuple[int | float, set[int]]: 291 """ 292 Function: 293 294 - Implements Algorithm 3: Bounded Multi-Source Shortest Path (BMSSP) 295 296 Required Arguments: 297 298 - recursion_depth: 299 - Type: int 300 - What: The depth of the recursion 301 - upper_bound: 302 - Type: float 303 - What: The upper bound for the search 304 - frontier: 305 - Type: set[int] 306 - What: The set of vertices to explore 307 308 Returns: 309 310 - new_upper_bound: 311 - Type: int | float 312 - What: The new upper bound for the search 313 - new_frontier: 314 - Type: set[int] 315 - What: Set of vertices v such that distance_matrix[v] < new_upper_bound 316 """ 317 # Base case 318 if recursion_depth == 0: 319 return self.base_case(upper_bound, frontier) 320 321 # Step 4: Find pivots and temporary frontier 322 pivots, temp_frontier = self.find_pivots(upper_bound, frontier) 323 324 # Step 5–6: initialize data_struct with pivots 325 # subset_size = 2^((l-1) * t) 326 subset_size = 2 ** ((recursion_depth - 1) * self.target_tree_depth) 327 data_struct = BmsspDataStructure( 328 subset_size=subset_size, upper_bound=upper_bound 329 ) 330 for pivot in pivots: 331 data_struct.insert_key_value(pivot, self.distance_matrix[pivot]) 332 333 # Track new_frontier and B' according to Algorithm 3 334 new_frontier = set() 335 min_pivot_distance = min( 336 (self.distance_matrix[p] for p in pivots), default=upper_bound 337 ) 338 last_min_pivot_distance = min_pivot_distance 339 340 # Work budget that scales with level: k^(2*l*t) 341 # k = self.pivot_relaxation_steps 342 # t = self.target_tree_depth 343 work_budget = self.pivot_relaxation_steps ** ( 344 2 * recursion_depth * self.target_tree_depth 345 ) 346 347 # Main loop 348 while len(new_frontier) < work_budget and not data_struct.is_empty(): 349 # Step 10: Pull from data_struct: get data_struct_frontier_i and upper_bound_i 350 data_struct_frontier_bound_i, data_struct_frontier_i = ( 351 data_struct.pull() 352 ) 353 if not data_struct_frontier_i: 354 # data_struct is empty -> success at this level 355 break 356 357 # Step 11: Recurse on (l-1, data_struct_frontier_bound_i, data_struct_frontier_i) 358 last_min_pivot_distance_i, new_frontier_i = self.recursive_bmssp( 359 recursion_depth - 1, 360 data_struct_frontier_bound_i, 361 data_struct_frontier_i, 362 ) 363 364 # Track results 365 new_frontier.update(new_frontier_i) 366 last_min_pivot_distance = last_min_pivot_distance_i # If we later stop due to budget, we must return the last_min_pivot_distance 367 368 # Step 13: Initialize intermediate_frontier to batch-prepend 369 intermediate_frontier = set() 370 371 # Step 14–20: relax edges from new_frontier_i and enqueue into D or intermediate_frontier per their interval 372 for new_frontier_idx in new_frontier_i: 373 new_frontier_distance = self.distance_matrix[new_frontier_idx] 374 for connection_idx, connection_distance in self.graph[ 375 new_frontier_idx 376 ].items(): 377 # Addition: Avoid self-loops 378 if connection_idx == new_frontier_idx: 379 continue 380 new_distance = new_frontier_distance + connection_distance 381 if new_distance <= self.distance_matrix[connection_idx]: 382 # Addition: Add predecessor tracking 383 if new_distance < self.distance_matrix[connection_idx]: 384 self.predecessor[connection_idx] = new_frontier_idx 385 self.distance_matrix[connection_idx] = new_distance 386 # Insert based on which interval the new distance falls into 387 if ( 388 data_struct_frontier_bound_i 389 <= new_distance 390 < upper_bound 391 ): 392 data_struct.insert_key_value( 393 connection_idx, new_distance 394 ) 395 elif ( 396 last_min_pivot_distance_i 397 <= new_distance 398 < data_struct_frontier_bound_i 399 ): 400 intermediate_frontier.add( 401 (connection_idx, new_distance) 402 ) 403 404 # Step 21: Batch prepend intermediate_frontier plus filtered data_struct_frontier_i in last_min_pivot_distance_i, data_struct_frontier_bound_i) 405 data_struct_frontier_i_filtered = { 406 (x, self.distance_matrix[x]) 407 for x in data_struct_frontier_i 408 if last_min_pivot_distance_i 409 <= self.distance_matrix[x] 410 < data_struct_frontier_bound_i 411 } 412 for frontier_idx, frontier_distance in ( 413 intermediate_frontier | data_struct_frontier_i_filtered 414 ): 415 data_struct.insert_key_value(frontier_idx, frontier_distance) 416 417 # Step 22: Final return 418 return min(last_min_pivot_distance, upper_bound), new_frontier | { 419 v 420 for v in temp_frontier 421 if self.distance_matrix[v] < last_min_pivot_distance 422 }
Function:
- Implements Algorithm 3: Bounded Multi-Source Shortest Path (BMSSP)
Required Arguments:
- recursion_depth:
- Type: int
- What: The depth of the recursion
- upper_bound:
- Type: float
- What: The upper bound for the search
- frontier:
- Type: set[int]
- What: The set of vertices to explore
Returns:
- new_upper_bound:
- Type: int | float
- What: The new upper bound for the search
- new_frontier:
- Type: set[int]
- What: Set of vertices v such that distance_matrix[v] < new_upper_bound