Add {all_pairs,single_source}_bellman_ford_path_length

graphblas_algorithms/algorithms/core.py

@@ -28,7 +28,7 @@ def k_truss(G: Graph, k) -> Graph:
             S = C
 
     # Remove isolate nodes
-    indices, _ = C.reduce_rowwise(monoid.any).to_coo()
+    indices, _ = C.reduce_rowwise(monoid.any).to_coo(values=False)
     Ktruss = C[indices, indices].new()
 
     # Convert back to networkx graph with correct node ids

graphblas_algorithms/algorithms/dag.py

@@ -9,7 +9,7 @@ def descendants(G, source):
     if source not in G._key_to_id:
         raise KeyError(f"The node {source} is not in the graph")
     index = G._key_to_id[source]
-    A = G._A
+    A = G.get_property("offdiag")
     q = Vector.from_coo(index, True, size=A.nrows, name="q")
     rv = q.dup(name="descendants")
     for _ in range(A.nrows):
@@ -25,7 +25,7 @@ def ancestors(G, source):
     if source not in G._key_to_id:
         raise KeyError(f"The node {source} is not in the graph")
     index = G._key_to_id[source]
-    A = G._A
+    A = G.get_property("offdiag")
     q = Vector.from_coo(index, True, size=A.nrows, name="q")
     rv = q.dup(name="descendants")
     for _ in range(A.nrows):

graphblas_algorithms/algorithms/exceptions.py

@@ -12,3 +12,7 @@ class EmptyGraphError(GraphBlasAlgorithmException):
 
 class PointlessConcept(GraphBlasAlgorithmException):
     pass
+
+
+class Unbounded(GraphBlasAlgorithmException):
+    pass

graphblas_algorithms/algorithms/shortest_paths/__init__.py

@@ -1,2 +1,3 @@
 from .dense import *
 from .generic import *
+from .weighted import *

graphblas_algorithms/algorithms/shortest_paths/generic.py

@@ -10,7 +10,7 @@ def has_path(G, source, target):
     dst = G._key_to_id[target]
     if src == dst:
         return True
-    A = G._A
+    A = G.get_property("offdiag")
     q_src = Vector.from_coo(src, True, size=A.nrows, name="q_src")
     seen_src = q_src.dup(name="seen_src")
     q_dst = Vector.from_coo(dst, True, size=A.nrows, name="q_dst")

graphblas_algorithms/algorithms/shortest_paths/weighted.py

@@ -0,0 +1,114 @@
+import numpy as np
+from graphblas import Matrix, Vector, binary, monoid, replace, select, unary
+from graphblas.semiring import any_pair, min_plus
+
+from ..exceptions import Unbounded
+
+__all__ = [
+    "single_source_bellman_ford_path_length",
+    "bellman_ford_path_lengths",
+]
+
+
+def single_source_bellman_ford_path_length(G, source):
+    # No need for `is_weighted=` keyword, b/c this is assumed to be weighted (I think)
+    index = G._key_to_id[source]
+    # Use `offdiag` instead of `A`, b/c self-loops don't contribute to the result,
+    # and negative self-loops are easy negative cycles to avoid.
+    # We check if we hit a self-loop negative cycle at the end.
+    A, has_negative_diagonal = G.get_properties("offdiag has_negative_diagonal")
+    if A.dtype == bool:
+        # Should we upcast e.g. INT8 to INT64 as well?
+        dtype = int
+    else:
+        dtype = A.dtype
+    n = A.nrows
+    d = Vector(dtype, n, name="single_source_bellman_ford_path_length")
+    d[index] = 0
+    cur = d.dup(name="cur")
+    mask = Vector(bool, n, name="mask")
+    one = unary.one[bool]
+    for _i in range(n - 1):
+        # This is a slightly modified Bellman-Ford algorithm.
+        # `cur` is the current frontier of values that improved in the previous iteration.
+        # This means that in this iteration we drop values from `cur` that are not better.
+        cur << min_plus(cur @ A)
+
+        # Mask is True where cur not in d or cur < d
+        mask << one(cur)
+        mask(binary.second) << binary.lt(cur & d)
+
+        # Drop values from `cur` that didn't improve
+        cur(mask.V, replace) << cur
+        if cur.nvals == 0:
+            break
+        # Update `d` with values that improved
+        d(cur.S) << cur
+    else:
+        # Check for negative cycle when for loop completes without breaking
+        cur << min_plus(cur @ A)
+        mask << binary.lt(cur & d)
+        if mask.reduce(monoid.lor):
+            raise Unbounded("Negative cycle detected.")
+    if has_negative_diagonal:
+        # We removed diagonal entries above, so check if we visited one with a negative weight
+        diag = G.get_property("diag")
+        cur << select.valuelt(diag, 0)
+        if any_pair(d @ cur):
+            raise Unbounded("Negative cycle detected.")
+    return d
+
+
+def bellman_ford_path_lengths(G, nodes=None, *, expand_output=False):
+    """
+
+    Parameters
+    ----------
+    expand_output : bool, default False
+        When False, the returned Matrix has one row per node in nodes.
+        When True, the returned Matrix has the same shape as the input Matrix.
+    """
+    # Same algorithms as in `single_source_bellman_ford_path_length`, but with
+    # `Cur` as a Matrix with each row corresponding to a source node.
+    A, has_negative_diagonal = G.get_properties("offdiag has_negative_diagonal")
+    if A.dtype == bool:
+        dtype = int
+    else:
+        dtype = A.dtype
+    n = A.nrows
+    if nodes is None:
+        # TODO: `D = Vector.from_iso_value(0, n, dtype).diag()`
+        D = Vector(dtype, n, name="bellman_ford_path_lengths_vector")
+        D << 0
+        D = D.diag(name="bellman_ford_path_lengths")
+    else:
+        ids = G.list_to_ids(nodes)
+        D = Matrix.from_coo(
+            np.arange(len(ids), dtype=np.uint64), ids, 0, dtype, nrows=len(ids), ncols=n
+        )
+    Cur = D.dup(name="Cur")
+    Mask = Matrix(bool, D.nrows, D.ncols, name="Mask")
+    one = unary.one[bool]
+    for _i in range(n - 1):
+        Cur << min_plus(Cur @ A)
+        Mask << one(Cur)
+        Mask(binary.second) << binary.lt(Cur & D)
+        Cur(Mask.V, replace) << Cur
+        if Cur.nvals == 0:
+            break
+        D(Cur.S) << Cur
+    else:
+        Cur << min_plus(Cur @ A)
+        Mask << binary.lt(Cur & D)
+        if Mask.reduce_scalar(monoid.lor):
+            raise Unbounded("Negative cycle detected.")
+    if has_negative_diagonal:
+        diag = G.get_property("diag")
+        cur = select.valuelt(diag, 0)
+        if any_pair(D @ cur).nvals > 0:
+            raise Unbounded("Negative cycle detected.")
+    if nodes is not None and expand_output:
+        rv = Matrix(D.dtype, n, n, name=D.name)
+        rv[ids, :] = D
+        return rv
+    return D

graphblas_algorithms/algorithms/simple_paths.py

@@ -8,7 +8,7 @@ def is_simple_path(G, nodes):
         return False
     if len(nodes) == 1:
         return nodes[0] in G
-    A = G._A
+    A = G._A  # offdiag instead?
     if A.nvals < len(nodes) - 1:
         return False
     key_to_id = G._key_to_id

graphblas_algorithms/classes/_caching.py

@@ -1,6 +1,8 @@
-from graphblas import op
+from graphblas import Scalar, dtypes, op
 from graphblas.core import operator
 
+NONNEGATIVE_DTYPES = {dtypes.BOOL, dtypes.UINT8, dtypes.UINT16, dtypes.UINT32, dtypes.UINT64}
+
 
 def get_reduce_to_vector(key, opname, methodname):
     op_ = op.from_string(opname)
@@ -142,7 +144,7 @@ def get_reduction(G, mask=None):
                 cache[f"{keybase}+"] = cache[key]
             return cache[key]
 
-    else:
+    elif key[-1] == "+":
 
         def get_reduction(G, mask=None):
             A = G._A
@@ -170,4 +172,18 @@ def get_reduction(G, mask=None):
                 cache[f"{keybase}-"] = cache[key]
             return cache[key]
 
+    elif key.endswith("_diagonal"):
+
+        def get_reduction(G, mask=None):
+            A = G._A
+            cache = G._cache
+            if key not in cache:
+                if not G.get_property("has_self_edges"):
+                    cache[key] = Scalar(op_[A.dtype].return_type, name=key)
+                else:
+                    cache[key] = G.get_property("diag").reduce(op_).new(name=key)
+            return cache[key]
+
+    else:  # pragma: no cover (sanity)
+        raise RuntimeError()
     return get_reduction

graphblas_algorithms/classes/_utils.py

@@ -65,14 +65,14 @@ def dict_to_vector(self, d, *, size=None, dtype=None, name=None):
     return Vector.from_coo(indices, values, size=size, dtype=dtype, name=name)
 
 
-def list_to_vector(self, nodes, dtype=bool, *, size=None, name=None):
+def list_to_vector(self, nodes, dtype=None, *, values=True, size=None, name=None):
     if nodes is None:
         return None
     if size is None:
         size = len(self)
     key_to_id = self._key_to_id
     index = [key_to_id[key] for key in nodes]
-    return Vector.from_coo(index, True, size=size, dtype=dtype, name=name)
+    return Vector.from_coo(index, values, size=size, dtype=dtype, name=name)
 
 
 def list_to_mask(self, nodes, *, size=None, name="mask"):

graphblas_algorithms/classes/digraph.py

@@ -6,7 +6,7 @@
 import graphblas_algorithms as ga
 
 from . import _utils
-from ._caching import get_reduce_to_scalar, get_reduce_to_vector
+from ._caching import NONNEGATIVE_DTYPES, get_reduce_to_scalar, get_reduce_to_vector
 from .graph import Graph
 
 
@@ -371,6 +371,65 @@ def get_total_recipm(G, mask=None):
     return cache["total_recip-"]
 
 
+def has_negative_diagonal(G, mask=None):
+    A = G._A
+    cache = G._cache
+    if "has_negative_diagonal" not in cache:
+        if A.dtype in NONNEGATIVE_DTYPES or A.dtype._is_udt or cache.get("has_self_edges") is False:
+            cache["has_negative_diagonal"] = False
+        elif (
+            cache.get("has_negative_edges+") is True
+            and cache.get("has_negative_edges-") is False
+            or cache.get("has_negative_edges+") is True
+            and cache.get("min_element-", 0) >= 0
+            or cache.get("min_element+", 0) < 0
+            and cache.get("min_element+", 0) < cache.get("min_element-", 0)
+        ):
+            cache["has_negative_diagonal"] = True
+        else:
+            cache["has_negative_diagonal"] = G.get_property("min_diagonal").get(0) < 0
+    return cache["has_negative_diagonal"]
+
+
+def has_negative_edgesp(G, mask=None):
+    A = G._A
+    cache = G._cache
+    if "has_negative_edges+" not in cache:
+        if A.dtype in NONNEGATIVE_DTYPES or A.dtype._is_udt:
+            cache["has_negative_edges+"] = False
+        elif (
+            cache.get("has_negative_edges-")
+            or cache.get("min_element+", 0) < 0
+            or cache.get("min_element-", 0) < 0
+            or cache.get("min_diagonal", 0) < 0
+            or cache.get("has_negative_diagonal")
+        ):
+            cache["has_negative_edges+"] = True
+        elif cache.get("has_negative_edges-") is False:
+            cache["has_negative_edges+"] = G.get_property("min_diagonal").get(0) < 0
+        else:
+            cache["has_negative_edges+"] = G.get_property("min_element+").get(0) < 0
+    return cache["has_negative_edges+"]
+
+
+def has_negative_edgesm(G, mask=None):
+    A = G._A
+    cache = G._cache
+    if "has_negative_edges-" not in cache:
+        if A.dtype in NONNEGATIVE_DTYPES or A.dtype._is_udt:
+            cache["has_negative_edges-"] = False
+        elif (
+            cache.get("has_negative_edges+")
+            and cache.get("has_self_edges") is False
+            or cache.get("has_negative_edges+")
+            and cache.get("has_negative_diagonal") is False
+        ):
+            cache["has_negative_edges-"] = True
+        else:
+            cache["has_negative_edges-"] = G.get_property("min_element-").get(0) < 0
+    return cache["has_negative_edges-"]
+
+
 def has_self_edges(G, mask=None):
     """A.diag().nvals > 0"""
     A = G._A
@@ -388,6 +447,8 @@ def has_self_edges(G, mask=None):
                 cache["has_self_edges"] = cache["U-"].nvals < cache["U+"].nvals
             else:
                 cache["has_self_edges"] = cache["U-"].nvals + cache["L-"].nvals < A.nvals
+        elif cache.get("has_negative_diagonal") is True:
+            cache["has_self_edges"] = True
         elif "total_recip-" in cache and "total_recip+" in cache:
             cache["has_self_edges"] = cache["total_recip+"] > cache["total_recip-"]
         elif "row_degrees-" in cache and "row_degrees+" in cache:
@@ -473,9 +534,14 @@ def __missing__(self, key):
                 get_reduction = get_reduce_to_scalar(key, opname)
             else:
                 get_reduction = get_reduce_to_vector(key, opname, methodname)
-            self[key] = get_reduction
-            return get_reduction
-        raise KeyError(key)
+        elif key.endswith("_diagonal"):
+            # e.g., min_diagonal
+            opname = key[: -len("_diagonal")]
+            get_reduction = get_reduce_to_scalar(key, opname)
+        else:
+            raise KeyError(key)
+        self[key] = get_reduction
+        return get_reduction
 
 
 class DiGraph(Graph):
@@ -507,6 +573,12 @@ class DiGraph(Graph):
                     "total_degrees-",
                     "total_recip+",  # scalar; I don't like this name
                     "total_recip-",  # scalar; I don't like this name
+                    "min_diagonal",
+                    "min_element+",
+                    "min_element-",
+                    "has_negative_diagonal",
+                    "has_negative_edges-",
+                    "has_negative_edges+",
                     "has_self_edges",
                 ]
             )
@@ -528,6 +600,9 @@ class DiGraph(Graph):
             "total_degrees-": get_total_degreesm,
             "total_recip+": get_total_recipp,
             "total_recip-": get_total_recipm,
+            "has_negative_diagonal": has_negative_diagonal,
+            "has_negative_edges-": has_negative_edgesm,
+            "has_negative_edges+": has_negative_edgesp,
             "has_self_edges": has_self_edges,
         }
     )