1// Copyright 2014 Sonia Keys 2// License MIT: http://opensource.org/licenses/MIT 3 4package graph 5 6import ( 7 "errors" 8 "fmt" 9 10 "github.com/soniakeys/bits" 11) 12 13// dir_RO.go is code generated from dir_cg.go by directives in graph.go. 14// Editing dir_cg.go is okay. It is the code generation source. 15// DO NOT EDIT dir_RO.go. 16// The RO means read only and it is upper case RO to slow you down a bit 17// in case you start to edit the file. 18 19// Balanced returns true if for every node in g, in-degree equals out-degree. 20// 21// There are equivalent labeled and unlabeled versions of this method. 22func (g Directed) Balanced() bool { 23 for n, in := range g.InDegree() { 24 if in != len(g.AdjacencyList[n]) { 25 return false 26 } 27 } 28 return true 29} 30 31// Copy makes a deep copy of g. 32// Copy also computes the arc size ma, the number of arcs. 33// 34// There are equivalent labeled and unlabeled versions of this method. 35func (g Directed) Copy() (c Directed, ma int) { 36 l, s := g.AdjacencyList.Copy() 37 return Directed{l}, s 38} 39 40// Cyclic determines if g contains a cycle, a non-empty path from a node 41// back to itself. 42// 43// Cyclic returns true if g contains at least one cycle. It also returns 44// an example of an arc involved in a cycle. 45// Cyclic returns false if g is acyclic. 46// 47// Also see Topological, which detects cycles. 48// 49// There are equivalent labeled and unlabeled versions of this method. 50func (g Directed) Cyclic() (cyclic bool, fr NI, to NI) { 51 a := g.AdjacencyList 52 fr, to = -1, -1 53 temp := bits.New(len(a)) 54 perm := bits.New(len(a)) 55 var df func(int) 56 df = func(n int) { 57 switch { 58 case temp.Bit(n) == 1: 59 cyclic = true 60 return 61 case perm.Bit(n) == 1: 62 return 63 } 64 temp.SetBit(n, 1) 65 for _, nb := range a[n] { 66 df(int(nb)) 67 if cyclic { 68 if fr < 0 { 69 fr, to = NI(n), nb 70 } 71 return 72 } 73 } 74 temp.SetBit(n, 0) 75 perm.SetBit(n, 1) 76 } 77 for n := range a { 78 if perm.Bit(n) == 1 { 79 continue 80 } 81 if df(n); cyclic { // short circuit as soon as a cycle is found 82 break 83 } 84 } 85 return 86} 87 88// DegreeCentralization returns out-degree centralization. 89// 90// Out-degree of a node is one measure of node centrality and is directly 91// available from the adjacency list representation. This allows degree 92// centralization for the graph to be very efficiently computed. 93// 94// The value returned is from 0 to 1 inclusive for simple directed graphs of 95// two or more nodes. As a special case, 0 is returned for graphs of 0 or 1 96// nodes. The value returned can be > 1 for graphs with loops or parallel 97// edges. 98// 99// In-degree centralization can be computed as DegreeCentralization of the 100// transpose. 101// 102// There are equivalent labeled and unlabeled versions of this method. 103func (g Directed) DegreeCentralization() float64 { 104 a := g.AdjacencyList 105 if len(a) <= 1 { 106 return 0 107 } 108 var max, sum int 109 for _, to := range a { 110 if len(to) > max { 111 max = len(to) 112 } 113 sum += len(to) 114 } 115 l1 := len(a) - 1 116 return float64(len(a)*max-sum) / float64(l1*l1) 117} 118 119// Dominators computes the immediate dominator for each node reachable from 120// start. 121// 122// The slice returned as Dominators.Immediate will have the length of 123// g.AdjacencyList. Nodes without a path to end will have a value of -1. 124// 125// See also the method Doms. Internally Dominators must construct the 126// transpose of g and also compute a postordering of a spanning tree of the 127// subgraph reachable from start. If you happen to have either of these 128// computed anyway, it can be more efficient to call Doms directly. 129func (g Directed) Dominators(start NI) Dominators { 130 a := g.AdjacencyList 131 l := len(a) 132 // ExampleDoms shows traditional depth-first postorder, but it works to 133 // generate a reverse preorder. Also breadth-first works instead of 134 // depth-first and may allow Doms to run a little faster by presenting 135 // a shallower tree. 136 post := make([]NI, l) 137 a.BreadthFirst(start, func(n NI) { 138 l-- 139 post[l] = n 140 }) 141 tr, _ := g.Transpose() 142 return g.Doms(tr, post[l:]) 143} 144 145// Doms computes either immediate dominators or postdominators. 146// 147// The slice returned as Dominators.Immediate will have the length of 148// g.AdjacencyList. Nodes without a path to end will have a value of -1. 149// 150// But see also the simpler methods Dominators and PostDominators. 151// 152// Doms requires argument tr to be the transpose graph of receiver g, 153// and requres argument post to be a post ordering of receiver g. More 154// specifically a post ordering of a spanning tree of the subgraph reachable 155// from some start node in g. The start node will always be the last node in 156// this postordering so it does not need to passed as a separate argument. 157// 158// Doms can be used to construct either dominators or postdominators. 159// To construct dominators on a graph f, generate a postordering p on f 160// and call f.Doms(f.Transpose(), p). To construct postdominators, generate 161// the transpose t first, then a postordering p on t (not f), and call 162// t.Doms(f, p). 163// 164// Caution: The argument tr is retained in the returned Dominators object 165// and is used by the method Dominators.Frontier. It is not deep-copied 166// so it is invalid to call Doms, modify the tr graph, and then call Frontier. 167func (g Directed) Doms(tr Directed, post []NI) Dominators { 168 a := g.AdjacencyList 169 dom := make([]NI, len(a)) 170 pi := make([]int, len(a)) 171 for i, n := range post { 172 pi[n] = i 173 } 174 intersect := func(b1, b2 NI) NI { 175 for b1 != b2 { 176 for pi[b1] < pi[b2] { 177 b1 = dom[b1] 178 } 179 for pi[b2] < pi[b1] { 180 b2 = dom[b2] 181 } 182 } 183 return b1 184 } 185 for n := range dom { 186 dom[n] = -1 187 } 188 start := post[len(post)-1] 189 dom[start] = start 190 for changed := false; ; changed = false { 191 for i := len(post) - 2; i >= 0; i-- { 192 b := post[i] 193 var im NI 194 fr := tr.AdjacencyList[b] 195 var j int 196 var fp NI 197 for j, fp = range fr { 198 if dom[fp] >= 0 { 199 im = fp 200 break 201 } 202 } 203 for _, p := range fr[j:] { 204 if dom[p] >= 0 { 205 im = intersect(im, p) 206 } 207 } 208 if dom[b] != im { 209 dom[b] = im 210 changed = true 211 } 212 } 213 if !changed { 214 return Dominators{dom, tr} 215 } 216 } 217} 218 219// PostDominators computes the immediate postdominator for each node that can 220// reach node end. 221// 222// The slice returned as Dominators.Immediate will have the length of 223// g.AdjacencyList. Nodes without a path to end will have a value of -1. 224// 225// See also the method Doms. Internally Dominators must construct the 226// transpose of g and also compute a postordering of a spanning tree of the 227// subgraph of the transpose reachable from end. If you happen to have either 228// of these computed anyway, it can be more efficient to call Doms directly. 229// 230// See the method Doms anyway for the caution note. PostDominators calls 231// Doms internally, passing receiver g as Doms argument tr. The caution means 232// that it is invalid to call PostDominators, modify the graph g, then call 233// Frontier. 234func (g Directed) PostDominators(end NI) Dominators { 235 tr, _ := g.Transpose() 236 a := tr.AdjacencyList 237 l := len(a) 238 post := make([]NI, l) 239 a.BreadthFirst(end, func(n NI) { 240 l-- 241 post[l] = n 242 }) 243 return tr.Doms(g, post[l:]) 244} 245 246// called from Dominators.Frontier via interface 247func (from Directed) domFrontiers(d Dominators) DominanceFrontiers { 248 im := d.Immediate 249 f := make(DominanceFrontiers, len(im)) 250 for i := range f { 251 if im[i] >= 0 { 252 f[i] = map[NI]struct{}{} 253 } 254 } 255 for b, fr := range from.AdjacencyList { 256 if len(fr) < 2 { 257 continue 258 } 259 imb := im[b] 260 for _, p := range fr { 261 for runner := p; runner != imb; runner = im[runner] { 262 f[runner][NI(b)] = struct{}{} 263 } 264 } 265 } 266 return f 267} 268 269// Eulerian scans a directed graph to determine if it is Eulerian. 270// 271// If the graph represents an Eulerian cycle, it returns -1, -1, nil. 272// 273// If the graph does not represent an Eulerian cycle but does represent an 274// Eulerian path, it returns the start and end nodes of the path, and nil. 275// 276// Otherwise it returns an error indicating a reason the graph is non-Eulerian. 277// Also in this case it returns a relevant node in either start or end. 278// 279// See also method EulerianStart, which short-circuits when it finds a start 280// node whereas this method completely validates a graph as Eulerian. 281// 282// There are equivalent labeled and unlabeled versions of this method. 283func (g Directed) Eulerian() (start, end NI, err error) { 284 ind := g.InDegree() 285 start = -1 286 end = -1 287 for n, to := range g.AdjacencyList { 288 switch { 289 case len(to) > ind[n]: 290 if start >= 0 { 291 return NI(n), -1, errors.New("multiple start candidates") 292 } 293 if len(to) > ind[n]+1 { 294 return NI(n), -1, errors.New("excessive out-degree") 295 } 296 start = NI(n) 297 case ind[n] > len(to): 298 if end >= 0 { 299 return -1, NI(n), errors.New("multiple end candidates") 300 } 301 if ind[n] > len(to)+1 { 302 return -1, NI(n), errors.New("excessive in-degree") 303 } 304 end = NI(n) 305 } 306 } 307 return start, end, nil 308} 309 310// EulerianCycle finds an Eulerian cycle in a directed multigraph. 311// 312// * If g has no nodes, result is nil, nil. 313// 314// * If g is Eulerian, result is an Eulerian cycle with err = nil. 315// The first element of the result represents only a start node. 316// The remaining elements represent the half arcs of the cycle. 317// 318// * Otherwise, result is nil, with a non-nil error giving a reason the graph 319// is not Eulerian. 320// 321// Internally, EulerianCycle copies the entire graph g. 322// See EulerianCycleD for a more space efficient version. 323// 324// There are nearly equivalent labeled and unlabeled versions of this method. 325// In the labeled version the first element of of the 326func (g Directed) EulerianCycle() ([]NI, error) { 327 c, m := g.Copy() 328 return c.EulerianCycleD(m) 329} 330 331// EulerianCycleD finds an Eulerian cycle in a directed multigraph. 332// 333// EulerianCycleD is destructive on its receiver g. See EulerianCycle for 334// a non-destructive version. 335// 336// Argument ma must be the correct arc size, or number of arcs in g. 337// 338// * If g has no nodes, result is nil, nil. 339// 340// * If g is Eulerian, result is an Eulerian cycle with err = nil. 341// The first element of the result represents only a start node. 342// The remaining elements represent the half arcs of the cycle. 343// 344// * Otherwise, result is nil, with a non-nil error giving a reason the graph 345// is not Eulerian. 346// 347// There are equivalent labeled and unlabeled versions of this method. 348func (g Directed) EulerianCycleD(ma int) ([]NI, error) { 349 // algorithm adapted from "Sketch of Eulerian Circuit Algorithm" by 350 // Carl Lee, accessed at http://www.ms.uky.edu/~lee/ma515fa10/euler.pdf. 351 if g.Order() == 0 { 352 return nil, nil 353 } 354 e := newEulerian(g.AdjacencyList, ma) 355 e.p[0] = 0 356 for e.s >= 0 { 357 v := e.top() // v is node that starts cycle 358 e.push() 359 // if Eulerian, we'll always come back to starting node 360 if e.top() != v { 361 return nil, errors.New("not Eulerian") 362 } 363 e.keep() 364 } 365 if !e.uv.AllZeros() { 366 return nil, errors.New("not strongly connected") 367 } 368 return e.p, nil 369} 370 371// EulerianPath finds an Eulerian path in a directed multigraph. 372// 373// * If g has no nodes, result is nil, nil. 374// 375// * If g has an Eulerian path, result is an Eulerian path with err = nil. 376// The first element of the result represents only a start node. 377// The remaining elements represent the half arcs of the path. 378// 379// * Otherwise, result is nil, with a non-nil error giving a reason the graph 380// is not Eulerian. 381// 382// Internally, EulerianPath copies the entire graph g. 383// See EulerianPathD for a more space efficient version. 384// 385// There are equivalent labeled and unlabeled versions of this method. 386func (g Directed) EulerianPath() ([]NI, error) { 387 c, m := g.Copy() 388 start, err := c.EulerianStart() 389 if err != nil { 390 return nil, err 391 } 392 if start < 0 { 393 start = 0 394 } 395 return c.EulerianPathD(m, start) 396} 397 398// EulerianPathD finds an Eulerian path in a directed multigraph. 399// 400// EulerianPathD is destructive on its receiver g. See EulerianPath for 401// a non-destructive version. 402// 403// Argument ma must be the correct arc size, or number of arcs in g. 404// Argument start must be a valid start node for the path. 405// 406// * If g has no nodes, result is nil, nil. 407// 408// * If g has an Eulerian path starting at start, result is an Eulerian path 409// with err = nil. 410// The first element of the result represents only a start node. 411// The remaining elements represent the half arcs of the path. 412// 413// * Otherwise, result is nil, with a non-nil error giving a reason the graph 414// is not Eulerian. 415// 416// There are equivalent labeled and unlabeled versions of this method. 417func (g Directed) EulerianPathD(ma int, start NI) ([]NI, error) { 418 if g.Order() == 0 { 419 return nil, nil 420 } 421 e := newEulerian(g.AdjacencyList, ma) 422 e.p[0] = start 423 // unlike EulerianCycle, the first path doesn't have to be a cycle. 424 e.push() 425 e.keep() 426 for e.s >= 0 { 427 start = e.top() 428 e.push() 429 // paths after the first must be cycles though 430 // (as long as there are nodes on the stack) 431 if e.top() != start { 432 return nil, errors.New("no Eulerian path") 433 } 434 e.keep() 435 } 436 if !e.uv.AllZeros() { 437 return nil, errors.New("no Eulerian path") 438 } 439 return e.p, nil 440} 441 442// EulerianStart finds a candidate start node for an Eulerian path. 443// 444// A candidate start node in the directed case has out-degree one greater then 445// in-degree. EulerianStart scans the graph returning immediately with the 446// node (and err == nil) when it finds such a candidate. 447// 448// EulerianStart also returns immediately with an error if it finds the graph 449// cannot contain an Eulerian path. In this case it also returns a relevant 450// node. 451// 452// If the scan completes without finding a candidate start node, the graph 453// represents an Eulerian cycle. In this case it returns -1, nil, and any 454// node can be chosen as a start node for an eulerian path. 455// 456// See also method Eulerian, which completely validates a graph as Eulerian 457// whereas this method short-curcuits when it finds a start node. 458// 459// There are equivalent labeled and unlabeled versions of this method. 460func (g Directed) EulerianStart() (start NI, err error) { 461 ind := g.InDegree() 462 end := -1 463 for n, to := range g.AdjacencyList { 464 switch { 465 case len(to) > ind[n]: 466 if len(to) == ind[n]+1 { 467 return NI(n), nil // candidate start 468 } 469 return -1, errors.New("excessive out-degree") 470 case ind[n] > len(to): 471 if end >= 0 { 472 return NI(n), errors.New("multiple end candidates") 473 } 474 if ind[n] > len(to)+1 { 475 return NI(n), errors.New("excessive in-degree") 476 } 477 end = n 478 } 479 } 480 return -1, nil // cycle 481} 482 483type eulerian struct { 484 g AdjacencyList // working copy of graph, it gets consumed 485 m int // number of arcs in g, updated as g is consumed 486 uv bits.Bits // unvisited 487 // low end of p is stack of unfinished nodes 488 // high end is finished path 489 p []NI // stack + path 490 s int // stack pointer 491} 492 493func newEulerian(g AdjacencyList, m int) *eulerian { 494 e := &eulerian{ 495 g: g, 496 m: m, 497 uv: bits.New(len(g)), 498 p: make([]NI, m+1), 499 } 500 e.uv.SetAll() 501 return e 502} 503 504// starting with the node on top of the stack, move nodes with no arcs. 505func (e *eulerian) keep() { 506 for e.s >= 0 { 507 n := e.top() 508 if len(e.g[n]) > 0 { 509 break 510 } 511 e.p[e.m] = n 512 e.s-- 513 e.m-- 514 } 515} 516 517func (e *eulerian) top() NI { 518 return e.p[e.s] 519} 520 521// MaximalNonBranchingPaths finds all paths in a directed graph that are 522// "maximal" and "non-branching". 523// 524// A non-branching path is one where path nodes other than the first and last 525// have exactly one arc leading to the node and one arc leading from the node, 526// thus there is no possibility to branch away to a different path. 527// 528// A maximal non-branching path cannot be extended to a longer non-branching 529// path by including another node at either end. 530// 531// In the case of a cyclic non-branching path, the first and last nodes 532// of the path will be the same node, indicating an isolated cycle. 533// 534// The method calls the emit argument for each path or isolated cycle in g, 535// as long as emit returns true. If emit returns false, 536// MaximalNonBranchingPaths returns immediately. 537// 538// There are equivalent labeled and unlabeled versions of this method. 539func (g Directed) MaximalNonBranchingPaths(emit func([]NI) bool) { 540 a := g.AdjacencyList 541 ind := g.InDegree() 542 uv := bits.New(g.Order()) 543 uv.SetAll() 544 for v, vTo := range a { 545 if !(ind[v] == 1 && len(vTo) == 1) { 546 for _, w := range vTo { 547 n := []NI{NI(v), w} 548 uv.SetBit(v, 0) 549 uv.SetBit(int(w), 0) 550 wTo := a[w] 551 for ind[w] == 1 && len(wTo) == 1 { 552 u := wTo[0] 553 n = append(n, u) 554 uv.SetBit(int(u), 0) 555 w = u 556 wTo = a[w] 557 } 558 if !emit(n) { // n is a path 559 return 560 } 561 } 562 } 563 } 564 // use uv.From rather than uv.Iterate. 565 // Iterate doesn't work here because we're modifying uv 566 for b := uv.OneFrom(0); b >= 0; b = uv.OneFrom(b + 1) { 567 v := NI(b) 568 n := []NI{v} 569 for w := v; ; { 570 w = a[w][0] 571 uv.SetBit(int(w), 0) 572 n = append(n, w) 573 if w == v { 574 break 575 } 576 } 577 if !emit(n) { // n is an isolated cycle 578 return 579 } 580 } 581} 582 583// InDegree computes the in-degree of each node in g 584// 585// There are equivalent labeled and unlabeled versions of this method. 586func (g Directed) InDegree() []int { 587 ind := make([]int, g.Order()) 588 for _, nbs := range g.AdjacencyList { 589 for _, nb := range nbs { 590 ind[nb]++ 591 } 592 } 593 return ind 594} 595 596// AddNode maps a node in a supergraph to a subgraph node. 597// 598// Argument p must be an NI in supergraph s.Super. AddNode panics if 599// p is not a valid node index of s.Super. 600// 601// AddNode is idempotent in that it does not add a new node to the subgraph if 602// a subgraph node already exists mapped to supergraph node p. 603// 604// The mapped subgraph NI is returned. 605func (s *DirectedSubgraph) AddNode(p NI) (b NI) { 606 if int(p) < 0 || int(p) >= s.Super.Order() { 607 panic(fmt.Sprint("AddNode: NI ", p, " not in supergraph")) 608 } 609 if b, ok := s.SubNI[p]; ok { 610 return b 611 } 612 a := s.Directed.AdjacencyList 613 b = NI(len(a)) 614 s.Directed.AdjacencyList = append(a, nil) 615 s.SuperNI = append(s.SuperNI, p) 616 s.SubNI[p] = b 617 return 618} 619 620// AddArc adds an arc to a subgraph. 621// 622// Arguments fr, to must be NIs in supergraph s.Super. As with AddNode, 623// AddArc panics if fr and to are not valid node indexes of s.Super. 624// 625// The arc specfied by fr, to must exist in s.Super. Further, the number of 626// parallel arcs in the subgraph cannot exceed the number of corresponding 627// parallel arcs in the supergraph. That is, each arc already added to the 628// subgraph counts against the arcs available in the supergraph. If a matching 629// arc is not available, AddArc returns an error. 630// 631// If a matching arc is available, subgraph nodes are added as needed, the 632// subgraph arc is added, and the method returns nil. 633func (s *DirectedSubgraph) AddArc(fr NI, to NI) error { 634 // verify supergraph NIs first, but without adding subgraph nodes just yet. 635 if int(fr) < 0 || int(fr) >= s.Super.Order() { 636 panic(fmt.Sprint("AddArc: NI ", fr, " not in supergraph")) 637 } 638 if int(to) < 0 || int(to) >= s.Super.Order() { 639 panic(fmt.Sprint("AddArc: NI ", to, " not in supergraph")) 640 } 641 // count existing matching arcs in subgraph 642 n := 0 643 a := s.Directed.AdjacencyList 644 if bf, ok := s.SubNI[fr]; ok { 645 if bt, ok := s.SubNI[to]; ok { 646 // both NIs already exist in subgraph, need to count arcs 647 bTo := to 648 bTo = bt 649 for _, t := range a[bf] { 650 if t == bTo { 651 n++ 652 } 653 } 654 } 655 } 656 // verify matching arcs are available in supergraph 657 for _, t := range (*s.Super).AdjacencyList[fr] { 658 if t == to { 659 if n > 0 { 660 n-- // match existing arc 661 continue 662 } 663 // no more existing arcs need to be matched. nodes can finally 664 // be added as needed and then the arc can be added. 665 bf := s.AddNode(fr) 666 to = s.AddNode(to) 667 s.Directed.AdjacencyList[bf] = 668 append(s.Directed.AdjacencyList[bf], to) 669 return nil // success 670 } 671 } 672 return errors.New("arc not available in supergraph") 673} 674 675// InduceList constructs a node-induced subgraph. 676// 677// The subgraph is induced on receiver graph g. Argument l must be a list of 678// NIs in receiver graph g. Receiver g becomes the supergraph of the induced 679// subgraph. 680// 681// Duplicate NIs are allowed in list l. The duplicates are effectively removed 682// and only a single corresponding node is created in the subgraph. Subgraph 683// NIs are mapped in the order of list l, execpt for ignoring duplicates. 684// NIs in l that are not in g will panic. 685// 686// Returned is the constructed Subgraph object containing the induced subgraph 687// and the mappings to the supergraph. 688func (g *Directed) InduceList(l []NI) *DirectedSubgraph { 689 sub, sup := mapList(l) 690 return &DirectedSubgraph{ 691 Super: g, 692 SubNI: sub, 693 SuperNI: sup, 694 Directed: Directed{ 695 g.AdjacencyList.induceArcs(sub, sup), 696 }} 697} 698 699// InduceBits constructs a node-induced subgraph. 700// 701// The subgraph is induced on receiver graph g. Argument t must be a bitmap 702// representing NIs in receiver graph g. Receiver g becomes the supergraph 703// of the induced subgraph. NIs in t that are not in g will panic. 704// 705// Returned is the constructed Subgraph object containing the induced subgraph 706// and the mappings to the supergraph. 707func (g *Directed) InduceBits(t bits.Bits) *DirectedSubgraph { 708 sub, sup := mapBits(t) 709 return &DirectedSubgraph{ 710 Super: g, 711 SubNI: sub, 712 SuperNI: sup, 713 Directed: Directed{ 714 g.AdjacencyList.induceArcs(sub, sup), 715 }} 716} 717 718// IsTree identifies trees in directed graphs. 719// 720// Return value isTree is true if the subgraph reachable from root is a tree. 721// Further, return value allTree is true if the entire graph g is reachable 722// from root. 723// 724// There are equivalent labeled and unlabeled versions of this method. 725func (g Directed) IsTree(root NI) (isTree, allTree bool) { 726 a := g.AdjacencyList 727 v := bits.New(len(a)) 728 v.SetAll() 729 var df func(NI) bool 730 df = func(n NI) bool { 731 if v.Bit(int(n)) == 0 { 732 return false 733 } 734 v.SetBit(int(n), 0) 735 for _, to := range a[n] { 736 if !df(to) { 737 return false 738 } 739 } 740 return true 741 } 742 isTree = df(root) 743 return isTree, isTree && v.AllZeros() 744} 745 746// PageRank computes a significance score for each node of a graph. 747// 748// The algorithm is credited to Google founders Brin and Lawrence. 749// 750// Argument d is a damping factor. Reportedly a value of .85 works well. 751// Argument n is a number of iterations. Reportedly values of 20 to 50 752// work well. 753// 754// Returned is the PageRank score for each node of g. 755// 756// There are equivalent labeled and unlabeled versions of this method. 757func (g Directed) PageRank(d float64, n int) []float64 { 758 // Following "PageRank Explained" by Ian Rogers, accessed at 759 // http://www.cs.princeton.edu/~chazelle/courses/BIB/pagerank.htm 760 a := g.AdjacencyList 761 p0 := make([]float64, len(a)) 762 p1 := make([]float64, len(a)) 763 for i := range p0 { 764 p0[i] = 1 765 } 766 d1 := 1 - d 767 for ; n > 0; n-- { 768 for i := range p1 { 769 p1[i] = d1 770 } 771 for fr, to := range a { 772 f := d / float64(len(to)) 773 for _, to := range to { 774 p1[to] += p0[fr] * f 775 } 776 } 777 p0, p1 = p1, p0 778 } 779 return p0 780} 781 782// StronglyConnectedComponents identifies strongly connected components in 783// a directed graph. 784// 785// The method calls the emit function for each component identified. The 786// argument to emit is the node list of a component. The emit function must 787// return true for the method to continue identifying components. If emit 788// returns false, the method returns immediately. 789// 790// Note well: The backing slice for the node list passed to emit is reused 791// across emit calls. If you need to retain the node list you must copy it. 792// 793// The components emitted represent a partition of the nodes in g. 794// So for example, if the first component emitted has the same length as g 795// then it will be the only component and it means the entire graph g is 796// strongly connected. 797// 798// See also Condensation which returns a condensation graph in addition 799// to the strongly connected components. 800// 801// There are equivalent labeled and unlabeled versions of this method. 802// 803// The algorithm here is by David Pearce. See also alt.SCCPathBased and 804// alt.SCCTarjan. 805func (g Directed) StronglyConnectedComponents(emit func([]NI) bool) { 806 // See Algorithm 3 PEA FIND SCC2(V,E) in "An Improved Algorithm for 807 // Finding the Strongly Connected Components of a Directed Graph" 808 // by David J. Pearce. 809 a := g.AdjacencyList 810 rindex := make([]int, len(a)) 811 var S, scc []NI 812 index := 1 813 c := len(a) - 1 814 var visit func(NI) bool 815 visit = func(v NI) bool { 816 root := true 817 rindex[v] = index 818 index++ 819 for _, w := range a[v] { 820 if rindex[w] == 0 { 821 if !visit(w) { 822 return false 823 } 824 } 825 if rindex[w] < rindex[v] { 826 rindex[v] = rindex[w] 827 root = false 828 } 829 } 830 if !root { 831 S = append(S, v) 832 return true 833 } 834 scc = scc[:0] 835 index-- 836 for last := len(S) - 1; last >= 0; last-- { 837 w := S[last] 838 if rindex[v] > rindex[w] { 839 break 840 } 841 S = S[:last] 842 rindex[w] = c 843 scc = append(scc, w) 844 index-- 845 } 846 rindex[v] = c 847 c-- 848 return emit(append(scc, v)) 849 } 850 for v := range a { 851 if rindex[v] == 0 && !visit(NI(v)) { 852 break 853 } 854 } 855} 856 857// Condensation returns strongly connected components and their 858// condensation graph. 859// 860// A condensation represents a directed acyclic graph. 861// Components are ordered in a reverse topological ordering. 862// 863// See also StronglyConnectedComponents, which returns the components only. 864// 865// There are equivalent labeled and unlabeled versions of this method. 866func (g Directed) Condensation() (scc [][]NI, cd AdjacencyList) { 867 a := g.AdjacencyList 868 b := make([]NI, len(a)) // backing slice for scc 869 g.StronglyConnectedComponents(func(c []NI) bool { 870 n := copy(b, c) 871 scc = append(scc, b[:n]) 872 b = b[n:] 873 return true 874 }) 875 cd = make(AdjacencyList, len(scc)) // return value 876 cond := make([]NI, len(a)) // mapping from g node to cd node 877 for cn, c := range scc { 878 for _, n := range c { 879 cond[n] = NI(cn) // map g node to cd node 880 } 881 var tos []NI // list of 'to' nodes 882 m := bits.New(len(cd)) // tos map 883 m.SetBit(cn, 1) 884 for _, n := range c { 885 for _, to := range a[n] { 886 if ct := cond[to]; m.Bit(int(ct)) == 0 { 887 m.SetBit(int(ct), 1) 888 tos = append(tos, ct) 889 } 890 } 891 } 892 cd[cn] = tos 893 } 894 return 895} 896 897// Topological computes a topological ordering of a directed acyclic graph. 898// 899// For an acyclic graph, return value ordering is a permutation of node numbers 900// in topologically sorted order and cycle will be nil. If the graph is found 901// to be cyclic, ordering will be nil and cycle will be the path of a found 902// cycle. 903// 904// There are equivalent labeled and unlabeled versions of this method. 905func (g Directed) Topological() (ordering, cycle []NI) { 906 i := -1 907 return g.dfTopo(func() NI { 908 i++ 909 if i < g.Order() { 910 return NI(i) 911 } 912 return -1 913 }) 914} 915 916func (g Directed) dfTopo(f func() NI) (ordering, cycle []NI) { 917 a := g.AdjacencyList 918 ordering = make([]NI, len(a)) 919 i := len(ordering) 920 temp := bits.New(len(a)) 921 perm := bits.New(len(a)) 922 var cycleFound bool 923 var cycleStart NI 924 var df func(NI) 925 df = func(n NI) { 926 switch { 927 case temp.Bit(int(n)) == 1: 928 cycleFound = true 929 cycleStart = n 930 return 931 case perm.Bit(int(n)) == 1: 932 return 933 } 934 temp.SetBit(int(n), 1) 935 for _, nb := range a[n] { 936 df(nb) 937 if cycleFound { 938 if cycleStart >= 0 { 939 // a little hack: orderng won't be needed so repurpose the 940 // slice as cycle. this is read out in reverse order 941 // as the recursion unwinds. 942 x := len(ordering) - 1 - len(cycle) 943 ordering[x] = n 944 cycle = ordering[x:] 945 if n == cycleStart { 946 cycleStart = -1 947 } 948 } 949 return 950 } 951 } 952 temp.SetBit(int(n), 0) 953 perm.SetBit(int(n), 1) 954 i-- 955 ordering[i] = n 956 } 957 for { 958 n := f() 959 if n < 0 { 960 return ordering[i:], nil 961 } 962 if perm.Bit(int(n)) == 1 { 963 continue 964 } 965 df(n) 966 if cycleFound { 967 return nil, cycle 968 } 969 } 970} 971 972// TopologicalKahn computes a topological ordering of a directed acyclic graph. 973// 974// For an acyclic graph, return value ordering is a permutation of node numbers 975// in topologically sorted order and cycle will be nil. If the graph is found 976// to be cyclic, ordering will be nil and cycle will be the path of a found 977// cycle. 978// 979// This function is based on the algorithm by Arthur Kahn and requires the 980// transpose of g be passed as the argument. 981// 982// There are equivalent labeled and unlabeled versions of this method. 983func (g Directed) TopologicalKahn(tr Directed) (ordering, cycle []NI) { 984 // code follows Wikipedia pseudocode. 985 var L, S []NI 986 // rem for "remaining edges," this function makes a local copy of the 987 // in-degrees and consumes that instead of consuming an input. 988 rem := make([]int, g.Order()) 989 for n, fr := range tr.AdjacencyList { 990 if len(fr) == 0 { 991 // accumulate "set of all nodes with no incoming edges" 992 S = append(S, NI(n)) 993 } else { 994 // initialize rem from in-degree 995 rem[n] = len(fr) 996 } 997 } 998 for len(S) > 0 { 999 last := len(S) - 1 // "remove a node n from S" 1000 n := S[last] 1001 S = S[:last] 1002 L = append(L, n) // "add n to tail of L" 1003 for _, m := range g.AdjacencyList[n] { 1004 // WP pseudo code reads "for each node m..." but it means for each 1005 // node m *remaining in the graph.* We consume rem rather than 1006 // the graph, so "remaining in the graph" for us means rem[m] > 0. 1007 if rem[m] > 0 { 1008 rem[m]-- // "remove edge from the graph" 1009 if rem[m] == 0 { // if "m has no other incoming edges" 1010 S = append(S, m) // "insert m into S" 1011 } 1012 } 1013 } 1014 } 1015 // "If graph has edges," for us means a value in rem is > 0. 1016 for c, in := range rem { 1017 if in > 0 { 1018 // recover cyclic nodes 1019 for _, nb := range g.AdjacencyList[c] { 1020 if rem[nb] > 0 { 1021 cycle = append(cycle, NI(c)) 1022 break 1023 } 1024 } 1025 } 1026 } 1027 if len(cycle) > 0 { 1028 return nil, cycle 1029 } 1030 return L, nil 1031} 1032 1033// TopologicalSubgraph computes a topological ordering of a subgraph of a 1034// directed acyclic graph. 1035// 1036// The subgraph considered is that reachable from the specified node list. 1037// 1038// For an acyclic subgraph, return value ordering is a permutation of reachable 1039// node numbers in topologically sorted order and cycle will be nil. If the 1040// subgraph is found to be cyclic, ordering will be nil and cycle will be 1041// the path of a found cycle. 1042// 1043// There are equivalent labeled and unlabeled versions of this method. 1044func (g Directed) TopologicalSubgraph(nodes []NI) (ordering, cycle []NI) { 1045 i := -1 1046 return g.dfTopo(func() NI { 1047 i++ 1048 if i < len(nodes) { 1049 return nodes[i] 1050 } 1051 return -1 1052 }) 1053} 1054 1055// TransitiveClosure returns the transitive closure of directed graph g. 1056// 1057// The algorithm is Warren's, which works most naturally with an adjacency 1058// matrix representation. The returned transitive closure is left in this 1059// adjacency matrix representation. For a graph g of order n, matrix tc 1060// is returned as a length n slice of length n bits.Bits values, where 1061// tc[from].Bit(to) == 1 represents an arc of the transitive closure. 1062func (g Directed) TransitiveClosure() []bits.Bits { 1063 // construct adjacency matrix 1064 a := g.AdjacencyList 1065 t := make([]bits.Bits, len(a)) 1066 for n := range t { 1067 tn := bits.New(len(a)) 1068 for _, to := range a[n] { 1069 tn.SetBit(int(to), 1) 1070 } 1071 t[n] = tn 1072 } 1073 // above diagonal 1074 for i := 1; i < len(a); i++ { 1075 ti := t[i] 1076 for k := 0; k < i; k++ { 1077 if ti.Bit(k) == 1 { 1078 ti.Or(ti, t[k]) 1079 } 1080 } 1081 } 1082 // below diagonal 1083 for i, ti := range t[:len(a)-1] { 1084 for k := i + 1; k < len(a); k++ { 1085 if ti.Bit(k) == 1 { 1086 ti.Or(ti, t[k]) 1087 } 1088 } 1089 } 1090 return t 1091} 1092