graph/path: use queue-based Bellman-Ford algorithm

See Sedgewick and Wayne, Algorithms 4th Edition from p672 onward.

graph/path/bellman_ford_moore.go

@@ -4,7 +4,10 @@
 
 package path
 
-import "gonum.org/v1/gonum/graph"
+import (
+	"gonum.org/v1/gonum/graph"
+	"gonum.org/v1/gonum/graph/internal/linear"
+)
 
 // BellmanFordFrom returns a shortest-path tree for a shortest path from u to all nodes in
 // the graph g, or false indicating that a negative cycle exists in the graph. If the graph
@@ -27,33 +30,24 @@ func BellmanFordFrom(u graph.Node, g graph.Graph) (path Shortest, ok bool) {
 	path = newShortestFrom(u, nodes)
 	path.dist[path.indexOf[u.ID()]] = 0
 
+	// Queue to keep track which nodes need to be relaxed.
+	// Only nodes whose vertex distance changed in the previous iterations need to be relaxed again.
+	queue := newBellmanFordQueue(path.indexOf)
+	queue.enqueue(u)
+
+	// The maximum of edges in a graph is |V| * (|V|-1) which is also the worst case complexity.
+	// If the queue-loop has more iterations than the amount of maximum edges
+	// it indicates that we have a negative cycle.
+	maxEdges := len(nodes) * (len(nodes) - 1)
+	var loops int
+
 	// TODO(kortschak): Consider adding further optimisations
 	// from http://arxiv.org/abs/1111.5414.
-	for i := 1; i < len(nodes); i++ {
-		changed := false
-		for j, u := range nodes {
-			uid := u.ID()
-			for _, v := range graph.NodesOf(g.From(uid)) {
-				vid := v.ID()
-				k := path.indexOf[vid]
-				w, ok := weight(uid, vid)
-				if !ok {
-					panic("bellman-ford: unexpected invalid weight")
-				}
-				joint := path.dist[j] + w
-				if joint < path.dist[k] {
-					path.set(k, joint, j)
-					changed = true
-				}
-			}
-		}
-		if !changed {
-			break
-		}
-	}
-
-	for j, u := range nodes {
+	for queue.len() != 0 {
+		u := queue.dequeue()
 		uid := u.ID()
+		j := path.indexOf[uid]
+
 		for _, v := range graph.NodesOf(g.From(uid)) {
 			vid := v.ID()
 			k := path.indexOf[vid]
@@ -61,12 +55,66 @@ func BellmanFordFrom(u graph.Node, g graph.Graph) (path Shortest, ok bool) {
 			if !ok {
 				panic("bellman-ford: unexpected invalid weight")
 			}
-			if path.dist[j]+w < path.dist[k] {
-				path.hasNegativeCycle = true
-				return path, false
+
+			joint := path.dist[j] + w
+			if joint < path.dist[k] {
+				path.set(k, joint, j)
+
+				if !queue.has(vid) {
+					queue.enqueue(v)
+				}
 			}
 		}
+
+		if loops > maxEdges {
+			path.hasNegativeCycle = true
+			return path, false
+		}
+		loops++
 	}
 
 	return path, true
 }
+
+// bellmanFordQueue is a queue for the Queue-based Bellman-Ford algorithm.
+type bellmanFordQueue struct {
+	// queue holds the nodes which need to be relaxed.
+	queue linear.NodeQueue
+
+	// onQueue keeps track whether a node is on the queue or not.
+	onQueue []bool
+
+	// indexOf contains a mapping holding the id of a node with its index in the onQueue array.
+	indexOf map[int64]int
+}
+
+// enqueue adds a node to the bellmanFordQueue.
+func (q *bellmanFordQueue) enqueue(n graph.Node) {
+	i := q.indexOf[n.ID()]
+	if q.onQueue[i] {
+		panic("bellman-ford: already queued")
+	}
+	q.onQueue[i] = true
+	q.queue.Enqueue(n)
+}
+
+// dequeue returns the first value of the bellmanFordQueue.
+func (q *bellmanFordQueue) dequeue() graph.Node {
+	n := q.queue.Dequeue()
+	q.onQueue[q.indexOf[n.ID()]] = false
+	return n
+}
+
+// len returns the number of nodes in the bellmanFordQueue.
+func (q *bellmanFordQueue) len() int { return q.queue.Len() }
+
+// has returns whether a node with the given id is in the queue.
+func (q bellmanFordQueue) has(id int64) bool { return q.onQueue[q.indexOf[id]] }
+
+// newBellmanFordQueue creates a new bellmanFordQueue.
+func newBellmanFordQueue(indexOf map[int64]int) bellmanFordQueue {
+	return bellmanFordQueue{
+		onQueue: make([]bool, len(indexOf)),
+		indexOf: indexOf,
+	}
+}

graph/path/bench_test.go

@@ -139,3 +139,46 @@ func BenchmarkAStarUndirectedmallWorld_100_5_20_2_Heur(b *testing.B) {
 	}
 	benchmarkAStarHeuristic(b, nswUndirected_100_5_20_2, h)
 }
+
+var (
+	gnpDirected_500_tenth  = gnpDirected(500, 0.1)
+	gnpDirected_1000_tenth = gnpDirected(1000, 0.1)
+	gnpDirected_2000_tenth = gnpDirected(2000, 0.1)
+	gnpDirected_500_half   = gnpDirected(500, 0.5)
+	gnpDirected_1000_half  = gnpDirected(1000, 0.5)
+	gnpDirected_2000_half  = gnpDirected(2000, 0.5)
+	gnpDirected_500_full   = gnpDirected(500, 1)
+	gnpDirected_1000_full  = gnpDirected(1000, 1)
+	gnpDirected_2000_full  = gnpDirected(2000, 1)
+)
+
+func gnpDirected(n int, p float64) graph.Directed {
+	g := simple.NewDirectedGraph()
+	gen.Gnp(g, n, p, nil)
+	return g
+}
+
+func BenchmarkBellmanFordFrom(b *testing.B) {
+	benchmarks := []struct {
+		name  string
+		graph graph.Directed
+	}{
+		{"500 tenth", gnpDirected_500_tenth},
+		{"1000 tenth", gnpDirected_1000_tenth},
+		{"2000 tenth", gnpDirected_2000_tenth},
+		{"500 half", gnpDirected_500_half},
+		{"1000 half", gnpDirected_1000_half},
+		{"2000 half", gnpDirected_2000_half},
+		{"500 full", gnpDirected_500_full},
+		{"1000 full", gnpDirected_1000_full},
+		{"2000 full", gnpDirected_2000_full},
+	}
+
+	for _, bm := range benchmarks {
+		b.Run(bm.name, func(b *testing.B) {
+			for i := 0; i < b.N; i++ {
+				BellmanFordFrom(bm.graph.Node(0), bm.graph)
+			}
+		})
+	}
+}