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.

BmsspDataStructure(subset_size: int, upper_bound: int | float)
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 = []
subset_size
upper_bound
best
heap
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:

Optional Arguments:

  • None
graph
original_graph_length
graph_length
distance_matrix
predecessor
pivot_relaxation_steps
target_tree_depth
max_tree_depth
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