aboutsummaryrefslogtreecommitdiffstats
path: root/vendor/github.com/soniakeys/graph/sssp.go
diff options
context:
space:
mode:
Diffstat (limited to 'vendor/github.com/soniakeys/graph/sssp.go')
-rw-r--r--vendor/github.com/soniakeys/graph/sssp.go881
1 files changed, 881 insertions, 0 deletions
diff --git a/vendor/github.com/soniakeys/graph/sssp.go b/vendor/github.com/soniakeys/graph/sssp.go
new file mode 100644
index 0000000..32cc192
--- /dev/null
+++ b/vendor/github.com/soniakeys/graph/sssp.go
@@ -0,0 +1,881 @@
+// Copyright 2013 Sonia Keys
+// License MIT: http://opensource.org/licenses/MIT
+
+package graph
+
+import (
+ "container/heap"
+ "fmt"
+ "math"
+)
+
+// rNode holds data for a "reached" node
+type rNode struct {
+ nx NI
+ state int8 // state constants defined below
+ f float64 // "g+h", path dist + heuristic estimate
+ fx int // heap.Fix index
+}
+
+// for rNode.state
+const (
+ unreached = 0
+ reached = 1
+ open = 1
+ closed = 2
+)
+
+type openHeap []*rNode
+
+// A Heuristic is defined on a specific end node. The function
+// returns an estimate of the path distance from node argument
+// "from" to the end node. Two subclasses of heuristics are "admissible"
+// and "monotonic."
+//
+// Admissible means the value returned is guaranteed to be less than or
+// equal to the actual shortest path distance from the node to end.
+//
+// An admissible estimate may further be monotonic.
+// Monotonic means that for any neighboring nodes A and B with half arc aB
+// leading from A to B, and for heuristic h defined on some end node, then
+// h(A) <= aB.ArcWeight + h(B).
+//
+// See AStarA for additional notes on implementing heuristic functions for
+// AStar search methods.
+type Heuristic func(from NI) float64
+
+// Admissible returns true if heuristic h is admissible on graph g relative to
+// the given end node.
+//
+// If h is inadmissible, the string result describes a counter example.
+func (h Heuristic) Admissible(g LabeledAdjacencyList, w WeightFunc, end NI) (bool, string) {
+ // invert graph
+ inv := make(LabeledAdjacencyList, len(g))
+ for from, nbs := range g {
+ for _, nb := range nbs {
+ inv[nb.To] = append(inv[nb.To],
+ Half{To: NI(from), Label: nb.Label})
+ }
+ }
+ // run dijkstra
+ // Dijkstra.AllPaths takes a start node but after inverting the graph
+ // argument end now represents the start node of the inverted graph.
+ f, dist, _ := inv.Dijkstra(end, -1, w)
+ // compare h to found shortest paths
+ for n := range inv {
+ if f.Paths[n].Len == 0 {
+ continue // no path, any heuristic estimate is fine.
+ }
+ if !(h(NI(n)) <= dist[n]) {
+ return false, fmt.Sprintf("h(%d) = %g, "+
+ "required to be <= found shortest path (%g)",
+ n, h(NI(n)), dist[n])
+ }
+ }
+ return true, ""
+}
+
+// Monotonic returns true if heuristic h is monotonic on weighted graph g.
+//
+// If h is non-monotonic, the string result describes a counter example.
+func (h Heuristic) Monotonic(g LabeledAdjacencyList, w WeightFunc) (bool, string) {
+ // precompute
+ hv := make([]float64, len(g))
+ for n := range g {
+ hv[n] = h(NI(n))
+ }
+ // iterate over all edges
+ for from, nbs := range g {
+ for _, nb := range nbs {
+ arcWeight := w(nb.Label)
+ if !(hv[from] <= arcWeight+hv[nb.To]) {
+ return false, fmt.Sprintf("h(%d) = %g, "+
+ "required to be <= arc weight + h(%d) (= %g + %g = %g)",
+ from, hv[from],
+ nb.To, arcWeight, hv[nb.To], arcWeight+hv[nb.To])
+ }
+ }
+ }
+ return true, ""
+}
+
+// AStarA finds a path between two nodes.
+//
+// AStarA implements both algorithm A and algorithm A*. The difference in the
+// two algorithms is strictly in the heuristic estimate returned by argument h.
+// If h is an "admissible" heuristic estimate, then the algorithm is termed A*,
+// otherwise it is algorithm A.
+//
+// Like Dijkstra's algorithm, AStarA with an admissible heuristic finds the
+// shortest path between start and end. AStarA generally runs faster than
+// Dijkstra though, by using the heuristic distance estimate.
+//
+// AStarA with an inadmissible heuristic becomes algorithm A. Algorithm A
+// will find a path, but it is not guaranteed to be the shortest path.
+// The heuristic still guides the search however, so a nearly admissible
+// heuristic is likely to find a very good path, if not the best. Quality
+// of the path returned degrades gracefully with the quality of the heuristic.
+//
+// The heuristic function h should ideally be fairly inexpensive. AStarA
+// may call it more than once for the same node, especially as graph density
+// increases. In some cases it may be worth the effort to memoize or
+// precompute values.
+//
+// Argument g is the graph to be searched, with arc weights returned by w.
+// As usual for AStar, arc weights must be non-negative.
+// Graphs may be directed or undirected.
+//
+// If AStarA finds a path it returns a FromList encoding the path, the arc
+// labels for path nodes, the total path distance, and ok = true.
+// Otherwise it returns ok = false.
+func (g LabeledAdjacencyList) AStarA(w WeightFunc, start, end NI, h Heuristic) (f FromList, labels []LI, dist float64, ok bool) {
+ // NOTE: AStarM is largely duplicate code.
+
+ f = NewFromList(len(g))
+ labels = make([]LI, len(g))
+ d := make([]float64, len(g))
+ r := make([]rNode, len(g))
+ for i := range r {
+ r[i].nx = NI(i)
+ }
+ // start node is reached initially
+ cr := &r[start]
+ cr.state = reached
+ cr.f = h(start) // total path estimate is estimate from start
+ rp := f.Paths
+ rp[start] = PathEnd{Len: 1, From: -1} // path length at start is 1 node
+ // oh is a heap of nodes "open" for exploration. nodes go on the heap
+ // when they get an initial or new "g" path distance, and therefore a
+ // new "f" which serves as priority for exploration.
+ oh := openHeap{cr}
+ for len(oh) > 0 {
+ bestPath := heap.Pop(&oh).(*rNode)
+ bestNode := bestPath.nx
+ if bestNode == end {
+ return f, labels, d[end], true
+ }
+ bp := &rp[bestNode]
+ nextLen := bp.Len + 1
+ for _, nb := range g[bestNode] {
+ alt := &r[nb.To]
+ ap := &rp[alt.nx]
+ // "g" path distance from start
+ g := d[bestNode] + w(nb.Label)
+ if alt.state == reached {
+ if g > d[nb.To] {
+ // candidate path to nb is longer than some alternate path
+ continue
+ }
+ if g == d[nb.To] && nextLen >= ap.Len {
+ // candidate path has identical length of some alternate
+ // path but it takes no fewer hops.
+ continue
+ }
+ // cool, we found a better way to get to this node.
+ // record new path data for this node and
+ // update alt with new data and make sure it's on the heap.
+ *ap = PathEnd{From: bestNode, Len: nextLen}
+ labels[nb.To] = nb.Label
+ d[nb.To] = g
+ alt.f = g + h(nb.To)
+ if alt.fx < 0 {
+ heap.Push(&oh, alt)
+ } else {
+ heap.Fix(&oh, alt.fx)
+ }
+ } else {
+ // bestNode being reached for the first time.
+ *ap = PathEnd{From: bestNode, Len: nextLen}
+ labels[nb.To] = nb.Label
+ d[nb.To] = g
+ alt.f = g + h(nb.To)
+ alt.state = reached
+ heap.Push(&oh, alt) // and it's now open for exploration
+ }
+ }
+ }
+ return // no path
+}
+
+// AStarAPath finds a shortest path using the AStarA algorithm.
+//
+// This is a convenience method with a simpler result than the AStarA method.
+// See documentation on the AStarA method.
+//
+// If a path is found, the non-nil node path is returned with the total path
+// distance. Otherwise the returned path will be nil.
+func (g LabeledAdjacencyList) AStarAPath(start, end NI, h Heuristic, w WeightFunc) ([]NI, float64) {
+ f, _, d, _ := g.AStarA(w, start, end, h)
+ return f.PathTo(end, nil), d
+}
+
+// AStarM is AStarA optimized for monotonic heuristic estimates.
+//
+// Note that this function requires a monotonic heuristic. Results will
+// not be meaningful if argument h is non-monotonic.
+//
+// See AStarA for general usage. See Heuristic for notes on monotonicity.
+func (g LabeledAdjacencyList) AStarM(w WeightFunc, start, end NI, h Heuristic) (f FromList, labels []LI, dist float64, ok bool) {
+ // NOTE: AStarM is largely code duplicated from AStarA.
+ // Differences are noted in comments in this method.
+
+ f = NewFromList(len(g))
+ labels = make([]LI, len(g))
+ d := make([]float64, len(g))
+ r := make([]rNode, len(g))
+ for i := range r {
+ r[i].nx = NI(i)
+ }
+ cr := &r[start]
+
+ // difference from AStarA:
+ // instead of a bit to mark a reached node, there are two states,
+ // open and closed. open marks nodes "open" for exploration.
+ // nodes are marked open as they are reached, then marked
+ // closed as they are found to be on the best path.
+ cr.state = open
+
+ cr.f = h(start)
+ rp := f.Paths
+ rp[start] = PathEnd{Len: 1, From: -1}
+ oh := openHeap{cr}
+ for len(oh) > 0 {
+ bestPath := heap.Pop(&oh).(*rNode)
+ bestNode := bestPath.nx
+ if bestNode == end {
+ return f, labels, d[end], true
+ }
+
+ // difference from AStarA:
+ // move nodes to closed list as they are found to be best so far.
+ bestPath.state = closed
+
+ bp := &rp[bestNode]
+ nextLen := bp.Len + 1
+ for _, nb := range g[bestNode] {
+ alt := &r[nb.To]
+
+ // difference from AStarA:
+ // Monotonicity means that f cannot be improved.
+ if alt.state == closed {
+ continue
+ }
+
+ ap := &rp[alt.nx]
+ g := d[bestNode] + w(nb.Label)
+
+ // difference from AStarA:
+ // test for open state, not just reached
+ if alt.state == open {
+
+ if g > d[nb.To] {
+ continue
+ }
+ if g == d[nb.To] && nextLen >= ap.Len {
+ continue
+ }
+ *ap = PathEnd{From: bestNode, Len: nextLen}
+ labels[nb.To] = nb.Label
+ d[nb.To] = g
+ alt.f = g + h(nb.To)
+
+ // difference from AStarA:
+ // we know alt was on the heap because we found it marked open
+ heap.Fix(&oh, alt.fx)
+ } else {
+ *ap = PathEnd{From: bestNode, Len: nextLen}
+ labels[nb.To] = nb.Label
+ d[nb.To] = g
+ alt.f = g + h(nb.To)
+
+ // difference from AStarA:
+ // nodes are opened when first reached
+ alt.state = open
+ heap.Push(&oh, alt)
+ }
+ }
+ }
+ return
+}
+
+// AStarMPath finds a shortest path using the AStarM algorithm.
+//
+// This is a convenience method with a simpler result than the AStarM method.
+// See documentation on the AStarM and AStarA methods.
+//
+// If a path is found, the non-nil node path is returned with the total path
+// distance. Otherwise the returned path will be nil.
+func (g LabeledAdjacencyList) AStarMPath(start, end NI, h Heuristic, w WeightFunc) ([]NI, float64) {
+ f, _, d, _ := g.AStarM(w, start, end, h)
+ return f.PathTo(end, nil), d
+}
+
+// implement container/heap
+func (h openHeap) Len() int { return len(h) }
+func (h openHeap) Less(i, j int) bool { return h[i].f < h[j].f }
+func (h openHeap) Swap(i, j int) {
+ h[i], h[j] = h[j], h[i]
+ h[i].fx = i
+ h[j].fx = j
+}
+func (p *openHeap) Push(x interface{}) {
+ h := *p
+ fx := len(h)
+ h = append(h, x.(*rNode))
+ h[fx].fx = fx
+ *p = h
+}
+
+func (p *openHeap) Pop() interface{} {
+ h := *p
+ last := len(h) - 1
+ *p = h[:last]
+ h[last].fx = -1
+ return h[last]
+}
+
+// BellmanFord finds shortest paths from a start node in a weighted directed
+// graph using the Bellman-Ford-Moore algorithm.
+//
+// WeightFunc w must translate arc labels to arc weights.
+// Negative arc weights are allowed but not negative cycles.
+// Loops and parallel arcs are allowed.
+//
+// If the algorithm completes without encountering a negative cycle the method
+// returns shortest paths encoded in a FromList, path distances indexed by
+// node, and return value end = -1.
+//
+// If it encounters a negative cycle reachable from start it returns end >= 0.
+// In this case the cycle can be obtained by calling f.BellmanFordCycle(end).
+//
+// Negative cycles are only detected when reachable from start. A negative
+// cycle not reachable from start will not prevent the algorithm from finding
+// shortest paths from start.
+//
+// See also NegativeCycle to find a cycle anywhere in the graph, and see
+// HasNegativeCycle for lighter-weight negative cycle detection,
+func (g LabeledDirected) BellmanFord(w WeightFunc, start NI) (f FromList, dist []float64, end NI) {
+ a := g.LabeledAdjacencyList
+ f = NewFromList(len(a))
+ dist = make([]float64, len(a))
+ inf := math.Inf(1)
+ for i := range dist {
+ dist[i] = inf
+ }
+ rp := f.Paths
+ rp[start] = PathEnd{Len: 1, From: -1}
+ dist[start] = 0
+ for _ = range a[1:] {
+ imp := false
+ for from, nbs := range a {
+ fp := &rp[from]
+ d1 := dist[from]
+ for _, nb := range nbs {
+ d2 := d1 + w(nb.Label)
+ to := &rp[nb.To]
+ // TODO improve to break ties
+ if fp.Len > 0 && d2 < dist[nb.To] {
+ *to = PathEnd{From: NI(from), Len: fp.Len + 1}
+ dist[nb.To] = d2
+ imp = true
+ }
+ }
+ }
+ if !imp {
+ break
+ }
+ }
+ for from, nbs := range a {
+ d1 := dist[from]
+ for _, nb := range nbs {
+ if d1+w(nb.Label) < dist[nb.To] {
+ // return nb as end of a path with negative cycle at root
+ return f, dist, NI(from)
+ }
+ }
+ }
+ return f, dist, -1
+}
+
+// BellmanFordCycle decodes a negative cycle detected by BellmanFord.
+//
+// Receiver f and argument end must be results returned from BellmanFord.
+func (f FromList) BellmanFordCycle(end NI) (c []NI) {
+ p := f.Paths
+ var b Bits
+ for b.Bit(end) == 0 {
+ b.SetBit(end, 1)
+ end = p[end].From
+ }
+ for b.Bit(end) == 1 {
+ c = append(c, end)
+ b.SetBit(end, 0)
+ end = p[end].From
+ }
+ for i, j := 0, len(c)-1; i < j; i, j = i+1, j-1 {
+ c[i], c[j] = c[j], c[i]
+ }
+ return
+}
+
+// HasNegativeCycle returns true if the graph contains any negative cycle.
+//
+// HasNegativeCycle uses a Bellman-Ford-like algorithm, but finds negative
+// cycles anywhere in the graph. Also path information is not computed,
+// reducing memory use somewhat compared to BellmanFord.
+//
+// See also NegativeCycle to obtain the cycle, and see BellmanFord for
+// single source shortest path searches.
+func (g LabeledDirected) HasNegativeCycle(w WeightFunc) bool {
+ a := g.LabeledAdjacencyList
+ dist := make([]float64, len(a))
+ for _ = range a[1:] {
+ imp := false
+ for from, nbs := range a {
+ d1 := dist[from]
+ for _, nb := range nbs {
+ d2 := d1 + w(nb.Label)
+ if d2 < dist[nb.To] {
+ dist[nb.To] = d2
+ imp = true
+ }
+ }
+ }
+ if !imp {
+ break
+ }
+ }
+ for from, nbs := range a {
+ d1 := dist[from]
+ for _, nb := range nbs {
+ if d1+w(nb.Label) < dist[nb.To] {
+ return true // negative cycle
+ }
+ }
+ }
+ return false
+}
+
+// NegativeCycle finds a negative cycle if one exists.
+//
+// NegativeCycle uses a Bellman-Ford-like algorithm, but finds negative
+// cycles anywhere in the graph. If a negative cycle exists, one will be
+// returned. The result is nil if no negative cycle exists.
+//
+// See also HasNegativeCycle for lighter-weight cycle detection, and see
+// BellmanFord for single source shortest paths.
+func (g LabeledDirected) NegativeCycle(w WeightFunc) (c []NI) {
+ a := g.LabeledAdjacencyList
+ f := NewFromList(len(a))
+ p := f.Paths
+ for n := range p {
+ p[n] = PathEnd{From: -1, Len: 1}
+ }
+ dist := make([]float64, len(a))
+ for _ = range a {
+ imp := false
+ for from, nbs := range a {
+ fp := &p[from]
+ d1 := dist[from]
+ for _, nb := range nbs {
+ d2 := d1 + w(nb.Label)
+ to := &p[nb.To]
+ if fp.Len > 0 && d2 < dist[nb.To] {
+ *to = PathEnd{From: NI(from), Len: fp.Len + 1}
+ dist[nb.To] = d2
+ imp = true
+ }
+ }
+ }
+ if !imp {
+ return nil
+ }
+ }
+ var vis Bits
+a:
+ for n := range a {
+ end := NI(n)
+ var b Bits
+ for b.Bit(end) == 0 {
+ if vis.Bit(end) == 1 {
+ continue a
+ }
+ vis.SetBit(end, 1)
+ b.SetBit(end, 1)
+ end = p[end].From
+ if end < 0 {
+ continue a
+ }
+ }
+ for b.Bit(end) == 1 {
+ c = append(c, end)
+ b.SetBit(end, 0)
+ end = p[end].From
+ }
+ for i, j := 0, len(c)-1; i < j; i, j = i+1, j-1 {
+ c[i], c[j] = c[j], c[i]
+ }
+ return c
+ }
+ return nil // no negative cycle
+}
+
+// A NodeVisitor is an argument to some graph traversal methods.
+//
+// Graph traversal methods call the visitor function for each node visited.
+// Argument n is the node being visited.
+type NodeVisitor func(n NI)
+
+// An OkNodeVisitor function is an argument to some graph traversal methods.
+//
+// Graph traversal methods call the visitor function for each node visited.
+// The argument n is the node being visited. If the visitor function
+// returns true, the traversal will continue. If the visitor function
+// returns false, the traversal will terminate immediately.
+type OkNodeVisitor func(n NI) (ok bool)
+
+// BreadthFirst2 traverses a graph breadth first using a direction
+// optimizing algorithm.
+//
+// The code is experimental and currently seems no faster than the
+// conventional breadth first code.
+//
+// Use AdjacencyList.BreadthFirst instead.
+func BreadthFirst2(g, tr AdjacencyList, ma int, start NI, f *FromList, v OkNodeVisitor) int {
+ if tr == nil {
+ var d Directed
+ d, ma = Directed{g}.Transpose()
+ tr = d.AdjacencyList
+ }
+ switch {
+ case f == nil:
+ e := NewFromList(len(g))
+ f = &e
+ case f.Paths == nil:
+ *f = NewFromList(len(g))
+ }
+ if ma <= 0 {
+ ma = g.ArcSize()
+ }
+ rp := f.Paths
+ level := 1
+ rp[start] = PathEnd{Len: level, From: -1}
+ if !v(start) {
+ f.MaxLen = level
+ return -1
+ }
+ nReached := 1 // accumulated for a return value
+ // the frontier consists of nodes all at the same level
+ frontier := []NI{start}
+ mf := len(g[start]) // number of arcs leading out from frontier
+ ctb := ma / 10 // threshold change from top-down to bottom-up
+ k14 := 14 * ma / len(g) // 14 * mean degree
+ cbt := len(g) / k14 // threshold change from bottom-up to top-down
+ // var fBits, nextb big.Int
+ fBits := make([]bool, len(g))
+ nextb := make([]bool, len(g))
+ zBits := make([]bool, len(g))
+ for {
+ // top down step
+ level++
+ var next []NI
+ for _, n := range frontier {
+ for _, nb := range g[n] {
+ if rp[nb].Len == 0 {
+ rp[nb] = PathEnd{From: n, Len: level}
+ if !v(nb) {
+ f.MaxLen = level
+ return -1
+ }
+ next = append(next, nb)
+ nReached++
+ }
+ }
+ }
+ if len(next) == 0 {
+ break
+ }
+ frontier = next
+ if mf > ctb {
+ // switch to bottom up!
+ } else {
+ // stick with top down
+ continue
+ }
+ // convert frontier representation
+ nf := 0 // number of vertices on the frontier
+ for _, n := range frontier {
+ // fBits.SetBit(&fBits, n, 1)
+ fBits[n] = true
+ nf++
+ }
+ bottomUpLoop:
+ level++
+ nNext := 0
+ for n := range tr {
+ if rp[n].Len == 0 {
+ for _, nb := range tr[n] {
+ // if fBits.Bit(nb) == 1 {
+ if fBits[nb] {
+ rp[n] = PathEnd{From: nb, Len: level}
+ if !v(nb) {
+ f.MaxLen = level
+ return -1
+ }
+ // nextb.SetBit(&nextb, n, 1)
+ nextb[n] = true
+ nReached++
+ nNext++
+ break
+ }
+ }
+ }
+ }
+ if nNext == 0 {
+ break
+ }
+ fBits, nextb = nextb, fBits
+ // nextb.SetInt64(0)
+ copy(nextb, zBits)
+ nf = nNext
+ if nf < cbt {
+ // switch back to top down!
+ } else {
+ // stick with bottom up
+ goto bottomUpLoop
+ }
+ // convert frontier representation
+ mf = 0
+ frontier = frontier[:0]
+ for n := range g {
+ // if fBits.Bit(n) == 1 {
+ if fBits[n] {
+ frontier = append(frontier, NI(n))
+ mf += len(g[n])
+ fBits[n] = false
+ }
+ }
+ // fBits.SetInt64(0)
+ }
+ f.MaxLen = level - 1
+ return nReached
+}
+
+// DAGMinDistPath finds a single shortest path.
+//
+// Shortest means minimum sum of arc weights.
+//
+// Returned is the path and distance as returned by FromList.PathTo.
+//
+// This is a convenience method. See DAGOptimalPaths for more options.
+func (g LabeledDirected) DAGMinDistPath(start, end NI, w WeightFunc) ([]NI, float64, error) {
+ return g.dagPath(start, end, w, false)
+}
+
+// DAGMaxDistPath finds a single longest path.
+//
+// Longest means maximum sum of arc weights.
+//
+// Returned is the path and distance as returned by FromList.PathTo.
+//
+// This is a convenience method. See DAGOptimalPaths for more options.
+func (g LabeledDirected) DAGMaxDistPath(start, end NI, w WeightFunc) ([]NI, float64, error) {
+ return g.dagPath(start, end, w, true)
+}
+
+func (g LabeledDirected) dagPath(start, end NI, w WeightFunc, longest bool) ([]NI, float64, error) {
+ o, _ := g.Topological()
+ if o == nil {
+ return nil, 0, fmt.Errorf("not a DAG")
+ }
+ f, dist, _ := g.DAGOptimalPaths(start, end, o, w, longest)
+ if f.Paths[end].Len == 0 {
+ return nil, 0, fmt.Errorf("no path from %d to %d", start, end)
+ }
+ return f.PathTo(end, nil), dist[end], nil
+}
+
+// DAGOptimalPaths finds either longest or shortest distance paths in a
+// directed acyclic graph.
+//
+// Path distance is the sum of arc weights on the path.
+// Negative arc weights are allowed.
+// Where multiple paths exist with the same distance, the path length
+// (number of nodes) is used as a tie breaker.
+//
+// Receiver g must be a directed acyclic graph. Argument o must be either nil
+// or a topological ordering of g. If nil, a topologcal ordering is
+// computed internally. If longest is true, an optimal path is a longest
+// distance path. Otherwise it is a shortest distance path.
+//
+// Argument start is the start node for paths, end is the end node. If end
+// is a valid node number, the method returns as soon as the optimal path
+// to end is found. If end is -1, all optimal paths from start are found.
+//
+// Paths and path distances are encoded in the returned FromList and dist
+// slice. The number of nodes reached is returned as nReached.
+func (g LabeledDirected) DAGOptimalPaths(start, end NI, ordering []NI, w WeightFunc, longest bool) (f FromList, dist []float64, nReached int) {
+ a := g.LabeledAdjacencyList
+ f = NewFromList(len(a))
+ dist = make([]float64, len(a))
+ if ordering == nil {
+ ordering, _ = g.Topological()
+ }
+ // search ordering for start
+ o := 0
+ for ordering[o] != start {
+ o++
+ }
+ var fBetter func(cand, ext float64) bool
+ var iBetter func(cand, ext int) bool
+ if longest {
+ fBetter = func(cand, ext float64) bool { return cand > ext }
+ iBetter = func(cand, ext int) bool { return cand > ext }
+ } else {
+ fBetter = func(cand, ext float64) bool { return cand < ext }
+ iBetter = func(cand, ext int) bool { return cand < ext }
+ }
+ p := f.Paths
+ p[start] = PathEnd{From: -1, Len: 1}
+ f.MaxLen = 1
+ leaves := &f.Leaves
+ leaves.SetBit(start, 1)
+ nReached = 1
+ for n := start; n != end; n = ordering[o] {
+ if p[n].Len > 0 && len(a[n]) > 0 {
+ nDist := dist[n]
+ candLen := p[n].Len + 1 // len for any candidate arc followed from n
+ for _, to := range a[n] {
+ leaves.SetBit(to.To, 1)
+ candDist := nDist + w(to.Label)
+ switch {
+ case p[to.To].Len == 0: // first path to node to.To
+ nReached++
+ case fBetter(candDist, dist[to.To]): // better distance
+ case candDist == dist[to.To] && iBetter(candLen, p[to.To].Len): // same distance but better path length
+ default:
+ continue
+ }
+ dist[to.To] = candDist
+ p[to.To] = PathEnd{From: n, Len: candLen}
+ if candLen > f.MaxLen {
+ f.MaxLen = candLen
+ }
+ }
+ leaves.SetBit(n, 0)
+ }
+ o++
+ if o == len(ordering) {
+ break
+ }
+ }
+ return
+}
+
+// Dijkstra finds shortest paths by Dijkstra's algorithm.
+//
+// Shortest means shortest distance where distance is the
+// sum of arc weights. Where multiple paths exist with the same distance,
+// a path with the minimum number of nodes is returned.
+//
+// As usual for Dijkstra's algorithm, arc weights must be non-negative.
+// Graphs may be directed or undirected. Loops and parallel arcs are
+// allowed.
+func (g LabeledAdjacencyList) Dijkstra(start, end NI, w WeightFunc) (f FromList, dist []float64, reached int) {
+ r := make([]tentResult, len(g))
+ for i := range r {
+ r[i].nx = NI(i)
+ }
+ f = NewFromList(len(g))
+ dist = make([]float64, len(g))
+ current := start
+ rp := f.Paths
+ rp[current] = PathEnd{Len: 1, From: -1} // path length at start is 1 node
+ cr := &r[current]
+ cr.dist = 0 // distance at start is 0.
+ cr.done = true // mark start done. it skips the heap.
+ nDone := 1 // accumulated for a return value
+ var t tent
+ for current != end {
+ nextLen := rp[current].Len + 1
+ for _, nb := range g[current] {
+ // d.arcVis++
+ hr := &r[nb.To]
+ if hr.done {
+ continue // skip nodes already done
+ }
+ dist := cr.dist + w(nb.Label)
+ vl := rp[nb.To].Len
+ visited := vl > 0
+ if visited {
+ if dist > hr.dist {
+ continue // distance is worse
+ }
+ // tie breaker is a nice touch and doesn't seem to
+ // impact performance much.
+ if dist == hr.dist && nextLen >= vl {
+ continue // distance same, but number of nodes is no better
+ }
+ }
+ // the path through current to this node is shortest so far.
+ // record new path data for this node and update tentative set.
+ hr.dist = dist
+ rp[nb.To].Len = nextLen
+ rp[nb.To].From = current
+ if visited {
+ heap.Fix(&t, hr.fx)
+ } else {
+ heap.Push(&t, hr)
+ }
+ }
+ //d.ndVis++
+ if len(t) == 0 {
+ return f, dist, nDone // no more reachable nodes. AllPaths normal return
+ }
+ // new current is node with smallest tentative distance
+ cr = heap.Pop(&t).(*tentResult)
+ cr.done = true
+ nDone++
+ current = cr.nx
+ dist[current] = cr.dist // store final distance
+ }
+ // normal return for single shortest path search
+ return f, dist, -1
+}
+
+// DijkstraPath finds a single shortest path.
+//
+// Returned is the path and distance as returned by FromList.PathTo.
+func (g LabeledAdjacencyList) DijkstraPath(start, end NI, w WeightFunc) ([]NI, float64) {
+ f, dist, _ := g.Dijkstra(start, end, w)
+ return f.PathTo(end, nil), dist[end]
+}
+
+// tent implements container/heap
+func (t tent) Len() int { return len(t) }
+func (t tent) Less(i, j int) bool { return t[i].dist < t[j].dist }
+func (t tent) Swap(i, j int) {
+ t[i], t[j] = t[j], t[i]
+ t[i].fx = i
+ t[j].fx = j
+}
+func (s *tent) Push(x interface{}) {
+ nd := x.(*tentResult)
+ nd.fx = len(*s)
+ *s = append(*s, nd)
+}
+func (s *tent) Pop() interface{} {
+ t := *s
+ last := len(t) - 1
+ *s = t[:last]
+ return t[last]
+}
+
+type tentResult struct {
+ dist float64 // tentative distance, sum of arc weights
+ nx NI // slice index, "node id"
+ fx int // heap.Fix index
+ done bool
+}
+
+type tent []*tentResult