|
1 | | -from graphblas import Matrix, Vector, binary |
2 | | -from graphblas.select import offdiag |
3 | | -from graphblas.semiring import any_plus |
| 1 | +from graphblas import Matrix, Vector, binary, indexunary, replace, select |
| 2 | +from graphblas.semiring import any_plus, any_second |
4 | 3 |
|
5 | | -__all__ = ["floyd_warshall"] |
| 4 | +__all__ = ["floyd_warshall", "floyd_warshall_predecessor_and_distance"] |
6 | 5 |
|
7 | 6 |
|
8 | 7 | def floyd_warshall(G, is_weighted=False): |
| 8 | + return floyd_warshall_predecessor_and_distance(G, is_weighted, compute_predecessors=False)[1] |
| 9 | + |
| 10 | + |
| 11 | +def floyd_warshall_predecessor_and_distance(G, is_weighted=False, *, compute_predecessors=True): |
9 | 12 | # By using `offdiag` instead of `G._A`, we ensure that D will not become dense. |
10 | 13 | # Dense D may be better at times, but not including the diagonal will result in less work. |
11 | 14 | # Typically, Floyd-Warshall algorithms sets the diagonal of D to 0 at the beginning. |
12 | 15 | # This is unnecessary with sparse matrices, and we set the diagonal to 0 at the end. |
13 | 16 | # We also don't iterate over index `i` if either row i or column i are empty. |
14 | | - if G.is_directed(): |
| 17 | + if is_directed:=G.is_directed(): |
15 | 18 | A, row_degrees, column_degrees = G.get_properties("offdiag row_degrees- column_degrees-") |
16 | 19 | nonempty_nodes = binary.pair(row_degrees & column_degrees).new(name="nonempty_nodes") |
17 | 20 | else: |
18 | | - A, nonempty_nodes = G.get_properties("offdiag degrees-") |
| 21 | + A, nonempty_nodes = G.get_properties("U- degrees-") |
19 | 22 |
|
20 | 23 | if A.dtype == bool or not is_weighted: |
21 | 24 | dtype = int |
22 | 25 | else: |
23 | 26 | dtype = A.dtype |
24 | 27 | n = A.nrows |
25 | | - D = Matrix(dtype, nrows=n, ncols=n, name="floyd_warshall") |
| 28 | + D = Matrix(dtype, nrows=n, ncols=n, name="floyd_warshall_dist") |
26 | 29 | if is_weighted: |
27 | 30 | D << A |
28 | 31 | else: |
29 | 32 | D(A.S) << 1 # Like `D << unary.one[int](A)` |
30 | 33 | del A |
31 | | - |
32 | 34 | Row = Matrix(dtype, nrows=1, ncols=n, name="Row") |
33 | | - Col = Matrix(dtype, nrows=n, ncols=1, name="Col") |
| 35 | + if is_directed: |
| 36 | + Col = Matrix(dtype, nrows=n, ncols=1, name="Col") |
| 37 | + else: |
| 38 | + Col = None |
34 | 39 | Outer = Matrix(dtype, nrows=n, ncols=n, name="Outer") |
| 40 | + if compute_predecessors: |
| 41 | + Mask = Matrix(bool, nrows=n, ncols=n, name="Mask") |
| 42 | + P = indexunary.rowindex(D).new(name="floyd_warshall_pred") |
| 43 | + if P.dtype == dtype: |
| 44 | + P_row = Row |
| 45 | + else: |
| 46 | + P_row = Matrix(P.dtype, nrows=1, ncols=n, name="P_row") |
| 47 | + else: |
| 48 | + Mask = P = P_row = None |
| 49 | + |
35 | 50 | for i in nonempty_nodes: |
36 | | - Col << D[:, [i]] |
37 | 51 | Row << D[[i], :] |
| 52 | + if is_directed: |
| 53 | + Col << D[:, [i]] |
| 54 | + else: |
| 55 | + Row(binary.any) << D.T[[i], :] |
| 56 | + Col = Row.T |
38 | 57 | Outer << any_plus(Col @ Row) # Like `col.outer(row, binary.plus)` |
39 | | - D(binary.min) << offdiag(Outer) |
| 58 | + |
| 59 | + if not compute_predecessors: |
| 60 | + # It is faster (approx 10%-30%) to use a mask as is done below when computing |
| 61 | + # predecessors, but we choose to use less memory here by not using a mask. |
| 62 | + if is_directed: |
| 63 | + D(binary.min) << select.offdiag(Outer) |
| 64 | + else: |
| 65 | + D(binary.min) << select.triu(Outer, 1) |
| 66 | + else: |
| 67 | + # Update Outer to only include off-diagonal values that will update D and P. |
| 68 | + if is_directed: |
| 69 | + Mask << indexunary.offdiag(Outer) |
| 70 | + else: |
| 71 | + Mask << indexunary.triu(Outer, 1) |
| 72 | + Mask(binary.second) << binary.lt(Outer & D) |
| 73 | + Outer(Mask.V, replace) << Outer |
| 74 | + |
| 75 | + # Update distances; like `D(binary.min) << offdiag(any_plus(Col @ Row))` |
| 76 | + D(Outer.S) << Outer |
| 77 | + |
| 78 | + # Broadcast predecessors in P_row to updated values |
| 79 | + P_row << P[[i], :] |
| 80 | + if not is_directed: |
| 81 | + P_row(binary.any) << P.T[[i], :] |
| 82 | + Col = P_row.T |
| 83 | + P(Outer.S) << any_second(Col @ P_row) |
| 84 | + del Outer, Mask, Col, Row, P_row |
| 85 | + |
| 86 | + if not is_directed: |
| 87 | + # Symmetrize the results. |
| 88 | + # It may be nice to be able to return these as upper-triangular. |
| 89 | + D(binary.any) << D.T |
| 90 | + if compute_predecessors: |
| 91 | + P(binary.any) << P.T |
40 | 92 |
|
41 | 93 | # Set diagonal values to 0 (this way seems fast). |
42 | 94 | # The missing values are implied to be infinity, so we set diagonals explicitly to 0. |
43 | | - mask = Vector(bool, size=n, name="mask") |
44 | | - mask << True |
45 | | - Mask = mask.diag(name="Mask") |
46 | | - D(Mask.S) << 0 |
47 | | - return D |
| 95 | + diag_mask = Vector(bool, size=n, name="diag_mask") |
| 96 | + diag_mask << True |
| 97 | + Diag_mask = diag_mask.diag(name="Diag_mask") |
| 98 | + D(Diag_mask.S) << 0 |
| 99 | + |
| 100 | + return P, D |
0 commit comments