1 /** 2 Some additional mathematical functions. 3 4 Copyright: © 2018 Arne Ludwig <arne.ludwig@posteo.de> 5 License: Subject to the terms of the MIT license, as written in the 6 included LICENSE file. 7 Authors: Arne Ludwig <arne.ludwig@posteo.de> 8 */ 9 module dentist.util.math; 10 11 import dentist.util.algorithm : cmpLexicographically, sliceBy; 12 import std.algorithm : 13 all, 14 copy, 15 countUntil, 16 cumulativeFold, 17 filter, 18 map, 19 max, 20 sort, 21 sum, 22 swap, 23 uniq; 24 import std.array : Appender, array; 25 import std.conv : to; 26 import std.exception : assertThrown; 27 import std.functional : binaryFun, unaryFun; 28 import std.range : 29 assumeSorted, 30 chain, 31 ElementType, 32 enumerate, 33 isForwardRange, 34 isRandomAccessRange, 35 retro, 36 save, 37 walkLength; 38 import std.traits : isCallable, isIntegral, isNumeric; 39 import std.typecons : Flag, No, Yes; 40 41 debug import std.stdio : writeln; 42 43 /// Calculate the mean of range. 44 ElementType!Range mean(Range)(Range values) if (isForwardRange!Range) 45 { 46 auto sum = values.save.sum; 47 auto length = values.walkLength.to!(ElementType!Range); 48 49 return sum / length; 50 } 51 52 unittest 53 { 54 { 55 auto values = [2, 4, 8]; 56 assert(values.mean == 4); 57 } 58 { 59 auto values = [1.0, 2.0, 3.0, 4.0]; 60 assert(values.mean == 2.5); 61 } 62 } 63 64 /// Calculate the median of range. 65 ElementType!Range median(Range)(Range values) if (__traits(compiles, sort(values))) 66 { 67 assert(values.length > 0, "median is undefined for empty set"); 68 auto middleIdx = values.length / 2; 69 auto useAverage = values.length > 1 && values.length % 2 == 0; 70 auto sortedValues = values.sort; 71 72 if (useAverage) 73 return (sortedValues[middleIdx - 1] + sortedValues[middleIdx]) / 2; 74 else 75 return values[middleIdx]; 76 } 77 78 unittest 79 { 80 { 81 auto values = [4, 2, 8]; 82 assert(values.median == 4); 83 } 84 { 85 auto values = [4, 3, 2, 8]; 86 assert(values.median == 3); 87 } 88 { 89 auto values = [4, 6, 2, 8]; 90 assert(values.median == 5); 91 } 92 { 93 auto values = [2, 1, 3, 0, 4, 9, 8, 5, 6, 3, 9]; 94 assert(values.median == 4); 95 } 96 { 97 auto values = [2.0, 1.0, 4.0, 3.0]; 98 assert(values.median == 2.5); 99 } 100 { 101 auto values = [2.0, 1.0, 4.0, 3.0, 5.0]; 102 assert(values.median == 3.0); 103 } 104 } 105 106 /// Calculate the Nxx (e.g. N50) of values. 107 ElementType!Range N(real xx, Range, Num)(Range values, Num totalSize) if (__traits(compiles, sort(values))) 108 { 109 static assert(0 < xx && xx < 100, "N" ~ xx.to!string ~ " is undefined for empty set"); 110 assert(values.length > 0, "N" ~ xx.to!string ~ " is undefined for empty set"); 111 auto xxPercentile = xx/100.0 * totalSize; 112 auto sortedValues = values.sort; 113 auto targetIndex = sortedValues 114 .retro 115 .cumulativeFold!"a + b"(cast(ElementType!Range) 0) 116 .countUntil!"a >= b"(xxPercentile); 117 118 if (targetIndex == values.length) 119 return 0; 120 else 121 return sortedValues[$ - targetIndex - 1]; 122 } 123 124 unittest 125 { 126 { 127 auto totalSize = 54; 128 auto values = [2, 3, 4, 5, 6, 7, 8, 9, 10]; 129 enum N50 = 8; 130 enum N10 = 10; 131 132 assert(N!50(values, totalSize) == N50); 133 assert(N!10(values, totalSize) == N10); 134 } 135 { 136 auto totalSize = 32; 137 auto values = [2, 2, 2, 3, 3, 4, 8, 8]; 138 enum N50 = 8; 139 140 assert(N!50(values, totalSize) == N50); 141 } 142 } 143 144 enum RoundingMode : byte 145 { 146 floor, 147 round, 148 ceil, 149 } 150 151 /** 152 Round x upward according to base, ie. returns the next integer larger or 153 equal to x which is divisible by base. 154 155 Returns: x rounded upward according to base. 156 */ 157 Integer ceil(Integer)(in Integer x, in Integer base) pure nothrow 158 if (isIntegral!Integer) 159 { 160 return x % base == 0 161 ? x 162 : (x / base + 1) * base; 163 } 164 165 /// 166 unittest 167 { 168 assert(ceil(8, 10) == 10); 169 assert(ceil(32, 16) == 32); 170 assert(ceil(101, 100) == 200); 171 } 172 173 /** 174 Round x downward according to base, ie. returns the next integer smaller or 175 equal to x which is divisible by base. 176 177 Returns: x rounded downward according to base. 178 */ 179 Integer floor(Integer)(in Integer x, in Integer base) pure nothrow 180 if (isIntegral!Integer) 181 { 182 return (x / base) * base; 183 } 184 185 /// 186 unittest 187 { 188 assert(floor(8, 10) == 0); 189 assert(floor(32, 16) == 32); 190 assert(floor(101, 100) == 100); 191 } 192 193 /// Returns the absolute difference between two numbers. 194 Num absdiff(Num)(in Num a, in Num b) pure nothrow if (isNumeric!Num) 195 { 196 return a > b 197 ? a - b 198 : b - a; 199 } 200 201 /// 202 unittest 203 { 204 assert(absdiff(2UL, 3UL) == 1UL); 205 assert(absdiff(-42, 13) == 55); 206 assert(absdiff(2.5, 5) == 2.5); 207 } 208 209 /// Returns the result of `ceil(a / b)` but uses integer arithmetic only. 210 Integer ceildiv(Integer)(in Integer a, in Integer b) pure nothrow if (isIntegral!Integer) 211 { 212 Integer resultSign = (a < 0) ^ (b < 0) ? -1 : 1; 213 214 return resultSign < 0 || a % b == 0 215 ? a / b 216 : a / b + resultSign; 217 } 218 219 /// 220 unittest 221 { 222 assert(ceildiv(0, 3) == 0); 223 assert(ceildiv(1UL, 3UL) == 1UL); 224 assert(ceildiv(2L, 3L) == 1L); 225 assert(ceildiv(3U, 3U) == 1U); 226 assert(ceildiv(4, 3) == 2); 227 assert(ceildiv(-4, 4) == -1); 228 assert(ceildiv(-4, 3) == -1); 229 } 230 231 class EdgeExistsException : Exception 232 { 233 pure nothrow @nogc @safe this( 234 string file = __FILE__, 235 size_t line = __LINE__, 236 Throwable nextInChain = null, 237 ) 238 { 239 super("edge cannot be inserted: edge already exists", file, line, nextInChain); 240 } 241 } 242 243 class MissingEdgeException : Exception 244 { 245 pure nothrow @nogc @safe this( 246 string file = __FILE__, 247 size_t line = __LINE__, 248 Throwable nextInChain = null, 249 ) 250 { 251 super("edge not found", file, line, nextInChain); 252 } 253 } 254 255 class MissingNodeException : Exception 256 { 257 pure nothrow @nogc @safe this( 258 string file = __FILE__, 259 size_t line = __LINE__, 260 Throwable nextInChain = null, 261 ) 262 { 263 super("node not found", file, line, nextInChain); 264 } 265 } 266 267 /// This structure represents a graph with optional edge 268 /// payloads. The graph is represented as a list of edges which is 269 /// particularly suited for sparse graphs. While the set of nodes is fixed the 270 /// set of edges is mutable. 271 struct Graph(Node, Weight = void, Flag!"isDirected" isDirected = No.isDirected, EdgePayload = void) 272 { 273 static enum isWeighted = !is(Weight == void); 274 static enum hasEdgePayload = !is(EdgePayload == void); 275 276 static struct Edge 277 { 278 protected Node _start; 279 protected Node _end; 280 281 static if (isWeighted) 282 Weight weight; 283 284 static if (hasEdgePayload) 285 EdgePayload payload; 286 287 /// Construct an edge. 288 this(Node start, Node end) 289 { 290 this._start = start; 291 this._end = end; 292 293 static if (!isDirected) 294 { 295 if (end < start) 296 { 297 swap(this._start, this._end); 298 } 299 } 300 } 301 302 static if (isWeighted) 303 { 304 /// ditto 305 this(Node start, Node end, Weight weight) 306 { 307 this(start, end); 308 this.weight = weight; 309 } 310 } 311 312 static if (hasEdgePayload && !is(EdgePayload : Weight)) 313 { 314 /// ditto 315 this(Node start, Node end, EdgePayload payload) 316 { 317 this(start, end); 318 this.payload = payload; 319 } 320 } 321 322 static if (isWeighted && hasEdgePayload) 323 { 324 /// ditto 325 this(Node start, Node end, Weight weight, EdgePayload payload) 326 { 327 this(start, end); 328 this.weight = weight; 329 this.payload = payload; 330 } 331 } 332 333 /// Get the start of this edge. For undirected graphs this is the 334 /// smaller of both incident nodes. 335 @property Node start() const pure nothrow 336 { 337 return _start; 338 } 339 340 /// Get the end of this edge. For undirected graphs this is the 341 /// larger of both incident nodes. 342 @property Node end() const pure nothrow 343 { 344 return _end; 345 } 346 347 /** 348 Get target of this edge beginning at node `from`. For undirected 349 graphs returns the other node of this edge. 350 351 Throws: MissingNodeException if this edge does not start in node `from`. 352 */ 353 Node target(Node from) const 354 { 355 static if (isDirected) 356 { 357 if (start == from) 358 { 359 return end; 360 } 361 else 362 { 363 throw new MissingNodeException(); 364 } 365 } 366 else 367 { 368 if (start == from) 369 { 370 return end; 371 } 372 else if (end == from) 373 { 374 return start; 375 } 376 else 377 { 378 throw new MissingNodeException(); 379 } 380 } 381 } 382 383 /** 384 Get source of this edge beginning at node `from`. For undirected 385 graphs returns the other node of this edge. 386 387 Throws: MissingNodeException if this edge does not end in node `from`. 388 */ 389 static if (isDirected) 390 { 391 Node source(Node from) const 392 { 393 if (end == from) 394 { 395 return start; 396 } 397 else 398 { 399 throw new MissingNodeException(); 400 } 401 } 402 } 403 else 404 { 405 alias source = target; 406 } 407 408 /// Two edges are equal iff their incident nodes (and weight) are the 409 /// same. 410 bool opEquals(in Edge other) const pure nothrow 411 { 412 static if (isWeighted) 413 { 414 return this.start == other.start && this.end == other.end 415 && this.weight == other.weight; 416 } 417 else 418 { 419 return this.start == other.start && this.end == other.end; 420 } 421 } 422 423 /// Orders edge lexicographically by `start`, `end`(, `weight`). 424 int opCmp(in Edge other) const pure nothrow 425 { 426 static if (isWeighted) 427 { 428 return cmpLexicographically!( 429 typeof(this), 430 "a.start", 431 "a.end", 432 "a.weight", 433 )(this, other); 434 } 435 else 436 { 437 return cmpLexicographically!( 438 typeof(this), 439 "a.start", 440 "a.end", 441 )(this, other); 442 } 443 } 444 445 private int compareNodes(in Edge other) const pure nothrow 446 { 447 return cmpLexicographically!( 448 typeof(this), 449 "a.start", 450 "a.end", 451 )(this, other); 452 } 453 454 /** 455 Returns the node that is connects this edge with other edge. In 456 case of undirected graphs this is just the common node of both 457 edges; in directed case this is the end node of this edge if it 458 matches the start node of other edge. 459 460 Throws: MissingNodeException if the connecting node is undefined. 461 */ 462 Node getConnectingNode(in Edge other) const 463 { 464 static if (isDirected) 465 { 466 if (this.end == other.start) 467 { 468 return this.end; 469 } 470 } 471 else 472 { 473 if (this.end == other.start || this.end == other.end) 474 { 475 return this.end; 476 } 477 else if (this.start == other.start || this.start == other.end) 478 { 479 return this.start; 480 } 481 } 482 483 throw new MissingNodeException(); 484 } 485 } 486 487 static bool orderByNodes(in Edge a, in Edge b) nothrow pure 488 { 489 return a.compareNodes(b) < 0; 490 } 491 492 static bool groupByNodes(in Edge a, in Edge b) nothrow pure 493 { 494 return a.compareNodes(b) == 0; 495 } 496 497 /// Construct an edge for this graph. 498 static Edge edge(T...)(T args) 499 { 500 return Edge(args); 501 } 502 503 protected Node[] _nodes; 504 protected Appender!(Edge[]) _edges; 505 506 /// The set (ordered list) of nodes. 507 @property const(Node[]) nodes() const nothrow pure 508 { 509 return _nodes; 510 } 511 512 private @property void nodes(Node[] nodes) 513 { 514 nodes.sort(); 515 516 this._nodes = nodes.uniq.array; 517 } 518 519 /// Get the set (ordered list) of edges in this graph. 520 @property auto edges() nothrow pure 521 { 522 return chain(_edges.data); 523 } 524 525 /// ditto 526 @property auto edges() const nothrow pure 527 { 528 return chain(_edges.data); 529 } 530 531 /** 532 Construct a graph from a set of nodes (and edges). Modifies `nodes` 533 while sorting but releases it after construction. 534 535 Throws: MissingNodeException if an edge has a node that is not present 536 in this graph . 537 Throws: EdgeExistsException if an edge already exists when trying 538 inserting it. 539 */ 540 this(Node[] nodes) 541 { 542 this.nodes = nodes; 543 } 544 545 /// ditto 546 this(Node[] nodes, Edge[] edges) 547 { 548 this(nodes); 549 550 _edges.reserve(edges.length); 551 foreach (edge; edges) 552 { 553 add(this, edge); 554 } 555 } 556 557 this(this) 558 { 559 _nodes = _nodes.dup; 560 } 561 562 /// Add a set of edges to this graph without any checks. 563 void bulkAddForce(R)(R edges) if (isForwardRange!R && is(ElementType!R == Edge)) 564 { 565 this._edges ~= edges; 566 _edges.data.sort; 567 } 568 569 /// Add an edge to this graph. 570 /// See_Also: `Edge add(Graph, Edge)` 571 void opOpAssign(string op)(Edge edge) if (op == "~") 572 { 573 add(this, edge); 574 } 575 576 /// Some pre-defined conflict handlers for `add`. 577 static struct ConflictStrategy 578 { 579 static if (isWeighted) 580 { 581 /// Return an edge with sum of both weights. If given payload will be 582 /// kept from existingEdge . 583 static Edge sumWeights(Edge existingEdge, Edge newEdge) 584 { 585 existingEdge.weight += newEdge.weight; 586 587 return existingEdge; 588 } 589 590 /// 591 unittest 592 { 593 auto g1 = Graph!(int, int)([1, 2]); 594 alias CS = g1.ConflictStrategy; 595 596 g1 ~= g1.edge(1, 2, 1); 597 598 auto addedEdge = g1.add!(CS.sumWeights)(g1.edge(1, 2, 1)); 599 600 assert(addedEdge.weight == 2); 601 } 602 } 603 604 /// Throw `EdgeExistsException`. 605 static inout(Edge) error(inout(Edge) existingEdge, inout(Edge) newEdge) 606 { 607 throw new EdgeExistsException(); 608 } 609 610 /// 611 unittest 612 { 613 auto g1 = Graph!int([1, 2]); 614 alias CS = g1.ConflictStrategy; 615 616 g1 ~= g1.edge(1, 2); 617 618 assertThrown!EdgeExistsException(g1.add!(CS.error)(g1.edge(1, 2))); 619 } 620 621 /// Replace the existingEdge by newEdge. 622 static inout(Edge) replace(inout(Edge) existingEdge, inout(Edge) newEdge) 623 { 624 return newEdge; 625 } 626 627 /// 628 unittest 629 { 630 auto g1 = Graph!(int, int)([1, 2]); 631 alias CS = g1.ConflictStrategy; 632 633 g1 ~= g1.edge(1, 2, 1); 634 635 auto addedEdge = g1.add!(CS.replace)(g1.edge(1, 2, 2)); 636 637 assert(addedEdge.weight == 2); 638 } 639 640 /// Keep existingEdge – discard newEdge. 641 static inout(Edge) keep(inout(Edge) existingEdge, inout(Edge) newEdge) 642 { 643 return existingEdge; 644 } 645 646 /// 647 unittest 648 { 649 auto g1 = Graph!(int, int)([1, 2]); 650 alias CS = g1.ConflictStrategy; 651 652 g1 ~= g1.edge(1, 2, 1); 653 654 auto addedEdge = g1.add!(CS.keep)(g1.edge(1, 2, 2)); 655 656 assert(addedEdge.weight == 1); 657 } 658 } 659 660 /// Forcibly add an edge to this graph. 661 protected Edge forceAdd(Edge edge) 662 { 663 _edges ~= edge; 664 _edges.data.sort; 665 666 return edge; 667 } 668 669 /// Replace an edge in this graph. 670 protected Edge replaceEdge(in size_t edgeIdx, Edge newEdge) 671 { 672 auto shouldSort = _edges.data[edgeIdx] != newEdge; 673 674 _edges.data[edgeIdx] = newEdge; 675 676 if (shouldSort) 677 { 678 _edges.data.sort; 679 } 680 681 return newEdge; 682 } 683 684 /// Check if edge/node exists in this graph. Ignores the weight if weighted. 685 bool opBinaryRight(string op)(Node node) const pure nothrow if (op == "in") 686 { 687 auto sortedNodes = assumeSorted(nodes); 688 689 return sortedNodes.contains(node); 690 } 691 692 /// ditto 693 bool has(Node node) const pure nothrow 694 { 695 return node in this; 696 } 697 698 /// Check if edge exists in this graph. Only the `start` and `end` node 699 /// will be compared. 700 bool opBinaryRight(string op)(Edge edge) const pure nothrow if (op == "in") 701 { 702 auto sortedEdges = assumeSorted!orderByNodes(edges); 703 704 return sortedEdges.contains(edge); 705 } 706 707 /// ditto 708 bool has(Edge edge) const pure nothrow 709 { 710 return edge in this; 711 } 712 713 /// Get the designated edge from this graph. Only the `start` and `end` 714 /// node will be compared. 715 auto get(Edge edge) 716 { 717 auto sortedEdges = assumeSorted!orderByNodes(edges); 718 auto existingEdges = sortedEdges.equalRange(edge); 719 720 if (existingEdges.empty) 721 { 722 throw new MissingEdgeException(); 723 } 724 else 725 { 726 return existingEdges.front; 727 } 728 } 729 730 /// 731 unittest 732 { 733 auto g1 = Graph!(int, int)([1, 2]); 734 735 auto e1 = g1.edge(1, 2, 1); 736 737 g1 ~= e1; 738 739 assert(g1.get(g1.edge(1, 2)) == e1); 740 assertThrown!MissingEdgeException(g1.get(g1.edge(1, 1))); 741 } 742 743 /// Returns the ndex of node `n` in the list of nodes. 744 size_t indexOf(Node n) const 745 { 746 auto sortedNodes = assumeSorted(nodes); 747 auto tristectedNodes = sortedNodes.trisect(n); 748 749 if (tristectedNodes[1].empty) 750 { 751 throw new MissingNodeException(); 752 } 753 754 return tristectedNodes[0].length; 755 } 756 757 /// 758 unittest 759 { 760 761 auto g1 = Graph!(int, int)([1, 2]); 762 763 assert(g1.indexOf(1) == 0); 764 assert(g1.indexOf(2) == 1); 765 assertThrown!MissingNodeException(g1.indexOf(3)); 766 } 767 768 static if (isDirected) 769 { 770 /// Returns a range of in/outgoing edges of node `n`. 771 auto inEdges(Node n) nothrow pure 772 { 773 return _edges.data[].filter!(e => e.end == n); 774 } 775 776 /// ditto 777 auto inEdges(Node n) const nothrow pure 778 { 779 return edges[].filter!(e => e.end == n); 780 } 781 782 /// ditto 783 auto outEdges(Node n) nothrow pure 784 { 785 return _edges.data[].filter!(e => e.start == n); 786 } 787 788 /// ditto 789 auto outEdges(Node n) const nothrow pure 790 { 791 return edges[].filter!(e => e.start == n); 792 } 793 794 /// 795 unittest 796 { 797 import std.algorithm : equal; 798 799 auto g1 = Graph!(int, void, Yes.isDirected)([1, 2, 3]); 800 801 g1 ~= g1.edge(1, 1); 802 g1 ~= g1.edge(1, 2); 803 g1 ~= g1.edge(2, 2); 804 g1 ~= g1.edge(2, 3); 805 806 assert(g1.inEdges(1).equal([ 807 g1.edge(1, 1), 808 ])); 809 assert(g1.outEdges(1).equal([ 810 g1.edge(1, 1), 811 g1.edge(1, 2), 812 ])); 813 assert(g1.inEdges(2).equal([ 814 g1.edge(1, 2), 815 g1.edge(2, 2), 816 ])); 817 assert(g1.outEdges(2).equal([ 818 g1.edge(2, 2), 819 g1.edge(2, 3), 820 ])); 821 assert(g1.inEdges(3).equal([ 822 g1.edge(2, 3), 823 ])); 824 assert(g1.outEdges(3).empty); 825 } 826 827 /// Get the in/out degree of node `n`. 828 size_t inDegree(Node n) const nothrow pure 829 { 830 return inEdges(n).walkLength; 831 } 832 833 /// ditto 834 size_t outDegree(Node n) const nothrow pure 835 { 836 return outEdges(n).walkLength; 837 } 838 839 /// 840 unittest 841 { 842 auto g1 = Graph!(int, void, Yes.isDirected)([1, 2, 3]); 843 844 g1 ~= g1.edge(1, 1); 845 g1 ~= g1.edge(1, 2); 846 g1 ~= g1.edge(2, 2); 847 g1 ~= g1.edge(2, 3); 848 849 assert(g1.inDegree(1) == 1); 850 assert(g1.outDegree(1) == 2); 851 assert(g1.inDegree(2) == 2); 852 assert(g1.outDegree(2) == 2); 853 assert(g1.inDegree(3) == 1); 854 assert(g1.outDegree(3) == 0); 855 } 856 } 857 else 858 { 859 /// Returns a range of all edges incident to node `n`. 860 auto incidentEdges(Node n) nothrow pure 861 { 862 return _edges.data[].filter!(e => e.start == n || e.end == n); 863 } 864 865 /// ditto 866 auto incidentEdges(Node n) const nothrow pure 867 { 868 return edges[].filter!(e => e.start == n || e.end == n); 869 } 870 871 /// ditto 872 alias inEdges = incidentEdges; 873 874 /// ditto 875 alias outEdges = incidentEdges; 876 877 /// 878 unittest 879 { 880 import std.algorithm : equal; 881 882 auto g1 = Graph!int([1, 2, 3]); 883 884 g1 ~= g1.edge(1, 1); 885 g1 ~= g1.edge(1, 2); 886 g1 ~= g1.edge(2, 2); 887 g1 ~= g1.edge(2, 3); 888 889 assert(g1.incidentEdges(1).equal([ 890 g1.edge(1, 1), 891 g1.edge(1, 2), 892 ])); 893 assert(g1.incidentEdges(2).equal([ 894 g1.edge(1, 2), 895 g1.edge(2, 2), 896 g1.edge(2, 3), 897 ])); 898 assert(g1.incidentEdges(3).equal([ 899 g1.edge(2, 3), 900 ])); 901 } 902 903 IncidentEdgesCache allIncidentEdges() 904 { 905 return IncidentEdgesCache(this); 906 } 907 908 static struct IncidentEdgesCache 909 { 910 alias G = Graph!(Node, Weight, isDirected, EdgePayload); 911 G graph; 912 Edge[][] incidentEdges; 913 914 this(G graph) 915 { 916 this.graph = graph; 917 collectAllIncidentEdges(); 918 } 919 920 private void collectAllIncidentEdges() 921 { 922 preallocateMemory(); 923 924 size_t startIdx; 925 size_t endIdx; 926 foreach (edge; graph._edges.data) 927 { 928 if (graph._nodes[startIdx] < edge.start) 929 endIdx = startIdx; 930 while (graph._nodes[startIdx] < edge.start) 931 ++startIdx; 932 if (endIdx < startIdx) 933 endIdx = startIdx; 934 while (graph._nodes[endIdx] < edge.end) 935 ++endIdx; 936 937 incidentEdges[startIdx] ~= edge; 938 // Avoid double-counting of loops 939 if (startIdx != endIdx) 940 incidentEdges[endIdx] ~= edge; 941 } 942 } 943 944 void preallocateMemory() 945 { 946 auto degreesCache = graph.allDegrees(); 947 Edge[] buffer; 948 buffer.length = degreesCache.degrees.sum; 949 incidentEdges.length = degreesCache.degrees.length; 950 951 size_t sliceBegin; 952 size_t startIdx; 953 foreach (degree; degreesCache) 954 { 955 incidentEdges[startIdx] = buffer[sliceBegin .. sliceBegin + degree]; 956 incidentEdges[startIdx].length = 0; 957 958 sliceBegin += degree; 959 ++startIdx; 960 } 961 } 962 963 ref inout(Edge[]) opIndex(in Node node) inout 964 { 965 return incidentEdges[graph.indexOf(node)]; 966 } 967 968 ref inout(Edge[]) opIndex(in size_t nodeIdx) inout 969 { 970 return incidentEdges[nodeIdx]; 971 } 972 973 int opApply(scope int delegate(Edge[]) yield) 974 { 975 int result = 0; 976 977 foreach (currentIncidentEdges; incidentEdges) 978 { 979 result = yield(currentIncidentEdges); 980 if (result) 981 break; 982 } 983 984 return result; 985 } 986 987 int opApply(scope int delegate(Node, Edge[]) yield) 988 { 989 int result = 0; 990 991 foreach (i, currentIncidentEdges; incidentEdges) 992 { 993 result = yield(graph._nodes[i], currentIncidentEdges); 994 if (result) 995 break; 996 } 997 998 return result; 999 } 1000 } 1001 1002 /// Get the `adjacencyList` of this graph where nodes are represented 1003 /// by their index in the nodes list. 1004 size_t[][] adjacencyList() const 1005 { 1006 size_t[][] _adjacencyList; 1007 _adjacencyList.length = nodes.length; 1008 size_t[] targetsBuffer; 1009 targetsBuffer.length = 2 * edges.length; 1010 1011 foreach (i, node; _nodes) 1012 { 1013 auto bufferRest = edges 1014 .filter!(e => e.start == node || e.end == node) 1015 .map!(edge => indexOf(edge.target(node))) 1016 .copy(targetsBuffer); 1017 _adjacencyList[i] = targetsBuffer[0 .. $ - bufferRest.length]; 1018 _adjacencyList[i].sort; 1019 targetsBuffer = bufferRest; 1020 } 1021 1022 return _adjacencyList; 1023 } 1024 1025 /// 1026 unittest 1027 { 1028 auto g1 = Graph!int([1, 2, 3, 4]); 1029 1030 g1 ~= g1.edge(1, 1); 1031 g1 ~= g1.edge(1, 2); 1032 g1 ~= g1.edge(2, 2); 1033 g1 ~= g1.edge(2, 3); 1034 g1 ~= g1.edge(2, 4); 1035 g1 ~= g1.edge(3, 4); 1036 1037 assert(g1.adjacencyList() == [ 1038 [0, 1], 1039 [0, 1, 2, 3], 1040 [1, 3], 1041 [1, 2], 1042 ]); 1043 } 1044 1045 /// Get the degree of node `n`. 1046 size_t degree(Node n) const nothrow pure 1047 { 1048 return incidentEdges(n).walkLength; 1049 } 1050 1051 /// ditto 1052 alias inDegree = degree; 1053 1054 /// ditto 1055 alias outDegree = degree; 1056 1057 DegreesCache allDegrees() const 1058 { 1059 return DegreesCache(this); 1060 } 1061 1062 static struct DegreesCache 1063 { 1064 alias G = Graph!(Node, Weight, isDirected, EdgePayload); 1065 const(G) graph; 1066 size_t[] degrees; 1067 1068 this(in G graph) 1069 { 1070 this.graph = graph; 1071 collectAllDegrees(); 1072 } 1073 1074 private void collectAllDegrees() 1075 { 1076 degrees.length = graph._nodes.length; 1077 1078 size_t startIdx; 1079 size_t endIdx; 1080 foreach (edge; graph._edges.data) 1081 { 1082 if (graph._nodes[startIdx] < edge.start) 1083 endIdx = startIdx; 1084 while (graph._nodes[startIdx] < edge.start) 1085 ++startIdx; 1086 if (endIdx < startIdx) 1087 endIdx = startIdx; 1088 while (graph._nodes[endIdx] < edge.end) 1089 ++endIdx; 1090 1091 ++degrees[startIdx]; 1092 // Avoid double-counting of loops 1093 if (startIdx != endIdx) 1094 ++degrees[endIdx]; 1095 } 1096 } 1097 1098 size_t opIndex(in Node node) const 1099 { 1100 return degrees[graph.indexOf(node)]; 1101 } 1102 1103 size_t opIndex(in size_t nodeIdx) const 1104 { 1105 return degrees[nodeIdx]; 1106 } 1107 1108 int opApply(scope int delegate(size_t) yield) const 1109 { 1110 int result = 0; 1111 1112 foreach (degree; degrees) 1113 { 1114 result = yield(degree); 1115 if (result) 1116 break; 1117 } 1118 1119 return result; 1120 } 1121 1122 int opApply(scope int delegate(Node, size_t) yield) const 1123 { 1124 int result = 0; 1125 1126 foreach (i, degree; degrees) 1127 { 1128 result = yield(graph._nodes[i], degree); 1129 if (result) 1130 break; 1131 } 1132 1133 return result; 1134 } 1135 } 1136 } 1137 } 1138 1139 /// 1140 unittest 1141 { 1142 // +-+ +-+ 1143 // \ / \ / 1144 // (1)--(2) 1145 auto g1 = Graph!int([1, 2]); 1146 1147 g1 ~= g1.edge(1, 1); 1148 g1 ~= g1.edge(1, 2); 1149 g1.add(g1.edge(2, 2)); 1150 1151 assert(g1.edge(1, 1) in g1); 1152 assert(g1.edge(1, 2) in g1); 1153 assert(g1.edge(2, 1) in g1); 1154 assert(g1.has(g1.edge(2, 2))); 1155 assert(g1.allDegrees().degrees == [2, 2]); 1156 assert(g1.allIncidentEdges().incidentEdges == [ 1157 [g1.edge(1, 1), g1.edge(1, 2)], 1158 [g1.edge(1, 2), g1.edge(2, 2)], 1159 ]); 1160 1161 // 0.5 0.5 1162 // +-+ +-+ 1163 // \ / \ / 1164 // (1)-----(2) 1165 // 1.0 1166 auto g2 = Graph!(int, double)([1, 2]); 1167 1168 g2 ~= g2.edge(1, 1, 0.5); 1169 g2 ~= g2.edge(1, 2, 1.0); 1170 g2.add(g2.edge(2, 2, 0.5)); 1171 1172 assert(g2.edge(1, 1) in g2); 1173 assert(g2.edge(1, 2) in g2); 1174 assert(g2.edge(2, 1) in g2); 1175 assert(g2.has(g2.edge(2, 2))); 1176 assert(g2.allDegrees().degrees == [2, 2]); 1177 assert(g2.allIncidentEdges().incidentEdges == [ 1178 [g2.edge(1, 1, 0.5), g2.edge(1, 2, 1.0)], 1179 [g2.edge(1, 2, 1.0), g2.edge(2, 2, 0.5)], 1180 ]); 1181 1182 // 0.5 0.5 1183 // +-+ +-+ 1184 // \ v v / 1185 // (1)---->(2) 1186 // 1.0 1187 auto g3 = Graph!(int, double, Yes.isDirected)([1, 2]); 1188 1189 g3 ~= g3.edge(1, 1, 0.5); 1190 g3 ~= g3.edge(1, 2, 1.0); 1191 g3.add(g3.edge(2, 2, 0.5)); 1192 1193 assert(g3.edge(1, 1) in g3); 1194 assert(g3.edge(1, 2) in g3); 1195 assert(!(g3.edge(2, 1) in g3)); 1196 assert(g3.has(g3.edge(2, 2))); 1197 1198 // +-+ +-+ 1199 // \ v v / 1200 // (1)-->(2) 1201 auto g4 = Graph!(int, void, Yes.isDirected)([1, 2]); 1202 1203 g4 ~= g4.edge(1, 1); 1204 g4 ~= g4.edge(1, 2); 1205 g4.add(g4.edge(2, 2)); 1206 1207 assert(g4.edge(1, 1) in g4); 1208 assert(g4.edge(1, 2) in g4); 1209 assert(!(g4.edge(2, 1) in g4)); 1210 assert(g4.has(g4.edge(2, 2))); 1211 1212 // +-+ +-+ 1213 // \ / \ / 1214 // (1)--(2) 1215 // 1216 // payload(1, 1) = [1]; 1217 // payload(1, 2) = [2]; 1218 // payload(2, 2) = [3]; 1219 auto g5 = Graph!(int, void, No.isDirected, int[])([1, 2]); 1220 1221 g5 ~= g5.edge(1, 1, [1]); 1222 g5 ~= g5.edge(1, 2, [2]); 1223 g5.add(g5.edge(2, 2, [3])); 1224 1225 assert(g5.edge(1, 1) in g5); 1226 assert(g5.get(g5.edge(1, 1)).payload == [1]); 1227 assert(g5.edge(1, 2) in g5); 1228 assert(g5.get(g5.edge(1, 2)).payload == [2]); 1229 assert(g5.edge(2, 1) in g5); 1230 assert(g5.get(g5.edge(2, 1)).payload == [2]); 1231 assert(g5.has(g5.edge(2, 2))); 1232 assert(g5.get(g5.edge(2, 2)).payload == [3]); 1233 assert(g5.allDegrees().degrees == [2, 2]); 1234 assert(g5.allIncidentEdges().incidentEdges == [ 1235 [g5.edge(1, 1), g5.edge(1, 2)], 1236 [g5.edge(1, 2), g5.edge(2, 2)], 1237 ]); 1238 } 1239 1240 /// 1241 unittest 1242 { 1243 // -1 1 1 1244 // (1)----(2)---(3) (4)---(5) (6) 1245 size_t[] contigs = [1, 2, 3, 4, 5, 6]; 1246 auto contigGraph = Graph!(size_t, int)([1, 2, 3, 4, 5, 6]); 1247 1248 contigGraph.add(contigGraph.edge(1, 2, -1)); 1249 contigGraph.add(contigGraph.edge(2, 3, 1)); 1250 contigGraph.add(contigGraph.edge(4, 5, 1)); 1251 1252 foreach (contig; contigs) 1253 { 1254 assert(contigGraph.degree(contig) <= 2); 1255 } 1256 assert(contigGraph.allDegrees().degrees == [1, 2, 1, 1, 1, 0]); 1257 assert(contigGraph.allIncidentEdges().incidentEdges == [ 1258 [contigGraph.edge(1, 2, -1)], 1259 [contigGraph.edge(1, 2, -1), contigGraph.edge(2, 3, 1)], 1260 [contigGraph.edge(2, 3, 1)], 1261 [contigGraph.edge(4, 5, 1)], 1262 [contigGraph.edge(4, 5, 1)], 1263 [], 1264 ]); 1265 } 1266 1267 /// Add a set of edges to this graph and merge mutli-edges using `merge`. 1268 void bulkAdd(alias merge, G, R)(ref G graph, R edges) 1269 if (is(G : Graph!Params, Params...) && isForwardRange!R && is(ElementType!R == G.Edge)) 1270 { 1271 alias Edge = G.Edge; 1272 alias ReturnTypeMerge = typeof(merge(new Edge[0])); 1273 static assert(is(ReturnTypeMerge == Edge), "expected `Edge merge(Edge[] multiEdge)`"); 1274 1275 graph.bulkAddForce(edges); 1276 1277 auto bufferRest = graph 1278 ._edges 1279 .data 1280 .sliceBy!(G.groupByNodes) 1281 .map!(unaryFun!merge) 1282 .copy(graph._edges.data); 1283 graph._edges.shrinkTo(graph._edges.data.length - bufferRest.length); 1284 } 1285 1286 /// 1287 unittest 1288 { 1289 auto g1 = Graph!(int, int)([1, 2]); 1290 1291 static g1.Edge sumWeights(g1.Edge[] multiEdge) 1292 { 1293 auto sumOfWeights = multiEdge.map!"a.weight".sum; 1294 auto mergedEdge = multiEdge[0]; 1295 mergedEdge.weight = sumOfWeights; 1296 1297 return mergedEdge; 1298 } 1299 1300 auto edges = [ 1301 g1.edge(1, 2, 1), 1302 g1.edge(1, 2, 1), 1303 g1.edge(1, 2, 1), 1304 g1.edge(2, 3, 2), 1305 g1.edge(2, 3, 2), 1306 g1.edge(3, 4, 3), 1307 ]; 1308 g1.bulkAdd!sumWeights(edges); 1309 assert(g1.edges == [ 1310 g1.edge(1, 2, 3), 1311 g1.edge(2, 3, 4), 1312 g1.edge(3, 4, 3), 1313 ]); 1314 } 1315 1316 /// Add an edge to this graph and handle existing edges with `handleConflict`. 1317 /// The handler must have this signature `Edge handleConflict(Edge, Edge)`. 1318 G.Edge add(alias handleConflict = 1337, G)(ref G graph, G.Edge edge) 1319 if (is(G : Graph!Params, Params...)) 1320 { 1321 static if (isCallable!handleConflict) 1322 alias handleConflict_ = binaryFun!handleConflict; 1323 else 1324 alias handleConflict_ = binaryFun!(G.ConflictStrategy.error); 1325 1326 if (!graph.has(edge.start) || !graph.has(edge.end)) 1327 { 1328 throw new MissingNodeException(); 1329 } 1330 1331 auto sortedEdges = assumeSorted!(G.orderByNodes)(graph._edges.data); 1332 auto trisectedEdges = sortedEdges.trisect(edge); 1333 auto existingEdges = trisectedEdges[1]; 1334 auto existingEdgeIdx = trisectedEdges[0].length; 1335 1336 if (existingEdges.empty) 1337 { 1338 return graph.forceAdd(edge); 1339 } 1340 else 1341 { 1342 auto newEdge = handleConflict_(existingEdges.front, edge); 1343 1344 return graph.replaceEdge(existingEdgeIdx, newEdge); 1345 } 1346 } 1347 1348 /// 1349 unittest 1350 { 1351 auto g1 = Graph!(int, int)([1, 2]); 1352 1353 auto e1 = g1.edge(1, 2, 1); 1354 auto e2 = g1.edge(1, 2, 2); 1355 1356 g1 ~= e1; 1357 1358 assertThrown!EdgeExistsException(g1.add(e2)); 1359 1360 with (g1.ConflictStrategy) 1361 { 1362 g1.add!replace(e2); 1363 1364 assert(g1.get(g1.edge(1, 2)) == e2); 1365 1366 g1.add!keep(e1); 1367 1368 assert(g1.get(g1.edge(1, 2)) == e2); 1369 1370 g1.add!sumWeights(e2); 1371 1372 assert(g1.get(g1.edge(1, 2)).weight == 2 * e2.weight); 1373 } 1374 } 1375 1376 void filterEdges(alias pred, G)(ref G graph) if (is(G : Graph!Params, Params...)) 1377 { 1378 auto bufferRest = graph 1379 ._edges 1380 .data 1381 .filter!pred 1382 .copy(graph._edges.data); 1383 graph._edges.shrinkTo(graph._edges.data.length - bufferRest.length); 1384 } 1385 1386 class EmptySetException : Exception 1387 { 1388 this(string msg) 1389 { 1390 super(msg); 1391 } 1392 } 1393 1394 struct NaturalNumberSet 1395 { 1396 private static enum partSize = 8 * size_t.sizeof; 1397 private static enum size_t firstBit = 1; 1398 private static enum size_t lastBit = firstBit << (partSize - 1); 1399 private static enum size_t emptyPart = 0; 1400 private static enum size_t fullPart = ~emptyPart; 1401 1402 private size_t[] parts; 1403 private size_t nMax; 1404 1405 this(size_t initialNumElements, Flag!"addAll" addAll = No.addAll) 1406 { 1407 reserveFor(initialNumElements); 1408 1409 if (addAll) 1410 { 1411 foreach (i; 0 .. initialNumElements / partSize) 1412 parts[i] = fullPart; 1413 foreach (i; initialNumElements / partSize .. initialNumElements) 1414 add(i); 1415 } 1416 } 1417 1418 this(this) 1419 { 1420 parts = parts.dup; 1421 } 1422 1423 private this(size_t[] parts) 1424 { 1425 this.parts = parts; 1426 } 1427 1428 private bool inBounds(in size_t n) const pure nothrow 1429 { 1430 return n < nMax; 1431 } 1432 1433 void reserveFor(in size_t n) 1434 { 1435 if (parts.length == 0) 1436 { 1437 parts.length = max(1, ceil(n, partSize) / partSize); 1438 } 1439 1440 while (!inBounds(n)) 1441 { 1442 parts.length *= 2; 1443 nMax = parts.length * partSize; 1444 } 1445 } 1446 1447 @property size_t capacity() pure const nothrow 1448 { 1449 return nMax; 1450 } 1451 1452 private size_t partIdx(in size_t n) const pure nothrow 1453 { 1454 return n / partSize; 1455 } 1456 1457 private size_t idxInPart(in size_t n) const pure nothrow 1458 { 1459 return n % partSize; 1460 } 1461 1462 private size_t itemMask(in size_t n) const pure nothrow 1463 { 1464 return firstBit << idxInPart(n); 1465 } 1466 1467 static size_t inverse(in size_t n) pure nothrow 1468 { 1469 return n ^ fullPart; 1470 } 1471 1472 void add(in size_t n) 1473 { 1474 reserveFor(n); 1475 1476 parts[partIdx(n)] |= itemMask(n); 1477 } 1478 1479 void remove(in size_t n) 1480 { 1481 if (!inBounds(n)) 1482 { 1483 return; 1484 } 1485 1486 parts[partIdx(n)] &= inverse(itemMask(n)); 1487 } 1488 1489 bool has(in size_t n) const pure nothrow 1490 { 1491 if (!inBounds(n)) 1492 { 1493 return false; 1494 } 1495 1496 return (parts[partIdx(n)] & itemMask(n)) != emptyPart; 1497 } 1498 1499 bool empty() const pure nothrow 1500 { 1501 return parts.all!(part => part == emptyPart); 1502 } 1503 1504 size_t minElement() const 1505 { 1506 foreach (i, part; parts) 1507 { 1508 if (part != emptyPart) 1509 { 1510 size_t j = 0; 1511 1512 while (((part >> j) & firstBit) != firstBit) 1513 { 1514 ++j; 1515 } 1516 1517 return i * partSize + j; 1518 } 1519 } 1520 1521 throw new EmptySetException("empty set has no minElement"); 1522 } 1523 1524 size_t maxElement() const 1525 { 1526 foreach (i, part; parts.retro.enumerate) 1527 { 1528 if (part != emptyPart) 1529 { 1530 size_t j = 0; 1531 1532 while (((part << j) & lastBit) != lastBit) 1533 { 1534 ++j; 1535 } 1536 1537 return (parts.length - i - 1) * partSize + (partSize - j - 1); 1538 } 1539 } 1540 1541 throw new EmptySetException("empty set has no maxElement"); 1542 } 1543 1544 unittest 1545 { 1546 foreach (i; 0 .. 2 * NaturalNumberSet.partSize) 1547 { 1548 NaturalNumberSet set; 1549 1550 set.add(i + 5); 1551 set.add(i + 7); 1552 1553 assert(set.minElement() == i + 5); 1554 assert(set.maxElement() == i + 7); 1555 } 1556 } 1557 1558 /// Returns a range of the elements in this set. The elements are ordered 1559 /// ascending. 1560 @property auto elements() const pure nothrow 1561 { 1562 static struct ElementsRange 1563 { 1564 const NaturalNumberSet* set; 1565 bool _empty = false; 1566 size_t i = 0; 1567 size_t part; 1568 size_t j = 0; 1569 1570 this(const NaturalNumberSet* set) pure nothrow 1571 { 1572 this.set = set; 1573 this._empty = set.empty; 1574 1575 if (!this.empty) 1576 { 1577 this.part = set.parts[i]; 1578 if (!set.has(front)) 1579 { 1580 popFront(); 1581 } 1582 } 1583 } 1584 1585 void popFront() pure nothrow 1586 { 1587 assert(!empty, "Attempting to popFront an empty elements range"); 1588 ++j; 1589 1590 while (shiftedPartEmpty) 1591 { 1592 nextPart(); 1593 1594 if (empty) 1595 { 1596 return; 1597 } 1598 } 1599 1600 while (((part >> j) & firstBit) != firstBit && !shiftedPartEmpty) 1601 { 1602 ++j; 1603 } 1604 1605 if (shiftedPartEmpty) 1606 { 1607 popFront(); 1608 } 1609 } 1610 1611 @property size_t front() const pure nothrow 1612 { 1613 assert(!empty, "Attempting to fetch the front of an empty elements range"); 1614 return i * partSize + j; 1615 } 1616 1617 @property bool empty() const pure nothrow 1618 { 1619 return _empty; 1620 } 1621 1622 private @property bool shiftedPartEmpty() const pure nothrow 1623 { 1624 return (part >> j) == emptyPart || j >= partSize; 1625 } 1626 1627 private void nextPart() pure nothrow 1628 { 1629 // move to start of next part 1630 ++i; 1631 j = 0; 1632 1633 if (i < set.parts.length) 1634 { 1635 part = set.parts[i]; 1636 } 1637 else 1638 { 1639 _empty = true; 1640 } 1641 } 1642 } 1643 1644 return ElementsRange(&this); 1645 } 1646 1647 /// 1648 unittest 1649 { 1650 import std.algorithm : equal; 1651 import std.range : iota; 1652 1653 NaturalNumberSet set; 1654 auto someNumbers = iota(set.partSize).filter!"a % 3 == 0"; 1655 1656 foreach (i; someNumbers) 1657 { 1658 set.add(i); 1659 } 1660 1661 assert(equal(someNumbers, set.elements)); 1662 } 1663 } 1664 1665 unittest 1666 { 1667 NaturalNumberSet set; 1668 1669 // add some numbers 1670 foreach (i; 0 .. set.partSize) 1671 { 1672 if (i % 2 == 0) 1673 { 1674 set.add(i); 1675 } 1676 } 1677 1678 // force extension of set 1679 foreach (i; set.partSize .. 2 * set.partSize) 1680 { 1681 if (i % 3 == 0) 1682 { 1683 set.add(i); 1684 } 1685 } 1686 1687 // validate presence 1688 foreach (i; 0 .. 2 * set.partSize) 1689 { 1690 if (i / set.partSize == 0 && i % 2 == 0) 1691 { 1692 assert(set.has(i)); 1693 } 1694 else if (i / set.partSize == 1 && i % 3 == 0) 1695 { 1696 assert(set.has(i)); 1697 } 1698 else 1699 { 1700 assert(!set.has(i)); 1701 } 1702 } 1703 } 1704 1705 /** 1706 Find all maximal cliques in a graph represented by `adjacencyList`. 1707 The implementation is based on version 1 of the Bron-Kerbosch algorithm [1]. 1708 1709 [1]: Bron, C.; Kerbosch, J. (1973), "Algorithm 457: finding all cliques 1710 of an undirected graph", Communications of the ACM, 16 (9): 575–577, 1711 doi:10.1145/362342.362367. 1712 1713 Returns: list of sets of nodes each representing a maximal clique 1714 */ 1715 auto findAllCliques(in size_t[][] adjacencyList) 1716 { 1717 return BronKerboschVersion1(adjacencyList); 1718 } 1719 1720 /// 1721 unittest 1722 { 1723 auto g = Graph!int([0, 1, 2, 3, 4, 5, 6, 7, 8, 9]); 1724 g.add(g.edge(0, 1)); 1725 g.add(g.edge(0, 2)); 1726 g.add(g.edge(1, 2)); 1727 g.add(g.edge(1, 7)); 1728 g.add(g.edge(1, 8)); 1729 g.add(g.edge(2, 3)); 1730 g.add(g.edge(3, 4)); 1731 g.add(g.edge(3, 5)); 1732 g.add(g.edge(3, 6)); 1733 g.add(g.edge(4, 5)); 1734 g.add(g.edge(4, 6)); 1735 g.add(g.edge(5, 6)); 1736 g.add(g.edge(6, 7)); 1737 g.add(g.edge(7, 8)); 1738 1739 auto cliques = array(findAllCliques(g.adjacencyList())); 1740 1741 assert(cliques == [ 1742 [0, 1, 2], 1743 [1, 7, 8], 1744 [2, 3], 1745 [3, 4, 5, 6], 1746 [6, 7], 1747 [9], 1748 ]); 1749 } 1750 1751 private struct BronKerboschVersion1 1752 { 1753 const size_t[][] adjacencyList; 1754 1755 int opApply(scope int delegate(size_t[]) yield) 1756 { 1757 size_t[] clique; 1758 clique.reserve(adjacencyList.length); 1759 1760 auto candidates = NaturalNumberSet(adjacencyList.length, Yes.addAll); 1761 auto not = NaturalNumberSet(adjacencyList.length); 1762 1763 return extendClique(clique, candidates, not, yield); 1764 } 1765 1766 private int extendClique( 1767 size_t[] clique, 1768 NaturalNumberSet candidates, 1769 NaturalNumberSet not, 1770 scope int delegate(size_t[]) yield, 1771 ) 1772 { 1773 import std.stdio; 1774 1775 if (not.empty && candidates.empty) 1776 return clique.length == 0 ? 0 : yield(clique); 1777 1778 int result; 1779 1780 foreach (candidate; candidates.elements) 1781 { 1782 clique ~= candidate; 1783 1784 auto reducedCandidates = NaturalNumberSet(adjacencyList.length); 1785 auto reducedNot = NaturalNumberSet(adjacencyList.length); 1786 1787 foreach (neighbourNode; adjacencyList[candidate]) 1788 { 1789 if (candidates.has(neighbourNode)) 1790 reducedCandidates.add(neighbourNode); 1791 if (not.has(neighbourNode)) 1792 reducedNot.add(neighbourNode); 1793 } 1794 1795 result = extendClique(clique, reducedCandidates, reducedNot, yield); 1796 1797 if (result) 1798 return result; 1799 1800 candidates.remove(candidate); 1801 not.add(candidate); 1802 --clique.length; 1803 } 1804 1805 return result; 1806 } 1807 } 1808 1809 /** 1810 Calculate a longest increasing subsequence of `sequence`. This subsequence 1811 is not necessarily contiguous, or unique. Given a `sequence` of `n` 1812 elements the algorithm uses `O(n log n)` evaluation of `pred`. 1813 1814 See_Also: https://en.wikipedia.org/wiki/Longest_increasing_subsequence 1815 */ 1816 auto longestIncreasingSubsequence(alias pred = "a < b", Range)(Range sequence) 1817 if (isRandomAccessRange!Range) 1818 { 1819 alias lessThan = binaryFun!pred; 1820 1821 size_t[] subseqEnds; 1822 subseqEnds.length = sequence.length; 1823 size_t[] predecessors; 1824 predecessors.length = sequence.length; 1825 size_t subseqLength; 1826 1827 foreach (i; 0 .. sequence.length) 1828 { 1829 // Binary search for the largest positive j < subseqLength 1830 // such that sequence[subseqEnds[j]] < sequence[i] 1831 long lo = 0; 1832 long hi = subseqLength - 1; 1833 auto pivot = sequence[i]; 1834 assert(!lessThan(pivot, pivot), "`pred` is not anti-symmetric"); 1835 1836 while (lo <= hi) 1837 { 1838 auto mid = ceildiv(lo + hi, 2); 1839 1840 if (lessThan(sequence[subseqEnds[mid]], pivot)) 1841 lo = mid + 1; 1842 else 1843 hi = mid - 1; 1844 } 1845 1846 // After searching, lo + 1 is the length of the longest prefix of 1847 // sequence[i] 1848 auto newSubseqLength = lo + 1; 1849 1850 // The predecessor of sequence[i] is the last index of 1851 // the subsequence of length newSubseqLength - 1 1852 subseqEnds[lo] = i; 1853 if (lo > 0) 1854 predecessors[i] = subseqEnds[lo - 1]; 1855 1856 if (newSubseqLength > subseqLength) 1857 // If we found a subsequence longer than any we've 1858 // found yet, update subseqLength 1859 subseqLength = newSubseqLength; 1860 } 1861 1862 auto subsequenceResult = subseqEnds[0 .. subseqLength]; 1863 1864 if (subseqLength > 0) 1865 { 1866 // Reconstruct the longest increasing subsequence 1867 // Note: reusing memory from now unused subseqEnds 1868 auto k = subseqEnds[subseqLength - 1]; 1869 foreach_reverse (i; 0 .. subseqLength) 1870 { 1871 subsequenceResult[i] = k; 1872 k = predecessors[k]; 1873 } 1874 } 1875 1876 return subsequenceResult.map!(i => sequence[i]); 1877 } 1878 1879 /// Example from Wikipedia 1880 unittest 1881 { 1882 import std.algorithm : equal; 1883 1884 auto inputSequence = [0, 8, 4, 12, 2, 10, 6, 14, 1, 9, 5, 13, 3, 11, 7, 15]; 1885 auto expectedOutput = [0, 2, 6, 9, 11, 15]; 1886 1887 assert(inputSequence.longestIncreasingSubsequence.equal(expectedOutput)); 1888 } 1889 1890 /// Example using a different `pred` 1891 unittest 1892 { 1893 import std.algorithm : equal; 1894 import std.range : retro; 1895 1896 auto inputSequence = [0, 8, 4, 12, 2, 10, 6, 14, 1, 9, 5, 13, 3, 11, 7, 15]; 1897 auto expectedOutput = [12, 10, 9, 5, 3]; 1898 1899 assert(inputSequence.longestIncreasingSubsequence!"a > b".equal(expectedOutput)); 1900 } 1901 1902 unittest 1903 { 1904 import std.algorithm : equal; 1905 1906 int[] inputSequence = []; 1907 int[] expectedOutput = []; 1908 1909 assert(inputSequence.longestIncreasingSubsequence.equal(expectedOutput)); 1910 } 1911 1912 unittest 1913 { 1914 import std.algorithm : equal; 1915 1916 auto inputSequence = [1, 2, 3, 4, 5]; 1917 auto expectedOutput = [1, 2, 3, 4, 5]; 1918 1919 assert(inputSequence.longestIncreasingSubsequence.equal(expectedOutput)); 1920 } 1921 1922 unittest 1923 { 1924 import std.algorithm : equal; 1925 1926 auto inputSequence = [2, 1, 3, 4, 5]; 1927 auto expectedOutput = [1, 3, 4, 5]; 1928 1929 assert(inputSequence.longestIncreasingSubsequence.equal(expectedOutput)); 1930 } 1931 1932 unittest 1933 { 1934 import std.algorithm : equal; 1935 1936 auto inputSequence = [1, 2, 3, 5, 4]; 1937 auto expectedOutput = [1, 2, 3, 4]; 1938 1939 assert(inputSequence.longestIncreasingSubsequence.equal(expectedOutput)); 1940 }