00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036 #ifndef _DIGRAPH_H_
00037 #define _DIGRAPH_H_
00038
00039 #include <iostream>
00040 #include <new>
00041 #include <stack>
00042 #include <set>
00043 #include <memory>
00044
00045 #include "chomp/struct/multitab.h"
00046 #include "chomp/struct/flatmatr.h"
00047 #include "chomp/struct/bitfield.h"
00048 #include "chomp/struct/bitsets.h"
00049 #include "chomp/struct/fibheap.h"
00050 #include "chomp/system/textfile.h"
00051 #include "chomp/system/timeused.h"
00052
00053 namespace chomp {
00054 namespace homology {
00055
00056
00057
00058
00059
00060
00061
00062
00063 template <class numType>
00064 class dummyRounding
00065 {
00066 public:
00067
00068 static inline numType add_down (const numType &x, const numType &y)
00069 { return x + y; }
00070
00071
00072 static inline numType add_up (const numType &x, const numType &y)
00073 { return x + y; }
00074
00075
00076 static inline numType sub_down (const numType &x, const numType &y)
00077 { return x - y; }
00078
00079
00080 static inline numType sub_up (const numType &x, const numType &y)
00081 { return x - y; }
00082
00083
00084 static inline numType mul_down (const numType &x, const numType &y)
00085 { return x * y; }
00086
00087
00088 static inline numType mul_up (const numType &x, const numType &y)
00089 { return x * y; }
00090
00091
00092 static inline numType div_down (const numType &x, const numType &y)
00093 { return x / y; }
00094
00095
00096 static inline numType div_down (const numType &x, int y)
00097 { return x / y; }
00098
00099
00100 static inline numType div_up (const numType &x, const numType &y)
00101 { return x / y; }
00102
00103 private:
00104 };
00105
00106
00107
00108
00109
00110
00111
00112 class dummyArray
00113 {
00114 public:
00115
00116 dummyArray () {n = 0;}
00117
00118
00119 int &operator [] (int) {return n;}
00120
00121 private:
00122
00123 int n;
00124
00125 };
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137 template <class wType = int>
00138 class diGraph
00139 {
00140 public:
00141
00142 typedef wType weight_type;
00143
00144
00145 diGraph ();
00146
00147
00148 ~diGraph ();
00149
00150
00151 void swap (diGraph<wType> &g);
00152
00153
00154 void addVertex (void);
00155
00156
00157
00158 void addEdge (int target);
00159
00160
00161
00162 void addEdge (int target, const wType &weight);
00163
00164
00165 void setWeight (int vertex, int i, const wType &weight);
00166
00167
00168 void setWeight (int edge, const wType &weight);
00169
00170
00171
00172 template <class Table>
00173 void setWeights (const Table &tab);
00174
00175
00176
00177 void removeVertex (void);
00178
00179
00180
00181
00182
00183 void removeVertex (int vertex, bool updateweights = false);
00184
00185
00186 int countVertices (void) const;
00187
00188
00189 int countEdges (void) const;
00190
00191
00192 int countEdges (int vertex) const;
00193
00194
00195 int getEdge (int vertex, int i) const;
00196
00197
00198 const wType &getWeight (int vertex, int i) const;
00199
00200
00201 const wType &getWeight (int edge) const;
00202
00203
00204
00205 template <class Table>
00206 void getWeights (Table &tab) const;
00207
00208
00209 void writeEdges (int (*table) [2]) const;
00210
00211
00212 template <class wType1>
00213 void transpose (diGraph<wType1> &result,
00214 bool copyweights = false) const;
00215
00216
00217
00218
00219 template <class Table, class wType1>
00220 void subgraph (diGraph<wType1> &result, const Table &tab,
00221 bool copyweights = false) const;
00222
00223
00224
00225
00226 template <class Table, class Color>
00227 void DFScolor (Table &tab, const Color &color, int vertex = 0) const;
00228
00229
00230 template <class Table, class Color>
00231 void DFScolorRecurrent (Table &tab, const Color &color,
00232 int vertex = 0) const;
00233
00234
00235 template <class Table, class Color>
00236 void DFScolorStack (Table &tab, const Color &color,
00237 int vertex = 0) const;
00238
00239
00240
00241 template <class Table>
00242 void DFSfinishTime (Table &tab) const;
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253 template <class Table1, class Table2, class Table3>
00254 int DFSforest (const Table1 &ordered, Table2 &compVertices,
00255 Table3 &compEnds, bool nontrivial = false,
00256 diGraph<wType> *sccGraph = 0) const;
00257
00258
00259
00260
00261
00262 int shortestPath (int source, int destination) const;
00263
00264
00265
00266
00267
00268 int shortestLoop (int origin) const;
00269
00270
00271
00272
00273
00274
00275
00276
00277 template <class lenTable, class weightsType, class roundType>
00278 void Dijkstra (const roundType &rounding, int source, lenTable &len,
00279 weightsType &edgeWeights) const;
00280
00281
00282 template <class lenTable, class roundType>
00283 void Dijkstra (const roundType &rounding, int source, lenTable &len)
00284 const;
00285
00286
00287 template <class lenTable>
00288 void Dijkstra (int source, lenTable &len) const;
00289
00290
00291
00292
00293
00294
00295
00296
00297
00298
00299
00300 template <class lenTable, class predTable, class roundType>
00301 bool BellmanFord (const roundType &rounding, int source,
00302 lenTable &len, wType *infinity, predTable pred) const;
00303
00304
00305 template <class lenTable, class predTable>
00306 bool BellmanFord (int source, lenTable &len, wType *infinity,
00307 predTable pred) const;
00308
00309
00310
00311
00312 template <class roundType>
00313 bool BellmanFord (const roundType &rounding, int source) const;
00314
00315
00316 bool BellmanFord (int source) const;
00317
00318
00319
00320
00321
00322
00323
00324 wType Edmonds () const;
00325
00326
00327 wType EdmondsOld () const;
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339
00340
00341
00342
00343 template <class arrayType, class roundType>
00344 wType FloydWarshall (const roundType &rounding, arrayType &arr,
00345 bool setInfinity = true, bool ignoreNegLoop = false) const;
00346
00347
00348 template <class arrayType>
00349 wType FloydWarshall (arrayType &arr, bool setInfinity = true,
00350 bool ignoreNegLoop = false) const;
00351
00352
00353
00354
00355
00356
00357
00358
00359
00360 template <class arrayType, class roundType>
00361 wType Johnson (const roundType &rounding, arrayType &arr,
00362 bool setInfinity = true, bool ignoreNegLoop = false) const;
00363
00364
00365 template <class arrayType>
00366 wType Johnson (arrayType &arr, bool setInfinity = true,
00367 bool ignoreNegLoop = false) const;
00368
00369
00370
00371
00372
00373
00374
00375
00376
00377
00378 template <class roundType>
00379 wType minPathWeight (const roundType &rounding,
00380 bool ignoreNegLoop = false, int sparseGraph = -1) const;
00381
00382
00383 wType minPathWeight (bool ignoreNegLoop = false,
00384 int sparseGraph = -1) const;
00385
00386
00387
00388
00389
00390 wType minMeanCycleWeight (diGraph<wType> *transposed = 0) const;
00391
00392
00393
00394
00395
00396
00397
00398
00399
00400 template <class roundType>
00401 wType minMeanCycleWeight (const roundType &rounding,
00402 diGraph<wType> *transposed) const;
00403
00404
00405
00406
00407
00408
00409 template <class arrayType, class roundType>
00410 wType minMeanPathWeight (const roundType &rounding,
00411 const arrayType &starting, int n) const;
00412
00413
00414 template <class arrayType>
00415 wType minMeanPathWeight (const arrayType &starting, int n) const;
00416
00417
00418
00419 template <class wType1, class wType2>
00420 friend bool operator == (const diGraph<wType1> &g1,
00421 const diGraph<wType2> &g2);
00422
00423
00424 template <class outType>
00425 outType &show (outType &out, bool showWeights = false) const;
00426
00427 protected:
00428
00429 int nVertices;
00430
00431
00432
00433 multitable<int> edgeEnds;
00434
00435
00436 multitable<int> edges;
00437
00438
00439 multitable<wType> weights;
00440
00441 private:
00442
00443 template <class Table>
00444 void DFSfinishTimeRecurrent (Table &tab, int vertex,
00445 int &counter) const;
00446
00447
00448 template <class Table>
00449 void DFSfinishTimeStack (Table &tab, int vertex,
00450 int &counter) const;
00451
00452
00453
00454 template <class Table1, class Table2>
00455 bool DFSforestRecurrent (Table1 &tab, Table1 &ntab, int vertex,
00456 int treeNumber, int countTrees, Table2 &compVertices,
00457 int &curVertex, diGraph *sccGraph, int *sccEdgeAdded)
00458 const;
00459
00460
00461 template <class Table1, class Table2>
00462 bool DFSforestStack (Table1 &tab, Table1 &ntab, int vertex,
00463 int treeNumber, int countTrees, Table2 &compVertices,
00464 int &curVertex, diGraph *sccGraph, int *sccEdgeAdded)
00465 const;
00466
00467
00468 struct edgeTriple
00469 {
00470
00471 int vert1, vert2;
00472
00473
00474 wType weight;
00475
00476
00477 friend bool operator < (const edgeTriple &x,
00478 const edgeTriple &y)
00479 {
00480 return (x. weight < y. weight);
00481 }
00482 };
00483
00484
00485
00486 class posWeight
00487 {
00488 public:
00489
00490 posWeight ()
00491 {
00492 value = 0;
00493 return;
00494 }
00495
00496
00497 explicit posWeight (const wType &_value)
00498 {
00499 if (_value < 0)
00500 value = 0;
00501 else
00502 value = _value;
00503 return;
00504 }
00505
00506
00507 void setInfinity ()
00508 {
00509 value = -1;
00510 return;
00511 }
00512
00513
00514 bool isInfinity () const
00515 {
00516 return (value < 0);
00517 }
00518
00519
00520 const wType &getValue () const
00521 {
00522 return value;
00523 }
00524
00525
00526 friend bool operator < (const posWeight &x,
00527 const posWeight &y)
00528 {
00529 if (y. isInfinity ())
00530 return !(x. isInfinity ());
00531 else if (x. isInfinity ())
00532 return false;
00533 return (x. value < y. value);
00534 }
00535
00536
00537 friend std::ostream &operator << (std::ostream &out,
00538 const posWeight &x)
00539 {
00540 if (x. isInfinity ())
00541 out << "+oo";
00542 else
00543 out << x. getValue ();
00544 return out;
00545 }
00546
00547 private:
00548
00549 wType value;
00550 };
00551
00552 };
00553
00554
00555
00556 template <class wType>
00557 inline diGraph<wType>::diGraph (): nVertices (0),
00558 edgeEnds (1024), edges (4096), weights (4096)
00559 {
00560 return;
00561 }
00562
00563 template <class wType>
00564 inline diGraph<wType>::~diGraph ()
00565 {
00566 return;
00567 }
00568
00569
00570
00571 template <class wType1, class wType2>
00572 inline bool operator == (const diGraph<wType1> &g1,
00573 const diGraph<wType2> &g2)
00574 {
00575 if (g1. nVertices != g2. nVertices)
00576 return false;
00577 if (!g1. nVertices)
00578 return true;
00579 for (int i = 0; i < g1. nVertices; ++ i)
00580 if (g1. edgeEnds [i] != g2. edgeEnds [i])
00581 return false;
00582 int nEdges = g1. edgeEnds [g1. nVertices - 1];
00583 for (int i = 0; i < nEdges; ++ i)
00584 if (g1. edges [i] != g2. edges [i])
00585 return false;
00586 return true;
00587 }
00588
00589 template <class wType1, class wType2>
00590 inline bool operator != (const diGraph<wType1> &g1,
00591 const diGraph<wType2> &g2)
00592 {
00593 return !(g1 == g2);
00594 }
00595
00596
00597
00598 template <class wType>
00599 inline void diGraph<wType>::swap (diGraph<wType> &g)
00600 {
00601 my_swap (nVertices, g. nVertices);
00602 edgeEnds. swap (g. edgeEnds);
00603 edges. swap (g. edges);
00604 weights. swap (g. weights);
00605 return;
00606 }
00607
00608 template <class wType>
00609 inline void diGraph<wType>::addVertex (void)
00610 {
00611 edgeEnds [nVertices] = nVertices ? edgeEnds [nVertices - 1] : 0;
00612 ++ nVertices;
00613 return;
00614 }
00615
00616 template <class wType>
00617 inline void diGraph<wType>::addEdge (int target)
00618 {
00619 if (!nVertices)
00620 throw "Trying to add an edge to an empty graph.";
00621
00622
00623 if (target < 0)
00624 throw "Trying to add an edge to a negative vertex.";
00625 int &offset = edgeEnds [nVertices - 1];
00626 if (offset + 1 <= 0)
00627 throw "Too many edges in a diGraph (limit = 2,147,483,647).";
00628 edges [offset] = target;
00629 ++ offset;
00630 return;
00631 }
00632
00633 template <class wType>
00634 inline void diGraph<wType>::addEdge (int target, const wType &weight)
00635 {
00636 if (!nVertices)
00637 throw "Trying to add an edge to an empty graph.";
00638
00639
00640 if (target < 0)
00641 throw "Trying to add an edge to a negative vertex.";
00642 int &offset = edgeEnds [nVertices - 1];
00643 if (offset + 1 <= 0)
00644 throw "Too many edges in a diGraph (limit = 2,147,483,647).";
00645 edges [offset] = target;
00646 weights [offset] = weight;
00647 ++ offset;
00648 return;
00649 }
00650
00651 template <class wType>
00652 inline void diGraph<wType>::setWeight (int vertex, int i,
00653 const wType &weight)
00654 {
00655 if (!vertex)
00656 weights [i] = weight;
00657 else
00658 weights [edgeEnds [vertex - 1] + i] = weight;
00659 return;
00660 }
00661
00662 template <class wType>
00663 inline void diGraph<wType>::setWeight (int edge, const wType &weight)
00664 {
00665 weights [edge] = weight;
00666 return;
00667 }
00668
00669 template <class wType> template <class Table>
00670 inline void diGraph<wType>::setWeights (const Table &tab)
00671 {
00672 if (!nVertices)
00673 return;
00674 int nEdges = edgeEnds [nVertices - 1];
00675 for (int i = 0; i < nEdges; ++ i)
00676 weights [i] = tab [i];
00677 return;
00678 }
00679
00680 template <class wType>
00681 inline void diGraph<wType>::removeVertex (void)
00682 {
00683 if (!nVertices)
00684 throw "Trying to remove a vertex from the empty graph.";
00685 -- nVertices;
00686 return;
00687 }
00688
00689 template <class wType>
00690 inline void diGraph<wType>::removeVertex (int vertex, bool updateweights)
00691 {
00692
00693 if ((vertex < 0) || (vertex >= nVertices))
00694 throw "Trying to remove a vertex that is not in the graph.";
00695
00696
00697 int curEdge = 0;
00698 int newEdge = 0;
00699 int skipped = 0;
00700 for (int v = 0; v < nVertices; ++ v)
00701 {
00702
00703 if (!skipped && (v == vertex))
00704 {
00705 curEdge = edgeEnds [v];
00706 ++ skipped;
00707 continue;
00708 }
00709
00710
00711 int maxEdge = edgeEnds [v];
00712 for (; curEdge < maxEdge; ++ curEdge)
00713 {
00714 if (edges [curEdge] == vertex)
00715 continue;
00716 int target = edges [curEdge];
00717 edges [newEdge] = (target < vertex) ? target :
00718 (target - 1);
00719 if (updateweights)
00720 weights [newEdge] = weights [curEdge];
00721 ++ newEdge;
00722 }
00723 edgeEnds [v - skipped] = newEdge;
00724 }
00725
00726
00727 nVertices -= skipped;
00728
00729 return;
00730 }
00731
00732 template <class wType>
00733 inline int diGraph<wType>::countVertices (void) const
00734 {
00735 return nVertices;
00736 }
00737
00738 template <class wType>
00739 inline int diGraph<wType>::countEdges (void) const
00740 {
00741 if (!nVertices)
00742 return 0;
00743 else
00744 return edgeEnds [nVertices - 1];
00745 }
00746
00747 template <class wType>
00748 inline int diGraph<wType>::countEdges (int vertex) const
00749 {
00750 if (!vertex)
00751 return edgeEnds [0];
00752 else
00753 return edgeEnds [vertex] - edgeEnds [vertex - 1];
00754 }
00755
00756 template <class wType>
00757 inline int diGraph<wType>::getEdge (int vertex, int i) const
00758 {
00759 if (!vertex)
00760 return edges [i];
00761 else
00762 return edges [edgeEnds [vertex - 1] + i];
00763 }
00764
00765 template <class wType>
00766 inline const wType &diGraph<wType>::getWeight (int vertex, int i) const
00767 {
00768 if (!vertex)
00769 return weights [i];
00770 else
00771 return weights [edgeEnds [vertex - 1] + i];
00772 }
00773
00774 template <class wType>
00775 inline const wType &diGraph<wType>::getWeight (int edge) const
00776 {
00777 return weights [edge];
00778 }
00779
00780 template <class wType> template <class Table>
00781 inline void diGraph<wType>::getWeights (Table &tab) const
00782 {
00783 if (!nVertices)
00784 return;
00785 int nEdges = edgeEnds [nVertices - 1];
00786 for (int i = 0; i < nEdges; ++ i)
00787 tab [i] = weights [i];
00788 return;
00789 }
00790
00791 template <class wType>
00792 inline void diGraph<wType>::writeEdges (int (*table) [2]) const
00793 {
00794 int curEdge = 0;
00795 for (int curVertex = 0; curVertex < nVertices; ++ curVertex)
00796 {
00797 int maxEdge = edgeEnds [curVertex];
00798 while (curEdge < maxEdge)
00799 {
00800 table [curEdge] [0] = curVertex;
00801 table [curEdge] [1] = edges [curEdge];
00802 ++ curEdge;
00803 }
00804 }
00805 return;
00806 }
00807
00808 template <class wType> template <class wType1>
00809 inline void diGraph<wType>::transpose (diGraph<wType1> &result,
00810 bool copyweights) const
00811 {
00812
00813 result. nVertices = nVertices;
00814 if (!nVertices)
00815 return;
00816
00817
00818 for (int i = 0; i < nVertices; ++ i)
00819 result. edgeEnds [i] = 0;
00820 int nEdges = edgeEnds [nVertices - 1];
00821 for (int i = 0; i < nEdges; ++ i)
00822 if (edges [i] < nVertices - 1)
00823 ++ result. edgeEnds [edges [i] + 1];
00824 for (int i = 2; i < nVertices; ++ i)
00825 result. edgeEnds [i] += result. edgeEnds [i - 1];
00826
00827
00828 int curEdge = 0;
00829 for (int i = 0; i < nVertices; ++ i)
00830 {
00831 for (; curEdge < edgeEnds [i]; ++ curEdge)
00832 {
00833 int j = edges [curEdge];
00834 int &offset = result. edgeEnds [j];
00835 result. edges [offset] = i;
00836 if (copyweights)
00837 result. weights [offset] = weights [curEdge];
00838 ++ offset;
00839 }
00840 }
00841 return;
00842 }
00843
00844 template <class wType> template <class Table, class wType1>
00845 inline void diGraph<wType>::subgraph (diGraph<wType1> &g,
00846 const Table &tab, bool copyweights) const
00847 {
00848
00849 int *numbers = new int [nVertices];
00850 int curNumber = 0;
00851 for (int i = 0; i < nVertices; ++ i)
00852 if (tab [i])
00853 numbers [i] = curNumber ++;
00854 else
00855 numbers [i] = -1;
00856
00857
00858 for (int i = 0; i < nVertices; ++ i)
00859 {
00860 if (numbers [i] < 0)
00861 continue;
00862 g. addVertex ();
00863 int firstEdge = i ? edgeEnds [i - 1] : 0;
00864 int endEdge = edgeEnds [i];
00865 for (int j = firstEdge; j < endEdge; ++ j)
00866 {
00867 int edgeEnd = edges [j];
00868 if (numbers [edgeEnd] >= 0)
00869 {
00870 if (copyweights)
00871 g. addEdge (numbers [edgeEnd],
00872 weights [j]);
00873 else
00874 g. addEdge (numbers [edgeEnd]);
00875 }
00876 }
00877 }
00878
00879
00880 delete [] numbers;
00881 return;
00882 }
00883
00884
00885
00886 template <class wType> template <class Table, class Color>
00887 inline void diGraph<wType>::DFScolorRecurrent (Table &tab,
00888 const Color &color, int vertex) const
00889 {
00890 tab [vertex] = color;
00891 int maxEdge = edgeEnds [vertex];
00892 for (int i = vertex ? edgeEnds [vertex - 1] : 0; i < maxEdge; ++ i)
00893 {
00894 int next = edges [i];
00895 if (tab [next] != color)
00896 DFScolorRecurrent (tab, color, next);
00897 }
00898 return;
00899 }
00900
00901 template <class wType> template <class Table, class Color>
00902 inline void diGraph<wType>::DFScolorStack (Table &tab, const Color &color,
00903 int vertex) const
00904 {
00905
00906 std::stack<int> s_vertex;
00907 std::stack<int> s_edge;
00908 std::stack<int> s_maxedge;
00909
00910
00911 tab [vertex] = color;
00912
00913
00914 int edge = vertex ? edgeEnds [vertex - 1] : 0;
00915 int maxedge = edgeEnds [vertex];
00916
00917 while (1)
00918 {
00919
00920
00921 if (edge >= maxedge)
00922 {
00923
00924 if (s_vertex. empty ())
00925 return;
00926
00927
00928 vertex = s_vertex. top ();
00929 s_vertex. pop ();
00930 edge = s_edge. top ();
00931 s_edge. pop ();
00932 maxedge = s_maxedge. top ();
00933 s_maxedge. pop ();
00934 continue;
00935 }
00936
00937
00938 int next = edges [edge ++];
00939 if (tab [next] != color)
00940 {
00941
00942 s_vertex. push (vertex);
00943 s_edge. push (edge);
00944 s_maxedge. push (maxedge);
00945
00946
00947 vertex = next;
00948
00949
00950 tab [vertex] = color;
00951
00952
00953 edge = vertex ? edgeEnds [vertex - 1] : 0;
00954 maxedge = edgeEnds [vertex];
00955 }
00956 }
00957 }
00958
00959 template <class wType> template <class Table, class Color>
00960 inline void diGraph<wType>::DFScolor (Table &tab, const Color &color,
00961 int vertex) const
00962 {
00963 DFScolorStack (tab, color, vertex);
00964 return;
00965 }
00966
00967
00968
00969 template <class wType> template <class Table>
00970 inline void diGraph<wType>::DFSfinishTimeRecurrent (Table &tab, int vertex,
00971 int &counter) const
00972 {
00973
00974 tab [vertex] = -1;
00975
00976
00977 for (int edge = vertex ? edgeEnds [vertex - 1] : 0;
00978 edge < edgeEnds [vertex]; ++ edge)
00979 {
00980 int next = edges [edge];
00981 if (!tab [next])
00982 DFSfinishTimeRecurrent (tab, next, counter);
00983 }
00984
00985
00986 tab [vertex] = ++ counter;
00987 return;
00988 }
00989
00990 template <class wType> template <class Table>
00991 inline void diGraph<wType>::DFSfinishTimeStack (Table &tab, int vertex,
00992 int &counter) const
00993 {
00994
00995 std::stack<int> s_vertex;
00996 std::stack<int> s_edge;
00997 std::stack<int> s_maxedge;
00998
00999
01000 tab [vertex] = -1;
01001
01002
01003 int edge = vertex ? edgeEnds [vertex - 1] : 0;
01004 int maxedge = edgeEnds [vertex];
01005
01006 while (1)
01007 {
01008
01009
01010 if (edge >= maxedge)
01011 {
01012
01013
01014 tab [vertex] = ++ counter;
01015
01016
01017 if (s_vertex. empty ())
01018 return;
01019
01020
01021 vertex = s_vertex. top ();
01022 s_vertex. pop ();
01023 edge = s_edge. top ();
01024 s_edge. pop ();
01025 maxedge = s_maxedge. top ();
01026 s_maxedge. pop ();
01027 continue;
01028 }
01029
01030
01031 int next = edges [edge ++];
01032 if (!tab [next])
01033 {
01034
01035 s_vertex. push (vertex);
01036 s_edge. push (edge);
01037 s_maxedge. push (maxedge);
01038
01039
01040 vertex = next;
01041
01042
01043 tab [vertex] = -1;
01044
01045
01046 edge = vertex ? edgeEnds [vertex - 1] : 0;
01047 maxedge = edgeEnds [vertex];
01048 }
01049 }
01050
01051 return;
01052 }
01053
01054 template <class wType> template <class Table>
01055 inline void diGraph<wType>::DFSfinishTime (Table &tab) const
01056 {
01057
01058 for (int i = 0; i < nVertices; ++ i)
01059 tab [i] = 0;
01060 int counter = 0;
01061
01062
01063 for (int i = 0; i < nVertices; ++ i)
01064 if (!tab [i])
01065 DFSfinishTimeStack (tab, i, counter);
01066
01067 #ifdef DIGRAPH_DEBUG
01068 int *tabdebug = new int [nVertices];
01069 for (int i = 0; i < nVertices; ++ i)
01070 tabdebug [i] = 0;
01071 int counterdebug = 0;
01072 for (int i = 0; i < nVertices; ++ i)
01073 if (!tabdebug [i])
01074 DFSfinishTimeRecurrent (tabdebug, i, counterdebug);
01075 if (counter != counterdebug)
01076 throw "DFSfinishTime: Wrong counter.";
01077 for (int i = 0; i < nVertices; ++ i)
01078 if (tab [i] != tabdebug [i])
01079 {
01080 sbug << "\nDFSfinishTime error: tabRec [" << i <<
01081 "] = " << tab [i] << ", tabStack [" << i <<
01082 "] = " << tabdebug [i] << ".\n";
01083 throw "DFSfinishTime: Wrong numbers.";
01084 }
01085 sbug << "DEBUG: DFSfinishTime OK. ";
01086 #endif // DIGRAPH_DEBUG
01087 return;
01088 }
01089
01090
01091
01092 template <class wType> template <class Table1, class Table2>
01093 inline bool diGraph<wType>::DFSforestRecurrent (Table1 &tab, Table1 &ntab,
01094 int vertex, int treeNumber, int countTrees,
01095 Table2 &compVertices, int &curVertex,
01096 diGraph<wType> *sccGraph, int *sccEdgeAdded) const
01097 {
01098
01099 compVertices [curVertex ++] = vertex;
01100
01101
01102 tab [vertex] = treeNumber;
01103
01104
01105
01106
01107 bool loop = false;
01108 for (int edge = vertex ? edgeEnds [vertex - 1] : 0;
01109 edge < edgeEnds [vertex]; ++ edge)
01110 {
01111 int next = edges [edge];
01112 if (!tab [next])
01113 loop |= DFSforestRecurrent (tab, ntab, next,
01114 treeNumber, countTrees, compVertices,
01115 curVertex, sccGraph, sccEdgeAdded);
01116 else if (tab [next] == treeNumber)
01117 {
01118 if (sccGraph)
01119 {
01120 int target = ntab [treeNumber - 1];
01121 if (sccEdgeAdded [target] != treeNumber)
01122 {
01123 sccGraph -> addEdge (target);
01124 sccEdgeAdded [target] = treeNumber;
01125 }
01126 }
01127 loop = true;
01128 }
01129 else if (sccGraph)
01130 {
01131 int target = ntab [tab [next] - 1];
01132 if ((target >= 0) &&
01133 (sccEdgeAdded [target] != treeNumber))
01134 {
01135 sccGraph -> addEdge (target);
01136 sccEdgeAdded [target] = treeNumber;
01137 }
01138 }
01139 }
01140
01141 return loop;
01142 }
01143
01144 template <class wType> template <class Table1, class Table2>
01145 inline bool diGraph<wType>::DFSforestStack (Table1 &tab, Table1 &ntab,
01146 int vertex, int treeNumber, int countTrees,
01147 Table2 &compVertices, int &curVertex,
01148 diGraph<wType> *sccGraph, int *sccEdgeAdded) const
01149 {
01150
01151 std::stack<int> s_vertex;
01152 std::stack<int> s_edge;
01153 std::stack<int> s_maxedge;
01154 std::stack<bool> s_loop;
01155
01156
01157 compVertices [curVertex ++] = vertex;
01158
01159
01160 tab [vertex] = treeNumber;
01161
01162
01163
01164
01165 bool loop = false;
01166 int edge = vertex ? edgeEnds [vertex - 1] : 0;
01167 int maxedge = edgeEnds [vertex];
01168 while (1)
01169 {
01170
01171
01172 if (edge >= maxedge)
01173 {
01174
01175 if (s_vertex. empty ())
01176 return loop;
01177
01178
01179 vertex = s_vertex. top ();
01180 s_vertex. pop ();
01181 edge = s_edge. top ();
01182 s_edge. pop ();
01183 maxedge = s_maxedge. top ();
01184 s_maxedge. pop ();
01185 loop |= s_loop. top ();
01186 s_loop. pop ();
01187 continue;
01188 }
01189
01190
01191 int next = edges [edge ++];
01192 if (!tab [next])
01193 {
01194
01195 s_vertex. push (vertex);
01196 s_edge. push (edge);
01197 s_maxedge. push (maxedge);
01198 s_loop. push (loop);
01199
01200
01201 vertex = next;
01202
01203
01204 compVertices [curVertex ++] = vertex;
01205
01206
01207 tab [vertex] = treeNumber;
01208
01209
01210 loop = false;
01211 edge = vertex ? edgeEnds [vertex - 1] : 0;
01212 maxedge = edgeEnds [vertex];
01213 }
01214 else if (tab [next] == treeNumber)
01215 {
01216 if (sccGraph)
01217 {
01218 int target = ntab [treeNumber - 1];
01219 if (sccEdgeAdded [target] != treeNumber)
01220 {
01221 sccGraph -> addEdge (target);
01222 sccEdgeAdded [target] = treeNumber;
01223 }
01224 }
01225 loop = true;
01226 }
01227 else if (sccGraph)
01228 {
01229 int target = ntab [tab [next] - 1];
01230 if ((target >= 0) &&
01231 (sccEdgeAdded [target] != treeNumber))
01232 {
01233 sccGraph -> addEdge (target);
01234 sccEdgeAdded [target] = treeNumber;
01235 }
01236 }
01237 }
01238
01239 return loop;
01240 }
01241
01242 template <class wType> template <class Table1, class Table2, class Table3>
01243 inline int diGraph<wType>::DFSforest (const Table1 &ordered,
01244 Table2 &compVertices, Table3 &compEnds, bool nontrivial,
01245 diGraph<wType> *sccGraph) const
01246 {
01247
01248
01249 int *tab = new int [nVertices];
01250 for (int i = 0; i < nVertices; ++ i)
01251 tab [i] = 0;
01252
01253
01254
01255 int *ntab = new int [nVertices];
01256
01257
01258
01259
01260 int *sccEdgeAdded = sccGraph ? new int [nVertices] : 0;
01261 if (sccGraph)
01262 {
01263 for (int n = 0; n < nVertices; ++ n)
01264 sccEdgeAdded [n] = -1;
01265 }
01266
01267
01268 int treeNumber = 0;
01269
01270
01271 int countTrees = 0;
01272 int curVertex = 0;
01273
01274
01275 for (int i = 0; i < nVertices; ++ i)
01276 {
01277
01278 int vertex = ordered [i];
01279
01280
01281 if (tab [vertex])
01282 continue;
01283
01284
01285 if (sccGraph)
01286 sccGraph -> addVertex ();
01287
01288
01289 int prevVertex = curVertex;
01290
01291
01292 if (sccGraph)
01293 ntab [treeNumber] = countTrees;
01294 ++ treeNumber;
01295 bool loop = DFSforestStack (tab, ntab, vertex,
01296 treeNumber, countTrees, compVertices,
01297 curVertex, sccGraph, sccEdgeAdded);
01298
01299
01300 compEnds [countTrees ++] = curVertex;
01301
01302
01303 if (nontrivial && !loop)
01304 {
01305 -- countTrees;
01306 curVertex = prevVertex;
01307 if (sccGraph)
01308 {
01309 ntab [treeNumber - 1] = -1;
01310 sccGraph -> removeVertex ();
01311 }
01312 }
01313 }
01314
01315 #ifdef DIGRAPH_DEBUG
01316 diGraph<wType> *sccGraphdebug = 0;
01317 if (sccGraph)
01318 sccGraphdebug = new diGraph<wType>;
01319
01320
01321 int *tabdebug = new int [nVertices];
01322 for (int i = 0; i < nVertices; ++ i)
01323 tabdebug [i] = 0;
01324
01325
01326
01327 int *ntabdebug = new int [nVertices];
01328
01329
01330
01331 int *sccEdgeAddeddebug = sccGraph ? new int [nVertices] : 0;
01332 if (sccGraph)
01333 {
01334 for (int n = 0; n < nVertices; ++ n)
01335 sccEdgeAddeddebug [n] = -1;
01336 }
01337
01338 int treeNumberdebug = 0;
01339
01340
01341 int countTreesdebug = 0;
01342 int curVertexdebug = 0;
01343
01344 int *compVerticesdebug = new int [nVertices];
01345 int *compEndsdebug = new int [nVertices];
01346
01347
01348 for (int i = 0; i < nVertices; ++ i)
01349 {
01350
01351 int vertex = ordered [i];
01352
01353
01354 if (tabdebug [vertex])
01355 continue;
01356
01357
01358 if (sccGraphdebug)
01359 sccGraphdebug -> addVertex ();
01360
01361
01362 int prevVertex = curVertexdebug;
01363
01364
01365 if (sccGraphdebug)
01366 ntabdebug [treeNumberdebug] = countTreesdebug;
01367 ++ treeNumberdebug;
01368 bool loop = DFSforestRecurrent (tabdebug, ntabdebug, vertex,
01369 treeNumberdebug, countTreesdebug, compVerticesdebug,
01370 curVertexdebug, sccGraphdebug, sccEdgeAddeddebug);
01371
01372
01373 compEndsdebug [countTreesdebug ++] = curVertexdebug;
01374
01375
01376 if (nontrivial && !loop)
01377 {
01378 -- countTreesdebug;
01379 curVertexdebug = prevVertex;
01380 if (sccGraphdebug)
01381 {
01382 ntabdebug [treeNumberdebug - 1] = -1;
01383 sccGraphdebug -> removeVertex ();
01384 }
01385 }
01386 }
01387 if (countTrees != countTreesdebug)
01388 throw "DFSforest: Wrong countTrees.";
01389 for (int i = 0; i < countTrees; ++ i)
01390 if (compEnds [i] != compEndsdebug [i])
01391 throw "DFSforest: Wrong compEnds.";
01392 for (int i = 0; i < compEndsdebug [countTrees - 1]; ++ i)
01393 if (compVertices [i] != compVerticesdebug [i])
01394 throw "DFSforest: Wrong vertices.";
01395 if (curVertex != curVertexdebug)
01396 throw "DFSforest: Wrong curVertex.";
01397 for (int i = 0; i < nVertices; ++ i)
01398 if (tab [i] != tabdebug [i])
01399 throw "DFSforest: Wrong tab.";
01400 if (sccGraph)
01401 {
01402 for (int i = 0; i < nVertices; ++ i)
01403 if (ntab [i] != ntabdebug [i])
01404 throw "DFSforest: Wrong ntab.";
01405 if (*sccGraph != *sccGraphdebug)
01406 throw "DFSforest: Wrong graph.";
01407 }
01408 if (sccEdgeAdded)
01409 {
01410 for (int i = 0; i < nVertices; ++ i)
01411 if (sccEdgeAdded [i] != sccEdgeAddeddebug [i])
01412 throw "DFSforest: Wrong sccEdgeAdded.";
01413 }
01414 if (treeNumber != treeNumberdebug)
01415 throw "DFSforest: Wrong treeNumber.";
01416 sbug << "DEBUG: DFSforest OK. ";
01417 if (!sccGraph)
01418 sbug << "(Graphs not compared.) ";
01419 delete [] compVerticesdebug;
01420 delete [] compEndsdebug;
01421 if (sccGraphdebug)
01422 delete sccGraphdebug;
01423 delete [] ntabdebug;
01424 delete [] tabdebug;
01425 if (sccEdgeAddeddebug)
01426 delete [] sccEdgeAddeddebug;
01427 #endif // DIGRAPH_DEBUG
01428
01429 if (sccEdgeAdded)
01430 delete [] sccEdgeAdded;
01431 delete [] ntab;
01432 delete [] tab;
01433 return countTrees;
01434 }
01435
01436 template <class wType>
01437 inline int diGraph<wType>::shortestPath (int source, int destination) const
01438 {
01439
01440 if ((source < 0) || (source >= nVertices) ||
01441 (destination < 0) || (destination >= nVertices))
01442 {
01443 throw "diGraph - shortest path: Wrong vertex number.";
01444 }
01445
01446 if (source >= 0)
01447 throw "Shortest path algorithm not implemented, yet.";
01448
01449
01450 std::auto_ptr<int> lengths_ptr (new int [nVertices]);
01451 int *lengths = lengths_ptr. get ();
01452 for (int i = 0; i < nVertices; ++ i)
01453 lengths [i] = -1;
01454
01455
01456 std::stack<int> s_vertex;
01457 std::stack<int> s_edge;
01458 std::stack<int> s_maxedge;
01459
01460
01461 int vertex = source;
01462 lengths [vertex] = 0;
01463
01464
01465 int edge = vertex ? edgeEnds [vertex - 1] : 0;
01466 int maxedge = edgeEnds [vertex];
01467
01468 while (1)
01469 {
01470
01471
01472 if (edge >= maxedge)
01473 {
01474
01475 if (s_vertex. empty ())
01476 return lengths [destination];
01477
01478
01479 vertex = s_vertex. top ();
01480 s_vertex. pop ();
01481 edge = s_edge. top ();
01482 s_edge. pop ();
01483 maxedge = s_maxedge. top ();
01484 s_maxedge. pop ();
01485 continue;
01486 }
01487
01488
01489 int next = edges [edge ++];
01490
01491
01492 if (lengths [next] != -1)
01493 {
01494 if ((lengths [next] > lengths [vertex] + 1) ||
01495 !lengths [next])
01496 {
01497 lengths [next] = lengths [vertex] + 1;
01498 }
01499 }
01500
01501
01502 else
01503 {
01504
01505 s_vertex. push (vertex);
01506 s_edge. push (edge);
01507 s_maxedge. push (maxedge);
01508
01509
01510 lengths [next] = lengths [vertex] + 1;
01511
01512
01513 vertex = next;
01514
01515
01516 edge = vertex ? edgeEnds [vertex - 1] : 0;
01517 maxedge = edgeEnds [vertex];
01518 }
01519 }
01520
01521 return 0;
01522 }
01523
01524 template <class wType>
01525 inline int diGraph<wType>::shortestLoop (int origin) const
01526 {
01527 return shortestPath (origin, origin);
01528 }
01529
01530 template <class wType>
01531 template <class lenTable, class weightsType, class roundType>
01532 inline void diGraph<wType>::Dijkstra (const roundType &rounding,
01533 int source, lenTable &len, weightsType &edgeWeights) const
01534 {
01535
01536 FibonacciHeap<posWeight> Q (nVertices);
01537
01538
01539
01540 for (int v = 0; v < nVertices; ++ v)
01541 {
01542 posWeight w (0);
01543 if (v != source)
01544 w. setInfinity ();
01545 int number = Q. Insert (w);
01546 if (number != v)
01547 {
01548 throw "Wrong implementation of Fibonacci heap "
01549 "for this version of Dijkstra.";
01550 }
01551 }
01552
01553
01554
01555 for (int i = 0; i < nVertices; ++ i)
01556 {
01557
01558 int minVertex = Q. Minimum ();
01559 posWeight minWeight = Q. ExtractMinimum ();
01560
01561 if (minWeight. isInfinity ())
01562 {
01563 len [minVertex] = -1;
01564 continue;
01565 }
01566 wType minValue = minWeight. getValue ();
01567 len [minVertex] = minValue;
01568
01569
01570
01571 int edge = minVertex ? edgeEnds [minVertex - 1] : 0;
01572 int maxEdge = edgeEnds [minVertex];
01573 for (; edge < maxEdge; ++ edge)
01574 {
01575
01576 int nextVertex = edges [edge];
01577
01578
01579
01580 const posWeight &nextWeight = Q. Value (nextVertex);
01581 wType newWeight = rounding. add_down (minValue,
01582 edgeWeights [edge]);
01583 if (newWeight < 0)
01584 newWeight = 0;
01585 if (nextWeight. isInfinity () ||
01586 (newWeight < nextWeight. getValue ()))
01587 {
01588 Q. DecreaseKey (nextVertex,
01589 posWeight (newWeight));
01590 }
01591 }
01592 }
01593 return;
01594 }
01595
01596 template <class wType>
01597 template <class lenTable, class roundType>
01598 inline void diGraph<wType>::Dijkstra (const roundType &rounding,
01599 int source, lenTable &len) const
01600 {
01601 this -> Dijkstra (rounding, source, len, this -> weights);
01602 return;
01603 }
01604
01605 template <class wType>
01606 template <class lenTable>
01607 inline void diGraph<wType>::Dijkstra (int source, lenTable &len) const
01608 {
01609 const dummyRounding<wType> rounding = dummyRounding<wType> ();
01610 this -> Dijkstra (rounding, source, len);
01611 return;
01612 }
01613
01614 template <class wType> template <class lenTable, class predTable,
01615 class roundType>
01616 inline bool diGraph<wType>::BellmanFord (const roundType &rounding,
01617 int source, lenTable &len, wType *infinity, predTable pred) const
01618 {
01619
01620 if ((source < 0) || (source >= nVertices))
01621 throw "Bellman-Ford: Wrong source vertex number.";
01622
01623
01624 BitField finite;
01625 finite. allocate (nVertices);
01626 finite. clearall (nVertices);
01627 finite. set (source);
01628 len [source] = 0;
01629
01630
01631 int countNegative = 0;
01632
01633
01634 for (int i = 0; i < nVertices; ++ i)
01635 pred [i] = -1;
01636
01637
01638 bool noNegativeLoop = false;
01639 int counter = 0;
01640 for (; counter <= nVertices; ++ counter)
01641 {
01642 bool modified = false;
01643 int curEdge = 0;
01644 for (int vertex = 0; vertex < nVertices; ++ vertex)
01645 {
01646 int maxEdge = edgeEnds [vertex];
01647 if (!finite. test (vertex))
01648 {
01649 curEdge = maxEdge;
01650 continue;
01651 }
01652 for (; curEdge < maxEdge; ++ curEdge)
01653 {
01654 int next = edges [curEdge];
01655 wType newlen = rounding. add_down
01656 (len [vertex], weights [curEdge]);
01657 if (!finite. test (next))
01658 {
01659 finite. set (next);
01660 modified = true;
01661 len [next] = newlen;
01662 pred [next] = vertex;
01663 if (newlen < 0)
01664 ++ countNegative;
01665 }
01666 else if (newlen < len [next])
01667 {
01668 modified = true;
01669 if (!(len [next] < 0) &&
01670 (newlen < 0))
01671 {
01672 ++ countNegative;
01673 }
01674 len [next] = newlen;
01675 pred [next] = vertex;
01676 }
01677 }
01678 }
01679 if (countNegative == nVertices)
01680 {
01681 noNegativeLoop = false;
01682 ++ counter;
01683 break;
01684 }
01685 if (!modified)
01686 {
01687 noNegativeLoop = true;
01688 ++ counter;
01689 break;
01690 }
01691 }
01692
01693
01694 if (false && chomp::homology::sbug. show)
01695 {
01696 chomp::homology::sbug << "Bellman-Ford: " <<
01697 counter << ((counter > 1) ? " loops (" : " loop (") <<
01698 nVertices << " vertices, " << countNegative <<
01699 " negative). " <<
01700 (noNegativeLoop ? "No negative loops.\n" :
01701 "A negative loop found.\n");
01702 }
01703
01704
01705 if (infinity && noNegativeLoop)
01706 {
01707 wType infty (0);
01708 bool first = true;
01709 for (int i = 0; i < nVertices; ++ i)
01710 {
01711 if (!finite. test (i))
01712 continue;
01713 if (first)
01714 {
01715 infty = len [i];
01716 first = false;
01717 }
01718 else if (infty < len [i])
01719 {
01720 infty = len [i];
01721 }
01722 }
01723 infty = infty + 1;
01724 for (int i = 0; i < nVertices; ++ i)
01725 {
01726 if (!finite. test (i))
01727 len [i] = infty;
01728 }
01729 *infinity = infty;
01730 }
01731
01732 finite. free ();
01733 return noNegativeLoop;
01734 }
01735
01736 template <class wType> template <class lenTable, class predTable>
01737 inline bool diGraph<wType>::BellmanFord (int source, lenTable &len,
01738 wType *infinity, predTable pred) const
01739 {
01740 const dummyRounding<wType> rounding = dummyRounding<wType> ();
01741 return this -> BellmanFord (rounding, source, len, infinity, pred);
01742 }
01743
01744 template <class wType> template <class roundType>
01745 inline bool diGraph<wType>::BellmanFord (const roundType &rounding,
01746 int source) const
01747 {
01748 std::auto_ptr<wType> len_ptr (new wType [nVertices]);
01749 wType *len = len_ptr. get ();
01750 wType *infinity = 0;
01751 dummyArray tab;
01752 return BellmanFord (rounding, source, len, infinity, tab);
01753 }
01754
01755 template <class wType>
01756 inline bool diGraph<wType>::BellmanFord (int source) const
01757 {
01758 const dummyRounding<wType> rounding = dummyRounding<wType> ();
01759 return BellmanFord (rounding, source);
01760 }
01761
01762 template <class wType>
01763 inline wType diGraph<wType>::Edmonds () const
01764 {
01765
01766 std::vector<edgeTriple> edgeTriples (countEdges ());
01767 int prevEdge = 0;
01768 int curEdge = 0;
01769 for (int vert = 0; vert < nVertices; ++ vert)
01770 {
01771 while (curEdge < edgeEnds [vert])
01772 {
01773 edgeTriple &e = edgeTriples [curEdge];
01774 e. vert1 = vert;
01775 e. vert2 = edges [curEdge];
01776 e. weight = weights [curEdge];
01777 ++ curEdge;
01778 }
01779 prevEdge = curEdge;
01780 }
01781 std::sort (edgeTriples. begin (), edgeTriples. end ());
01782
01783
01784 std::auto_ptr<int> root_ptr (new int [nVertices]);
01785 int *root = root_ptr. get ();
01786 for (int vert = 0; vert < nVertices; ++ vert)
01787 {
01788 root [vert] = -1;
01789 }
01790
01791
01792 wType totalWeight = 0;
01793 int nEdges = countEdges ();
01794 for (int curEdge = 0; curEdge < nEdges; ++ curEdge)
01795 {
01796
01797
01798 edgeTriple &e = edgeTriples [curEdge];
01799 int root1 = e. vert1;
01800 while (root [root1] >= 0)
01801 {
01802 root1 = root [root1];
01803 }
01804 int vert1 = e. vert1;
01805 while (root [vert1] >= 0)
01806 {
01807 int next = root [vert1];
01808 root [vert1] = root1;
01809 vert1 = next;
01810 }
01811 int root2 = e. vert2;
01812 while (root [root2] >= 0)
01813 {
01814 root2 = root [root2];
01815 }
01816 int vert2 = e. vert2;
01817 while (root [vert2] >= 0)
01818 {
01819 int next = root [vert2];
01820 root [vert2] = root2;
01821 vert2 = next;
01822 }
01823
01824
01825 if (root1 == root2)
01826 continue;
01827
01828
01829 totalWeight += e. weight;
01830
01831
01832 root [root1] = root2;
01833 }
01834 return totalWeight;
01835 }
01836
01837 template <class wType>
01838 inline wType diGraph<wType>::EdmondsOld () const
01839 {
01840
01841 std::vector<edgeTriple> edgeTriples (countEdges ());
01842 int prevEdge = 0;
01843 int curEdge = 0;
01844 for (int vert = 0; vert < nVertices; ++ vert)
01845 {
01846 while (curEdge < edgeEnds [vert])
01847 {
01848 edgeTriple &e = edgeTriples [curEdge];
01849 e. vert1 = vert;
01850 e. vert2 = edges [curEdge];
01851 e. weight = weights [curEdge];
01852 ++ curEdge;
01853 }
01854 prevEdge = curEdge;
01855 }
01856 std::sort (edgeTriples. begin (), edgeTriples. end ());
01857
01858
01859 std::auto_ptr<int> forest_ptr (new int [nVertices]);
01860 int *forest = forest_ptr. get ();
01861 std::auto_ptr<int> next_ptr (new int [nVertices]);
01862 int *next = next_ptr. get ();
01863 std::auto_ptr<int> prev_ptr (new int [nVertices]);
01864 int *prev = prev_ptr. get ();
01865 for (int vert = 0; vert < nVertices; ++ vert)
01866 {
01867 forest [vert] = vert;
01868 next [vert] = -1;
01869 prev [vert] = -1;
01870 }
01871
01872
01873 wType totalWeight = 0;
01874 int nEdges = countEdges ();
01875 for (int curEdge = 0; curEdge < nEdges; ++ curEdge)
01876 {
01877
01878 edgeTriple &e = edgeTriples [curEdge];
01879 if (forest [e. vert1] == forest [e. vert2])
01880 continue;
01881
01882
01883 totalWeight += e. weight;
01884
01885
01886 int vert1 = e. vert1;
01887 while (next [vert1] >= 0)
01888 {
01889 vert1 = next [vert1];
01890 }
01891
01892
01893 int vert2 = e. vert2;
01894 while (prev [vert2] >= 0)
01895 {
01896 vert2 = prev [vert2];
01897 }
01898
01899
01900 next [vert1] = vert2;
01901 prev [vert2] = vert1;
01902 int treeNumber = forest [vert1];
01903 while (vert2 >= 0)
01904 {
01905 forest [vert2] = treeNumber;
01906 vert2 = next [vert2];
01907 }
01908 }
01909 return totalWeight;
01910 }
01911
01912 template <class wType>
01913 template <class arrayType, class roundType>
01914 inline wType diGraph<wType>::FloydWarshall (const roundType &rounding,
01915 arrayType &arr, bool setInfinity, bool ignoreNegLoop) const
01916 {
01917
01918 if (!nVertices)
01919 return 0;
01920
01921
01922 BitField *finite = new BitField [nVertices];
01923 for (int i = 0; i < nVertices; ++ i)
01924 {
01925 finite [i]. allocate (nVertices);
01926 finite [i]. clearall (nVertices);
01927 }
01928
01929
01930 int curEdge = 0;
01931 for (int source = 0; source < nVertices; ++ source)
01932 {
01933 bool diagset = false;
01934 while (curEdge < edgeEnds [source])
01935 {
01936 int target = edges [curEdge];
01937 const wType &weight = weights [curEdge];
01938 if (source == target)
01939 {
01940 if (weight < 0)
01941 {
01942 arr [source] [target] = weight;
01943 diagset = true;
01944 }
01945 }
01946 else
01947 {
01948 arr [source] [target] = weight;
01949 finite [source]. set (target);
01950 }
01951 ++ curEdge;
01952 }
01953 if (!diagset)
01954 arr [source] [source] = 0;
01955 finite [source]. set (source);
01956 }
01957
01958
01959 for (int k = 0; k < nVertices; ++ k)
01960 {
01961 for (int i = 0; i < nVertices; ++ i)
01962 {
01963 if (!finite [i]. test (k))
01964 continue;
01965 for (int j = 0; j < nVertices; ++ j)
01966 {
01967 if (!finite [k]. test (j))
01968 continue;
01969 const wType sum = rounding. add_down
01970 (arr [i] [k], arr [k] [j]);
01971 if (finite [i]. test (j))
01972 {
01973 if (sum < arr [i] [j])
01974 arr [i] [j] = sum;
01975 }
01976 else
01977 {
01978 arr [i] [j] = sum;
01979 finite [i]. set (j);
01980 }
01981 }
01982 }
01983 }
01984
01985
01986
01987 if (!ignoreNegLoop)
01988 {
01989 for (int i = 0; i < nVertices; ++ i)
01990 {
01991 if (arr [i] [i] < 0)
01992 throw "Negative loop in Floyd-Warshall.";
01993 }
01994 }
01995
01996
01997 wType result = 0;
01998
01999
02000
02001 if (setInfinity)
02002 {
02003 wType &infinity = result;
02004 for (int i = 0; i < nVertices; ++ i)
02005 {
02006 for (int j = 0; j < nVertices; ++ j)
02007 {
02008 if (finite [i]. test (j) &&
02009 (infinity <= arr [i] [j]))
02010 {
02011 infinity = rounding. add_up
02012 (arr [i] [j], 1);
02013 }
02014 }
02015 }
02016 for (int i = 0; i < nVertices; ++ i)
02017 {
02018 for (int j = 0; j < nVertices; ++ j)
02019 {
02020 if (!finite [i]. test (j))
02021 arr [i] [j] = infinity;
02022 }
02023 }
02024 }
02025
02026
02027 else
02028 {
02029 wType &minWeight = result;
02030 bool firstTime = true;
02031 for (int i = 0; i < nVertices; ++ i)
02032 {
02033 for (int j = 0; j < nVertices; ++ j)
02034 {
02035 if (finite [i]. test (j))
02036 {
02037 if (firstTime)
02038 {
02039 minWeight = arr [i] [j];
02040 firstTime = false;
02041 }
02042 else if (arr [i] [j] < minWeight)
02043 {
02044 minWeight = arr [i] [j];
02045 }
02046 }
02047 }
02048 }
02049 }
02050
02051
02052 for (int i = 0; i < nVertices; ++ i)
02053 finite [i]. free ();
02054 delete [] finite;
02055
02056 return result;
02057 }
02058
02059 template <class wType>
02060 template <class arrayType>
02061 inline wType diGraph<wType>::FloydWarshall (arrayType &arr,
02062 bool setInfinity, bool ignoreNegLoop) const
02063 {
02064 const dummyRounding<wType> rounding = dummyRounding<wType> ();
02065 return FloydWarshall (rounding, arr, setInfinity, ignoreNegLoop);
02066 }
02067
02068 template <class wType>
02069 template <class arrayType, class roundType>
02070 inline wType diGraph<wType>::Johnson (const roundType &rounding,
02071 arrayType &arr, bool setInfinity, bool ignoreNegLoop) const
02072 {
02073
02074 if (false && sbug. show)
02075 {
02076 timeused stopWatch;
02077 wType res = FloydWarshall (rounding,
02078 arr, setInfinity, ignoreNegLoop);
02079 chomp::homology::sbug << "\n[Floyd-Warshall: " << res <<
02080 ", " << (double) stopWatch << " sec]\n";
02081 }
02082
02083
02084
02085
02086 if (!nVertices)
02087 return 0;
02088
02089
02090
02091
02092 wType *len = new wType [nVertices];
02093 for (int i = 0; i < nVertices; ++ i)
02094 len [i] = 0;
02095
02096
02097 bool noNegativeLoop = false;
02098 int counter = 0;
02099 for (; counter <= nVertices; ++ counter)
02100 {
02101 bool modified = false;
02102 int curEdge = 0;
02103 for (int vertex = 0; vertex < nVertices; ++ vertex)
02104 {
02105 int maxEdge = edgeEnds [vertex];
02106 for (; curEdge < maxEdge; ++ curEdge)
02107 {
02108 int next = edges [curEdge];
02109 wType newlen = rounding. add_down
02110 (len [vertex], weights [curEdge]);
02111 if (newlen < len [next])
02112 {
02113
02114
02115 if (counter > nVertices)
02116 {
02117 std::cout << vertex;
02118 }
02119 modified = true;
02120 len [next] = newlen;
02121 }
02122 }
02123 }
02124 if (!modified)
02125 {
02126 noNegativeLoop = true;
02127 ++ counter;
02128 break;
02129 }
02130 }
02131 if (!ignoreNegLoop && !noNegativeLoop)
02132 throw "Negative loop found in Johnson's algorithm.";
02133
02134
02135 wType *newWeights = new wType [edgeEnds [nVertices - 1]];
02136 int edge = 0;
02137 for (int vertex = 0; vertex < nVertices; ++ vertex)
02138 {
02139 int maxEdge = edgeEnds [vertex];
02140 for (; edge < maxEdge; ++ edge)
02141 {
02142 wType w = weights [edge];
02143 w = rounding. add_down (w, len [vertex]);
02144 w = rounding. sub_down (w, len [edges [edge]]);
02145 newWeights [edge] = (w < 0) ? 0 : w;
02146 }
02147 }
02148
02149
02150
02151
02152 for (int u = 0; u < nVertices; ++ u)
02153 {
02154 this -> Dijkstra (rounding, u, arr [u], newWeights);
02155 }
02156 delete [] newWeights;
02157
02158
02159
02160
02161 wType result = 0;
02162 if (setInfinity)
02163 {
02164 wType &infinity = result;
02165 for (int u = 0; u < nVertices; ++ u)
02166 {
02167 for (int v = 0; v < nVertices; ++ v)
02168 {
02169 wType w = arr [u] [v];
02170 if (w < 0)
02171 continue;
02172 w = rounding. add_down (w, len [v]);
02173 w = rounding. sub_down (w, len [u]);
02174 if (w < infinity)
02175 continue;
02176 infinity = rounding. add_up (w, 1);
02177 }
02178 }
02179 for (int u = 0; u < nVertices; ++ u)
02180 {
02181 for (int v = 0; v < nVertices; ++ v)
02182 {
02183 wType w = arr [u] [v];
02184 if (w < 0)
02185 {
02186 arr [u] [v] = infinity;
02187 continue;
02188 }
02189 w = rounding. add_down (w, len [v]);
02190 arr [u] [v] =
02191 rounding. sub_down (w, len [u]);
02192 }
02193 }
02194 }
02195 else
02196 {
02197 wType &minWeight = result;
02198 bool firstTime = true;
02199 for (int u = 0; u < nVertices; ++ u)
02200 {
02201 for (int v = 0; v < nVertices; ++ v)
02202 {
02203 wType w = arr [u] [v];
02204 if (w < 0)
02205 continue;
02206 w = rounding. add_down (w, len [v]);
02207 w = rounding. sub_down (w, len [u]);
02208 if (firstTime)
02209 {
02210 minWeight = w;
02211 firstTime = false;
02212 }
02213 else if (w < minWeight)
02214 minWeight = w;
02215 }
02216 }
02217 }
02218 delete [] len;
02219
02220
02221 if (false && sbug. show)
02222 {
02223
02224
02225 }
02226
02227 return result;
02228 }
02229
02230 template <class wType>
02231 template <class arrayType>
02232 inline wType diGraph<wType>::Johnson (arrayType &arr,
02233 bool setInfinity, bool ignoreNegLoop) const
02234 {
02235 const dummyRounding<wType> rounding = dummyRounding<wType> ();
02236 return Johnson (rounding, arr, setInfinity, ignoreNegLoop);
02237 }
02238
02239 template <class wType>
02240 template <class roundType>
02241 wType diGraph<wType>::minPathWeight (const roundType &rounding,
02242 bool ignoreNegLoop, int sparseGraph) const
02243 {
02244
02245 if (nVertices <= 0)
02246 return 0;
02247
02248
02249 wType **arr = new wType * [nVertices];
02250 for (int i = 0; i < nVertices; ++ i)
02251 arr [i] = new wType [nVertices];
02252
02253
02254
02255 bool sparse = false;
02256 if (sparseGraph == 1)
02257 sparse = true;
02258 else if (sparseGraph != 0)
02259 {
02260 double nEdgesD = this -> countEdges ();
02261 double nVerticesD = this -> countVertices ();
02262 if ((nVerticesD > 3000) && (nEdgesD < nVerticesD * 1000) &&
02263 (nEdgesD < nVerticesD * nVerticesD / 50))
02264 {
02265 sparse = true;
02266 }
02267 }
02268
02269
02270 wType minWeight = sparse ?
02271 this -> Johnson (rounding, arr, false, ignoreNegLoop) :
02272 this -> FloydWarshall (rounding, arr, false, ignoreNegLoop);
02273
02274
02275
02276
02277
02278
02279
02280
02281
02282
02283
02284
02285
02286
02287 for (int i = 0; i < nVertices; ++ i)
02288 delete [] (arr [i]);
02289 delete [] arr;
02290
02291 return minWeight;
02292 }
02293
02294 template <class wType>
02295 wType diGraph<wType>::minPathWeight (bool ignoreNegLoop, int sparseGraph)
02296 const
02297 {
02298 const dummyRounding<wType> rounding = dummyRounding<wType> ();
02299 return this -> minPathWeight (rounding, ignoreNegLoop, sparseGraph);
02300 }
02301
02302 template <class wType> template <class outType>
02303 inline outType &diGraph<wType>::show (outType &out, bool showWeights) const
02304 {
02305 out << "; Directed graph: " << nVertices << " vertices.\n";
02306 int curEdge = 0;
02307 for (int i = 0; i < nVertices; ++ i)
02308 {
02309 for (; curEdge < edgeEnds [i]; ++ curEdge)
02310 {
02311 out << i << " -> " << edges [curEdge];
02312 if (showWeights)
02313 out << " [" << weights [curEdge] << "]\n";
02314 else
02315 out << "\n";
02316 }
02317 }
02318 out << "; This is the end of the graph.\n";
02319 return out;
02320 }
02321
02322
02323
02324
02325
02326 template <class wType>
02327 inline std::ostream &operator << (std::ostream &out, const diGraph<wType> &g)
02328 {
02329 return g. show (out, false);
02330 }
02331
02332
02333
02334
02335
02336
02337
02338
02339
02340
02341 template <class wType, class Table1, class Table2>
02342 inline int SCC (const diGraph<wType> &g, Table1 &compVertices,
02343 Table2 &compEnds, diGraph<wType> *scc = 0,
02344 diGraph<wType> *transposed = 0, bool copyweights = false)
02345 {
02346
02347 int nVert = g. countVertices ();
02348 int *ordered = new int [nVert];
02349 int *tab = new int [nVert];
02350
02351
02352 g. DFSfinishTime (tab);
02353 for (int i = 0; i < nVert; ++ i)
02354 ordered [nVert - tab [i]] = i;
02355 delete [] tab;
02356
02357
02358 diGraph<wType> gT;
02359 if (!transposed)
02360 transposed = &gT;
02361 g. transpose (*transposed, copyweights);
02362
02363
02364 int n = transposed -> DFSforest (ordered, compVertices, compEnds,
02365 true, scc);
02366
02367
02368 delete [] ordered;
02369 return n;
02370 }
02371
02372
02373
02374 template <class wType>
02375 wType diGraph<wType>::minMeanCycleWeight (diGraph<wType> *transposed) const
02376 {
02377
02378 multitable<int> compVertices, compEnds;
02379 bool copyweights = !!transposed;
02380 int countSCC = SCC (*this, compVertices, compEnds,
02381 static_cast<diGraph<wType> *> (0), transposed, copyweights);
02382 if (countSCC <= 0)
02383 return 0;
02384
02385
02386 int maxCompSize = compEnds [0];
02387 for (int comp = 1; comp < countSCC; ++ comp)
02388 {
02389 int compSize = compEnds [comp] - compEnds [comp - 1];
02390 if (maxCompSize < compSize)
02391 maxCompSize = compSize;
02392 }
02393
02394
02395 wType **F = new wType * [maxCompSize + 1];
02396 BitField *finite = new BitField [maxCompSize + 1];
02397 for (int i = 0; i <= maxCompSize; ++ i)
02398 {
02399 F [i] = new wType [maxCompSize];
02400 finite [i]. allocate (maxCompSize);
02401 }
02402
02403
02404 int *numbers = new int [nVertices];
02405 int *components = new int [nVertices];
02406 for (int i = 0; i < nVertices; ++ i)
02407 components [i] = -1;
02408 int offset = 0;
02409 for (int comp = 0; comp < countSCC; ++ comp)
02410 {
02411 int maxOffset = compEnds [comp];
02412 int pos = offset;
02413 for (; pos < maxOffset; ++ pos)
02414 {
02415 numbers [compVertices [pos]] = pos - offset;
02416 components [compVertices [pos]] = comp;
02417 }
02418 offset = pos;
02419 }
02420
02421
02422 wType minWeight (0);
02423 for (int comp = 0; comp < countSCC; ++ comp)
02424 {
02425 int offset = comp ? compEnds [comp - 1] : 0;
02426 int compSize = compEnds [comp] - offset;
02427 for (int i = 0; i <= compSize; ++ i)
02428 finite [i]. clearall (compSize);
02429 F [0] [0] = 0;
02430 finite [0]. set (0);
02431
02432 for (int len = 1; len <= compSize; ++ len)
02433 {
02434
02435 for (int i = 0; i < compSize; ++ i)
02436 {
02437 if (!finite [len - 1]. test (i))
02438 continue;
02439 wType prevWeight = F [len - 1] [i];
02440 int source = compVertices [offset + i];
02441
02442
02443 int edgeOffset = source ?
02444 edgeEnds [source - 1] : 0;
02445 int edgeEnd = edgeEnds [source];
02446 for (; edgeOffset < edgeEnd; ++ edgeOffset)
02447 {
02448 int target = edges [edgeOffset];
02449 if (components [target] != comp)
02450 continue;
02451 int j = numbers [target];
02452 wType newWeight = prevWeight +
02453 weights [edgeOffset];
02454 if (!finite [len]. test (j))
02455 {
02456 finite [len]. set (j);
02457 F [len] [j] = newWeight;
02458 }
02459 else if (newWeight < F [len] [j])
02460 F [len] [j] = newWeight;
02461 }
02462 }
02463 }
02464
02465
02466 wType minCompWeight (0);
02467 bool firstMin = true;
02468 for (int vert = 0; vert < compSize; ++ vert)
02469 {
02470 if (!finite [compSize]. test (vert))
02471 continue;
02472 bool firstAverage = true;
02473 wType maxAverage = 0;
02474 for (int first = 0; first < compSize; ++ first)
02475 {
02476 if (!finite [first]. test (vert))
02477 continue;
02478 wType average = (F [compSize] [vert] -
02479 F [first] [vert]) /
02480 (compSize - first);
02481 if (firstAverage)
02482 {
02483 maxAverage = average;
02484 firstAverage = false;
02485 }
02486 else if (maxAverage < average)
02487 maxAverage = average;
02488 }
02489 if (firstMin || (maxAverage < minCompWeight))
02490 {
02491 if (firstAverage)
02492 throw "Bug in Karp's algorithm";
02493 minCompWeight = maxAverage;
02494 firstMin = false;
02495 }
02496 }
02497
02498
02499 if (!comp || (minCompWeight < minWeight))
02500 minWeight = minCompWeight;
02501 }
02502
02503
02504 delete [] components;
02505 delete [] numbers;
02506 for (int i = 0; i < maxCompSize; ++ i)
02507 {
02508 finite [i]. free ();
02509 delete [] (F [i]);
02510 }
02511 delete [] finite;
02512 delete [] F;
02513
02514
02515 return minWeight;
02516 }
02517
02518 template <class wType>
02519 template <class roundType>
02520 wType diGraph<wType>::minMeanCycleWeight (const roundType &rounding,
02521 diGraph<wType> *transposed) const
02522 {
02523
02524 multitable<int> compVertices, compEnds;
02525 bool copyweights = !!transposed;
02526 int countSCC = SCC (*this, compVertices, compEnds,
02527 static_cast<diGraph<wType> *> (0), transposed, copyweights);
02528 if (countSCC <= 0)
02529 return 0;
02530
02531
02532 int maxCompSize = compEnds [0];
02533 for (int comp = 1; comp < countSCC; ++ comp)
02534 {
02535 int compSize = compEnds [comp] - compEnds [comp - 1];
02536 if (maxCompSize < compSize)
02537 maxCompSize = compSize;
02538 }
02539
02540
02541 wType **FUpper = new wType * [maxCompSize + 1];
02542 wType **FLower = new wType * [maxCompSize + 1];
02543 BitField *finite = new BitField [maxCompSize + 1];
02544 for (int i = 0; i <= maxCompSize; ++ i)
02545 {
02546 FUpper [i] = new wType [maxCompSize];
02547 FLower [i] = new wType [maxCompSize];
02548 finite [i]. allocate (maxCompSize);
02549 }
02550
02551
02552 int *numbers = new int [nVertices];
02553 int *components = new int [nVertices];
02554 for (int i = 0; i < nVertices; ++ i)
02555 components [i] = -1;
02556 int offset = 0;
02557 for (int comp = 0; comp < countSCC; ++ comp)
02558 {
02559 int maxOffset = compEnds [comp];
02560 int pos = offset;
02561 for (; pos < maxOffset; ++ pos)
02562 {
02563 numbers [compVertices [pos]] = pos - offset;
02564 components [compVertices [pos]] = comp;
02565 }
02566 offset = pos;
02567 }
02568
02569
02570 wType minWeight (0);
02571 for (int comp = 0; comp < countSCC; ++ comp)
02572 {
02573 int offset = comp ? compEnds [comp - 1] : 0;
02574 int compSize = compEnds [comp] - offset;
02575 for (int i = 0; i <= compSize; ++ i)
02576 finite [i]. clearall (compSize);
02577 FUpper [0] [0] = 0;
02578 FLower [0] [0] = 0;
02579 finite [0]. set (0);
02580
02581 for (int len = 1; len <= compSize; ++ len)
02582 {
02583
02584 for (int i = 0; i < compSize; ++ i)
02585 {
02586 if (!finite [len - 1]. test (i))
02587 continue;
02588 wType prevUpper = FUpper [len - 1] [i];
02589 wType prevLower = FLower [len - 1] [i];
02590 int source = compVertices [offset + i];
02591
02592
02593 int edgeOffset = source ?
02594 edgeEnds [source - 1] : 0;
02595 int edgeEnd = edgeEnds [source];
02596 for (; edgeOffset < edgeEnd; ++ edgeOffset)
02597 {
02598 int target = edges [edgeOffset];
02599 if (components [target] != comp)
02600 continue;
02601 int j = numbers [target];
02602 wType newUpper = rounding. add_up
02603 (prevUpper,
02604 weights [edgeOffset]);
02605 wType newLower = rounding. add_down
02606 (prevLower,
02607 weights [edgeOffset]);
02608 if (!finite [len]. test (j))
02609 {
02610 finite [len]. set (j);
02611 FUpper [len] [j] = newUpper;
02612 FLower [len] [j] = newLower;
02613 }
02614 else
02615 {
02616 wType &curUpper =
02617 FUpper [len] [j];
02618 if (newUpper < curUpper)
02619 curUpper = newUpper;
02620 wType &curLower =
02621 FLower [len] [j];
02622 if (newLower < curLower)
02623 curLower = newLower;
02624 }
02625 }
02626 }
02627 }
02628
02629
02630 wType minCompWeight (0);
02631 bool firstMin = true;
02632 for (int vert = 0; vert < compSize; ++ vert)
02633 {
02634 if (!finite [compSize]. test (vert))
02635 continue;
02636 bool firstAverage = true;
02637 wType maxAverage = 0;
02638 for (int first = 0; first < compSize; ++ first)
02639 {
02640 if (!finite [first]. test (vert))
02641 continue;
02642 const wType diff = rounding. sub_down
02643 (FLower [compSize] [vert],
02644 FUpper [first] [vert]);
02645 wType average = rounding. div_down
02646 (diff, compSize - first);
02647 if (firstAverage)
02648 {
02649 maxAverage = average;
02650 firstAverage = false;
02651 }
02652 else if (maxAverage < average)
02653 maxAverage = average;
02654 }
02655 if (firstMin || (maxAverage < minCompWeight))
02656 {
02657 if (firstAverage)
02658 throw "Bug in Karp's algorithm";
02659 minCompWeight = maxAverage;
02660 firstMin = false;
02661 }
02662 }
02663
02664
02665 if (!comp || (minCompWeight < minWeight))
02666 minWeight = minCompWeight;
02667 }
02668
02669
02670 delete [] components;
02671 delete [] numbers;
02672 for (int i = 0; i < maxCompSize; ++ i)
02673 {
02674 finite [i]. free ();
02675 delete [] (FUpper [i]);
02676 delete [] (FLower [i]);
02677 }
02678 delete [] finite;
02679 delete [] FUpper;
02680 delete [] FLower;
02681
02682
02683 return minWeight;
02684 }
02685
02686 template <class wType>
02687 template <class arrayType, class roundType>
02688 wType diGraph<wType>::minMeanPathWeight (const roundType &rounding,
02689 const arrayType &starting, int n) const
02690 {
02691
02692 int nIndices = 2;
02693 wType **F = new wType * [nIndices];
02694 BitField *finite = new BitField [nIndices];
02695 for (int i = 0; i < nIndices; ++ i)
02696 {
02697 F [i] = new wType [nVertices];
02698 finite [i]. allocate (nVertices);
02699 }
02700
02701
02702 for (int i = 0; i < n; ++ i)
02703 {
02704 int v = starting [i];
02705 if ((v < 0) || (v >= nVertices))
02706 throw "Starting vertex out of range.";
02707 F [0] [v] = 0;
02708 finite [0]. set (v);
02709 }
02710
02711
02712 wType minWeight (0);
02713 bool firstWeight = true;
02714 for (int len = 1; len <= nVertices; ++ len)
02715 {
02716
02717 int prevIndex = (len - 1) & 1;
02718 int curIndex = len & 1;
02719 finite [curIndex]. clearall (nVertices);
02720
02721
02722 for (int source = 0; source < nVertices; ++ source)
02723 {
02724 if (!finite [prevIndex]. test (source))
02725 continue;
02726 wType prevWeight = F [prevIndex] [source];
02727
02728
02729 int edgeOffset = source ?
02730 edgeEnds [source - 1] : 0;
02731 int edgeEnd = edgeEnds [source];
02732 for (; edgeOffset < edgeEnd; ++ edgeOffset)
02733 {
02734 int target = edges [edgeOffset];
02735 wType newWeight = rounding. add_down
02736 (prevWeight, weights [edgeOffset]);
02737 if (!finite [curIndex]. test (target))
02738 {
02739 finite [curIndex]. set (target);
02740 F [curIndex] [target] = newWeight;
02741 }
02742 else if (newWeight < F [curIndex] [target])
02743 F [curIndex] [target] = newWeight;
02744 }
02745 }
02746
02747
02748 for (int vert = 0; vert < nVertices; ++ vert)
02749 {
02750 if (!finite [curIndex]. test (vert))
02751 continue;
02752 wType average = rounding. div_down
02753 (F [curIndex] [vert], len);
02754 if (firstWeight)
02755 {
02756 minWeight = average;
02757 firstWeight = false;
02758 }
02759 else if (average < minWeight)
02760 minWeight = average;
02761 }
02762 }
02763
02764
02765 for (int i = 0; i < nIndices; ++ i)
02766 {
02767 finite [i]. free ();
02768 delete [] (F [i]);
02769 }
02770 delete [] finite;
02771 delete [] F;
02772
02773
02774 return minWeight;
02775 }
02776
02777 template <class wType>
02778 template <class arrayType>
02779 wType diGraph<wType>::minMeanPathWeight (const arrayType &starting, int n)
02780 const
02781 {
02782 const dummyRounding<wType> rounding = dummyRounding<wType> ();
02783 return minMeanPathWeight (rounding, starting, n);
02784 }
02785
02786
02787
02788
02789
02790
02791
02792 template <class wType>
02793 inline int addRenumEdges (const diGraph<wType> &g, int vertex,
02794 const int *newNum, int cur, int *srcVert, diGraph<wType> &result)
02795 {
02796 int nEdges = g. countEdges (vertex);
02797
02798 for (int edge = 0; edge < nEdges; ++ edge)
02799 {
02800
02801 int dest = newNum [g. getEdge (vertex, edge)];
02802
02803 if (dest == cur)
02804 continue;
02805
02806
02807 if (srcVert [dest] == cur)
02808 continue;
02809
02810
02811 srcVert [dest] = cur;
02812
02813 result. addEdge (dest);
02814 }
02815 return 0;
02816 }
02817
02818
02819
02820
02821
02822
02823 template <class wType, class Table1, class Table2>
02824 inline int collapseVertices (const diGraph<wType> &g,
02825 int compNum, Table1 &compVertices, Table2 &compEnds,
02826 diGraph<wType> &result, int *newNumbers = 0)
02827 {
02828
02829 int nVert = g. countVertices ();
02830 if (!nVert)
02831 return 0;
02832
02833
02834
02835 int *newNum = newNumbers ? newNumbers : new int [nVert];
02836 for (int i = 0; i < nVert; ++ i)
02837 newNum [i] = -1;
02838 int cur = 0;
02839 int pos = 0;
02840 while (cur < compNum)
02841 {
02842 int posEnd = compEnds [cur];
02843 while (pos < posEnd)
02844 newNum [compVertices [pos ++]] = cur;
02845 ++ cur;
02846 }
02847 for (int i = 0; i < nVert; ++ i)
02848 if (newNum [i] < 0)
02849 newNum [i] = cur ++;
02850
02851
02852
02853 int newVert = nVert - compEnds [compNum - 1] + compNum;
02854
02855 if (cur != newVert)
02856 throw "DEBUG: Wrong numbers.";
02857 int *srcVert = new int [newVert];
02858 for (int i = 0; i < newVert; ++ i)
02859 srcVert [i] = -1;
02860
02861
02862
02863 cur = 0, pos = 0;
02864 while (cur < compNum)
02865 {
02866
02867 result. addVertex ();
02868
02869 int posEnd = compEnds [cur];
02870 while (pos < posEnd)
02871 {
02872
02873 int vertex = compVertices [pos ++];
02874
02875 addRenumEdges (g, vertex, newNum, cur,
02876 srcVert, result);
02877 }
02878 ++ cur;
02879 }
02880
02881 for (int vertex = 0; vertex < nVert; ++ vertex)
02882 {
02883
02884 if (newNum [vertex] > cur)
02885 throw "DEBUG: Wrong order.";
02886
02887 if (newNum [vertex] != cur)
02888 continue;
02889
02890 result. addVertex ();
02891 addRenumEdges (g, vertex, newNum, cur, srcVert, result);
02892 ++ cur;
02893 }
02894
02895
02896 delete [] srcVert;
02897 if (!newNumbers)
02898 delete [] newNum;
02899 return 0;
02900 }
02901
02902
02903
02904
02905
02906
02907 template <class wType>
02908 int inclusionGraph (const diGraph<wType> &g, int nVert,
02909 diGraph<wType> &result)
02910 {
02911
02912 int nVertG = g. countVertices ();
02913 if (!nVertG)
02914 return 0;
02915
02916
02917 BitField visited;
02918 visited. allocate (nVertG);
02919 visited. clearall (nVertG);
02920
02921
02922 BitSets lists (nVertG, nVert);
02923
02924
02925 std::stack<int> s_vertex;
02926 std::stack<int> s_edge;
02927
02928
02929 int startVertex = 0;
02930 int vertex = 0;
02931 int edge = 0;
02932 visited. set (vertex);
02933 while (1)
02934 {
02935
02936
02937 if (edge >= g. countEdges (vertex))
02938 {
02939
02940
02941 if (s_vertex. empty ())
02942 {
02943 while ((startVertex < nVertG) &&
02944 (visited. test (startVertex)))
02945 ++ startVertex;
02946 if (startVertex >= nVertG)
02947 break;
02948 vertex = startVertex;
02949 edge = 0;
02950 visited. set (vertex);
02951 continue;
02952 }
02953
02954
02955 int prev = s_vertex. top ();
02956
02957
02958
02959
02960 if (vertex >= nVert)
02961 lists. addset (prev, vertex);
02962
02963
02964 vertex = prev;
02965 s_vertex. pop ();
02966 edge = s_edge. top ();
02967 s_edge. pop ();
02968 continue;
02969 }
02970
02971
02972 int next = g. getEdge (vertex, edge ++);
02973
02974
02975 if (next < nVert)
02976 lists. add (vertex, next);
02977
02978
02979 if (!visited. test (next))
02980 {
02981
02982 s_vertex. push (vertex);
02983 s_edge. push (edge);
02984
02985
02986 vertex = next;
02987 edge = 0;
02988 visited. set (vertex);
02989 }
02990
02991
02992
02993
02994 else if (next >= nVert)
02995 {
02996 lists. addset (vertex, next);
02997 }
02998 }
02999
03000
03001 for (int vertex = 0; vertex < nVert; ++ vertex)
03002 {
03003 result. addVertex ();
03004 int edge = lists. get (vertex);
03005 while (edge >= 0)
03006 {
03007 result. addEdge (edge);
03008 edge = lists. get (vertex, edge + 1);
03009 }
03010 }
03011
03012
03013 visited. free ();
03014 return 0;
03015 }
03016
03017
03018
03019
03020
03021
03022
03023
03024
03025
03026
03027 template<class wType, class TConn>
03028 int inclusionGraph (const diGraph<wType> &g, int nVert,
03029 diGraph<wType> &result, TConn &connections)
03030 {
03031
03032 int nVertG = g. countVertices ();
03033 if (!nVertG)
03034 return 0;
03035
03036
03037 BitField visited;
03038 visited. allocate (nVertG);
03039 visited. clearall (nVertG);
03040
03041
03042 BitSets forwardlists (nVertG, nVert);
03043
03044
03045 std::stack<int> s_vertex;
03046 std::stack<int> s_edge;
03047
03048
03049 int startVertex = 0;
03050 int vertex = 0;
03051 int edge = 0;
03052 visited. set (vertex);
03053 while (1)
03054 {
03055
03056
03057 if (edge >= g. countEdges (vertex))
03058 {
03059
03060
03061 if (s_vertex. empty ())
03062 {
03063 while ((startVertex < nVertG) &&
03064 (visited. test (startVertex)))
03065 ++ startVertex;
03066 if (startVertex >= nVertG)
03067 break;
03068 vertex = startVertex;
03069 edge = 0;
03070 visited. set (vertex);
03071 continue;
03072 }
03073
03074
03075 int prev = s_vertex. top ();
03076
03077
03078
03079
03080 if (vertex >= nVert)
03081 forwardlists. addset (prev, vertex);
03082
03083
03084 vertex = prev;
03085 s_vertex. pop ();
03086 edge = s_edge. top ();
03087 s_edge. pop ();
03088 continue;
03089 }
03090
03091
03092 int next = g. getEdge (vertex, edge ++);
03093
03094
03095 if (next < nVert)
03096 forwardlists. add (vertex, next);
03097
03098
03099 if (!visited. test (next))
03100 {
03101
03102 s_vertex. push (vertex);
03103 s_edge. push (edge);
03104
03105
03106 vertex = next;
03107 edge = 0;
03108 visited. set (vertex);
03109 }
03110
03111
03112
03113
03114 else if (next >= nVert)
03115 {
03116 forwardlists. addset (vertex, next);
03117 }
03118 }
03119
03120
03121
03122
03123 BitSets backlists (nVertG, nVert);
03124
03125 diGraph<wType> gT;
03126 g. transpose (gT);
03127
03128
03129 if (!s_vertex. empty () || !s_edge. empty ())
03130 throw "DEBUG: Nonempty stacks in 'inclusionGraph'.";
03131
03132
03133 visited. clearall (nVertG);
03134
03135
03136 startVertex = 0;
03137 vertex = 0;
03138 edge = 0;
03139 visited. set (vertex);
03140 while (1)
03141 {
03142
03143
03144 if (edge >= gT. countEdges (vertex))
03145 {
03146
03147
03148 if (s_vertex. empty ())
03149 {
03150 while ((startVertex < nVertG) &&
03151 (visited. test (startVertex)))
03152 ++ startVertex;
03153 if (startVertex >= nVertG)
03154 break;
03155 vertex = startVertex;
03156 edge = 0;
03157 visited. set (vertex);
03158 continue;
03159 }
03160
03161
03162 int prev = s_vertex. top ();
03163
03164
03165
03166
03167 if (vertex >= nVert)
03168 backlists. addset (prev, vertex);
03169
03170
03171 vertex = prev;
03172 s_vertex. pop ();
03173 edge = s_edge. top ();
03174 s_edge. pop ();
03175 continue;
03176 }
03177
03178
03179 int next = gT. getEdge (vertex, edge ++);
03180
03181
03182 if (next < nVert)
03183 backlists. add (vertex, next);
03184
03185
03186 if (!visited. test (next))
03187 {
03188
03189 s_vertex. push (vertex);
03190 s_edge. push (edge);
03191
03192
03193 vertex = next;
03194 edge = 0;
03195 visited. set (vertex);
03196 }
03197
03198
03199
03200
03201 else if (next >= nVert)
03202 {
03203 backlists. addset (vertex, next);
03204 }
03205 }
03206
03207
03208
03209
03210 for (int v = nVert; v < nVertG; ++ v)
03211 {
03212 int bvertex = backlists. get (v);
03213 while (bvertex >= 0)
03214 {
03215 int fvertex = forwardlists. get (v);
03216 while (fvertex >= 0)
03217 {
03218 connections (bvertex, fvertex, v);
03219 fvertex = forwardlists. get (v, fvertex + 1);
03220 }
03221 bvertex = backlists. get (v, bvertex + 1);
03222 }
03223 }
03224
03225
03226
03227
03228 for (int vertex = 0; vertex < nVert; ++ vertex)
03229 {
03230 result. addVertex ();
03231 int edge = forwardlists. get (vertex);
03232 while (edge >= 0)
03233 {
03234 result. addEdge (edge);
03235 edge = forwardlists. get (vertex, edge + 1);
03236 }
03237 }
03238
03239
03240 visited. free ();
03241 return 0;
03242 }
03243
03244
03245
03246
03247
03248
03249 template <class wType, class matrix>
03250 void graph2matrix (const diGraph<wType> &g, matrix &m)
03251 {
03252 int nVert = g. countVertices ();
03253 for (int v = 0; v < nVert; ++ v)
03254 {
03255 int nEdges = g. countEdges (v);
03256 for (int e = 0; e < nEdges; ++ e)
03257 {
03258 int w = g. getEdge (v, e);
03259 m [v] [w] = 1;
03260 }
03261 }
03262 return;
03263 }
03264
03265
03266
03267
03268
03269 template <class wType, class matrix>
03270 void matrix2graph (const matrix &m, int n, diGraph<wType> &g)
03271 {
03272 for (int v = 0; v < n; ++ v)
03273 {
03274 g. addVertex ();
03275 for (int w = 0; w < n; ++ w)
03276 if (m [v] [w])
03277 g. addEdge (w);
03278 }
03279 return;
03280 }
03281
03282
03283
03284
03285
03286
03287 template <class matrix>
03288 void transitiveClosure (matrix &m, int n)
03289 {
03290 for (int k = 0; k < n; ++ k)
03291 {
03292 for (int i = 0; i < n; ++ i)
03293 {
03294 for (int j = 0; j < n; ++ j)
03295 {
03296 if (m [i] [k] && m [k] [j])
03297 m [i] [j] = 1;
03298 }
03299 }
03300 }
03301 return;
03302 }
03303
03304
03305
03306
03307
03308
03309
03310
03311 template <class matrix>
03312 void transitiveReduction (matrix &m, int n)
03313 {
03314 for (int k = n - 1; k >= 0; -- k)
03315 {
03316 for (int i = 0; i < n; ++ i)
03317 {
03318 for (int j = 0; j < n; ++ j)
03319 {
03320 if (m [i] [k] && m [k] [j])
03321 m [i] [j] = 0;
03322 }
03323 }
03324 }
03325 return;
03326 }
03327
03328
03329
03330 template <class wType>
03331 void transitiveReduction (const diGraph<wType> &g, diGraph<wType> &gRed)
03332 {
03333 int nVert = g. countVertices ();
03334 if (nVert <= 0)
03335 return;
03336 flatMatrix<char> m (nVert);
03337 graph2matrix (g, m);
03338 transitiveClosure (m, nVert);
03339 transitiveReduction (m, nVert);
03340 matrix2graph (m, nVert, gRed);
03341 return;
03342 }
03343
03344
03345 }
03346 }
03347
03348 #endif // _DIGRAPH_H_
03349
03350
03351