Add floyd_warshall_predecessor_and_distance (#43)

* Add floyd_warshall_predecessor_and_distance

* Use upper triangle for undirected graphs

* Save memory when we don't need to compute predecessors

.pre-commit-config.yaml

@@ -20,12 +20,12 @@ repos:
       - id: mixed-line-ending
       - id: trailing-whitespace
   - repo: https://github.com/abravalheri/validate-pyproject
-    rev: v0.11
+    rev: v0.12.1
     hooks:
       - id: validate-pyproject
         name: Validate pyproject.toml
   - repo: https://github.com/myint/autoflake
-    rev: v2.0.0
+    rev: v2.0.1
     hooks:
       - id: autoflake
         args: [--in-place]
@@ -44,7 +44,7 @@ repos:
       - id: auto-walrus
         args: [--line-length, "100"]
   - repo: https://github.com/psf/black
-    rev: 22.12.0
+    rev: 23.1.0
     hooks:
       - id: black
         args: [--target-version=py38]

graphblas_algorithms/algorithms/shortest_paths/dense.py

@@ -1,47 +1,100 @@
-from graphblas import Matrix, Vector, binary
-from graphblas.select import offdiag
-from graphblas.semiring import any_plus
+from graphblas import Matrix, Vector, binary, indexunary, replace, select
+from graphblas.semiring import any_plus, any_second
 
-__all__ = ["floyd_warshall"]
+__all__ = ["floyd_warshall", "floyd_warshall_predecessor_and_distance"]
 
 
 def floyd_warshall(G, is_weighted=False):
+    return floyd_warshall_predecessor_and_distance(G, is_weighted, compute_predecessors=False)[1]
+
+
+def floyd_warshall_predecessor_and_distance(G, is_weighted=False, *, compute_predecessors=True):
     # By using `offdiag` instead of `G._A`, we ensure that D will not become dense.
     # Dense D may be better at times, but not including the diagonal will result in less work.
     # Typically, Floyd-Warshall algorithms sets the diagonal of D to 0 at the beginning.
     # This is unnecessary with sparse matrices, and we set the diagonal to 0 at the end.
     # We also don't iterate over index `i` if either row i or column i are empty.
-    if G.is_directed():
+    if is_directed := G.is_directed():
         A, row_degrees, column_degrees = G.get_properties("offdiag row_degrees- column_degrees-")
         nonempty_nodes = binary.pair(row_degrees & column_degrees).new(name="nonempty_nodes")
     else:
-        A, nonempty_nodes = G.get_properties("offdiag degrees-")
+        A, nonempty_nodes = G.get_properties("U- degrees-")
 
     if A.dtype == bool or not is_weighted:
         dtype = int
     else:
         dtype = A.dtype
     n = A.nrows
-    D = Matrix(dtype, nrows=n, ncols=n, name="floyd_warshall")
+    D = Matrix(dtype, nrows=n, ncols=n, name="floyd_warshall_dist")
     if is_weighted:
         D << A
     else:
         D(A.S) << 1  # Like `D << unary.one[int](A)`
     del A
-
     Row = Matrix(dtype, nrows=1, ncols=n, name="Row")
-    Col = Matrix(dtype, nrows=n, ncols=1, name="Col")
+    if is_directed:
+        Col = Matrix(dtype, nrows=n, ncols=1, name="Col")
+    else:
+        Col = None
     Outer = Matrix(dtype, nrows=n, ncols=n, name="Outer")
+    if compute_predecessors:
+        Mask = Matrix(bool, nrows=n, ncols=n, name="Mask")
+        P = indexunary.rowindex(D).new(name="floyd_warshall_pred")
+        if P.dtype == dtype:
+            P_row = Row
+        else:
+            P_row = Matrix(P.dtype, nrows=1, ncols=n, name="P_row")
+    else:
+        Mask = P = P_row = None
+
     for i in nonempty_nodes:
-        Col << D[:, [i]]
         Row << D[[i], :]
+        if is_directed:
+            Col << D[:, [i]]
+        else:
+            Row(binary.any) << D.T[[i], :]
+            Col = Row.T
         Outer << any_plus(Col @ Row)  # Like `col.outer(row, binary.plus)`
-        D(binary.min) << offdiag(Outer)
+
+        if not compute_predecessors:
+            # It is faster (approx 10%-30%) to use a mask as is done below when computing
+            # predecessors, but we choose to use less memory here by not using a mask.
+            if is_directed:
+                D(binary.min) << select.offdiag(Outer)
+            else:
+                D(binary.min) << select.triu(Outer, 1)
+        else:
+            # Update Outer to only include off-diagonal values that will update D and P.
+            if is_directed:
+                Mask << indexunary.offdiag(Outer)
+            else:
+                Mask << indexunary.triu(Outer, 1)
+            Mask(binary.second) << binary.lt(Outer & D)
+            Outer(Mask.V, replace) << Outer
+
+            # Update distances; like `D(binary.min) << offdiag(any_plus(Col @ Row))`
+            D(Outer.S) << Outer
+
+            # Broadcast predecessors in P_row to updated values
+            P_row << P[[i], :]
+            if not is_directed:
+                P_row(binary.any) << P.T[[i], :]
+                Col = P_row.T
+            P(Outer.S) << any_second(Col @ P_row)
+    del Outer, Mask, Col, Row, P_row
+
+    if not is_directed:
+        # Symmetrize the results.
+        # It may be nice to be able to return these as upper-triangular.
+        D(binary.any) << D.T
+        if compute_predecessors:
+            P(binary.any) << P.T
 
     # Set diagonal values to 0 (this way seems fast).
     # The missing values are implied to be infinity, so we set diagonals explicitly to 0.
-    mask = Vector(bool, size=n, name="mask")
-    mask << True
-    Mask = mask.diag(name="Mask")
-    D(Mask.S) << 0
-    return D
+    diag_mask = Vector(bool, size=n, name="diag_mask")
+    diag_mask << True
+    Diag_mask = diag_mask.diag(name="Diag_mask")
+    D(Diag_mask.S) << 0
+
+    return P, D

graphblas_algorithms/classes/_utils.py

@@ -109,26 +109,27 @@ def set_to_vector(self, nodes, dtype=bool, *, ignore_extra=False, size=None, nam
     return Vector.from_coo(index, True, size=size, dtype=dtype, name=name)
 
 
-def vector_to_dict(self, v, *, mask=None, fillvalue=None):
+def vector_to_dict(self, v, *, mask=None, fill_value=None):
     if mask is not None:
-        if fillvalue is not None and v.nvals < mask.parent.nvals:
-            v(mask, binary.first) << fillvalue
-    elif fillvalue is not None and v.nvals < v.size:
-        v(mask=~v.S) << fillvalue
+        if fill_value is not None and v.nvals < mask.parent.nvals:
+            v(mask, binary.first) << fill_value
+    elif fill_value is not None and v.nvals < v.size:
+        v(mask=~v.S) << fill_value
     id_to_key = self.id_to_key
     return {id_to_key[index]: value for index, value in zip(*v.to_coo(sort=False))}
 
 
-def vector_to_nodemap(self, v, *, mask=None, fillvalue=None):
+def vector_to_nodemap(self, v, *, mask=None, fill_value=None, values_are_keys=False):
     from .nodemap import NodeMap
 
     if mask is not None:
-        if fillvalue is not None and v.nvals < mask.parent.nvals:
-            v(mask, binary.first) << fillvalue
-    elif fillvalue is not None and v.nvals < v.size:
-        v(mask=~v.S) << fillvalue
+        if fill_value is not None and v.nvals < mask.parent.nvals:
+            v(mask, binary.first) << fill_value
+        fill_value = None
 
-    rv = NodeMap(v, key_to_id=self._key_to_id)
+    rv = NodeMap(
+        v, fill_value=fill_value, values_are_keys=values_are_keys, key_to_id=self._key_to_id
+    )
     rv._id_to_key = self._id_to_key
     return rv
 
@@ -147,7 +148,25 @@ def vector_to_set(self, v):
     return {id_to_key[index] for index in indices}
 
 
-def matrix_to_dicts(self, A, *, use_row_index=False, use_column_index=False):
+def matrix_to_nodenodemap(self, A, *, fill_value=None, values_are_keys=False):
+    from .nodemap import NodeNodeMap
+
+    rv = NodeNodeMap(
+        A, fill_value=fill_value, values_are_keys=values_are_keys, key_to_id=self._key_to_id
+    )
+    rv._id_to_key = self._id_to_key
+    return rv
+
+
+def matrix_to_vectornodemap(self, A):
+    from .nodemap import VectorNodeMap
+
+    rv = VectorNodeMap(A, key_to_id=self._key_to_id)
+    rv._id_to_key = self._id_to_key
+    return rv
+
+
+def matrix_to_dicts(self, A, *, use_row_index=False, use_column_index=False, values_are_keys=False):
     """Convert a Matrix to a dict of dicts of the form ``{row: {col: val}}``
 
     Use ``use_row_index=True`` to return the row index as keys in the dict,
@@ -167,6 +186,8 @@ def matrix_to_dicts(self, A, *, use_row_index=False, use_column_index=False):
     indptr = d["indptr"]
     values = d["values"].tolist()
     id_to_key = self.id_to_key
+    if values_are_keys:
+        values = [id_to_key[val] for val in values]
     it = zip(rows, np.lib.stride_tricks.sliding_window_view(indptr, 2).tolist())
     if use_row_index and use_column_index:
         return {

graphblas_algorithms/classes/digraph.py

@@ -569,6 +569,8 @@ def __init__(self, incoming_graph_data=None, *, key_to_id=None, **attr):
     list_to_mask = _utils.list_to_mask
     list_to_ids = _utils.list_to_ids
     matrix_to_dicts = _utils.matrix_to_dicts
+    matrix_to_nodenodemap = _utils.matrix_to_nodenodemap
+    matrix_to_vectornodemap = _utils.matrix_to_vectornodemap
     set_to_vector = _utils.set_to_vector
     to_networkx = _utils.to_networkx
     vector_to_dict = _utils.vector_to_dict

graphblas_algorithms/classes/graph.py

@@ -275,6 +275,8 @@ def __init__(self, incoming_graph_data=None, *, key_to_id=None, **attr):
     list_to_ids = _utils.list_to_ids
     list_to_keys = _utils.list_to_keys
     matrix_to_dicts = _utils.matrix_to_dicts
+    matrix_to_nodenodemap = _utils.matrix_to_nodenodemap
+    matrix_to_vectornodemap = _utils.matrix_to_vectornodemap
     set_to_vector = _utils.set_to_vector
     to_networkx = _utils.to_networkx
     vector_to_dict = _utils.vector_to_dict