• Main Page
  • Modules
  • Namespaces
  • Classes
  • Files
  • File List
  • File Members

chomp/struct/digraph.h

Go to the documentation of this file.
00001 /// @addtogroup struct
00002 /// @{
00003 
00004 /////////////////////////////////////////////////////////////////////////////
00005 ///
00006 /// @file digraph.h
00007 ///
00008 /// This header file contains the definition of a weighted directed graph
00009 /// class and several algorithms on this graph, including some minimal path
00010 /// algorithms with rounding control to compute rigorous results.
00011 ///
00012 /// @author Pawel Pilarczyk
00013 ///
00014 /////////////////////////////////////////////////////////////////////////////
00015 
00016 // Copyright (C) 1997-2013 by Pawel Pilarczyk.
00017 //
00018 // This file is part of the Homology Library.  This library is free software;
00019 // you can redistribute it and/or modify it under the terms of the GNU
00020 // General Public License as published by the Free Software Foundation;
00021 // either version 2 of the License, or (at your option) any later version.
00022 //
00023 // This library is distributed in the hope that it will be useful,
00024 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00025 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00026 // GNU General Public License for more details.
00027 //
00028 // You should have received a copy of the GNU General Public License along
00029 // with this software; see the file "license.txt".  If not, write to the
00030 // Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
00031 // MA 02111-1307, USA.
00032 
00033 // Started in January 2006. Last revision: May 29, 2010.
00034 
00035 
00036 #ifndef _CHOMP_STRUCT_DIGRAPH_H_
00037 #define _CHOMP_STRUCT_DIGRAPH_H_
00038 
00039 #include <iostream>
00040 #include <new>
00041 #include <stack>
00042 #include <queue>
00043 #include <set>
00044 #include <memory>
00045 #include <vector>
00046 #include <algorithm>
00047 
00048 #include "chomp/struct/multitab.h"
00049 #include "chomp/struct/hashsets.h"
00050 #include "chomp/struct/flatmatr.h"
00051 #include "chomp/struct/bitfield.h"
00052 #include "chomp/struct/bitsets.h"
00053 #include "chomp/struct/fibheap.h"
00054 #include "chomp/struct/autoarray.h"
00055 #include "chomp/system/textfile.h"
00056 #include "chomp/system/timeused.h"
00057 
00058 namespace chomp {
00059 namespace homology {
00060 
00061 
00062 // --------------------------------------------------
00063 // ----------------- DUMMY ROUNDING -----------------
00064 // --------------------------------------------------
00065 
00066 /// A dummy class for rounding operations which does not actually do
00067 /// any rounding. Please, use it as a template for your own rounding classes.
00068 template <class numType>
00069 class dummyRounding
00070 {
00071 public:
00072         /// Adds two numbers with the result rounded downwards.
00073         static inline numType add_down (const numType &x, const numType &y)
00074         { return x + y; }
00075 
00076         /// Adds two numbers with the result rounded upwards.
00077         static inline numType add_up (const numType &x, const numType &y)
00078         { return x + y; }
00079 
00080         /// Subtracts two numbers with the result rounded downwards.
00081         static inline numType sub_down (const numType &x, const numType &y)
00082         { return x - y; }
00083 
00084         /// Subtracts two numbers with the result rounded upwards.
00085         static inline numType sub_up (const numType &x, const numType &y)
00086         { return x - y; }
00087 
00088         /// Multiplies two numbers with the result rounded downwards.
00089         static inline numType mul_down (const numType &x, const numType &y)
00090         { return x * y; }
00091 
00092         /// Multiplies two numbers with the result rounded upwards.
00093         static inline numType mul_up (const numType &x, const numType &y)
00094         { return x * y; }
00095 
00096         /// Divides two numbers with the result rounded downwards.
00097         static inline numType div_down (const numType &x, const numType &y)
00098         { return x / y; }
00099 
00100         /// Divides a number by an integer with the result rounded downwards.
00101         static inline numType div_down (const numType &x, int_t y)
00102         { return x / y; }
00103 
00104         /// Divides two numbers with the result rounded upwards.
00105         static inline numType div_up (const numType &x, const numType &y)
00106         { return x / y; }
00107 
00108 private:
00109 }; /* class dummyRounding */
00110 
00111 
00112 // --------------------------------------------------
00113 // ------------------ DUMMY ARRAY -------------------
00114 // --------------------------------------------------
00115 
00116 /// A dummy array of integers that ignores all the assigned values.
00117 class dummyArray
00118 {
00119 public:
00120         /// The constructor of a dummy array.
00121         dummyArray () {n = 0;}
00122 
00123         /// Operator [] which returns a reference to a dummy value.
00124         int_t &operator [] (int_t) {return n;}
00125 
00126 private:
00127         /// A dummy number used as a dump box for assigning values.
00128         int_t n;
00129 
00130 }; /* class dummyArray */
00131 
00132 
00133 // --------------------------------------------------
00134 // ----------------- DIRECTED GRAPH -----------------
00135 // --------------------------------------------------
00136 
00137 // #define DIGRAPH_DEBUG
00138 
00139 /// This class defines a directed graph with very limited number of
00140 /// operations, and a few specific algorithms implemented on it, like DFS.
00141 /// The graph can be treated as weighted if necessary.
00142 template <class wType = int>
00143 class diGraph
00144 {
00145 public:
00146         /// The type of the weight of edges.
00147         typedef wType weight_type;
00148         
00149         /// The default constructor of an empty graph.
00150         /// Note: The default copy constructor and assignment operator
00151         /// generated by the compiler are good for copying graphs.
00152         diGraph ();
00153 
00154         /// The destructor.
00155         ~diGraph ();
00156 
00157         /// Swaps the data with another graph.
00158         void swap (diGraph<wType> &g);
00159 
00160         /// Adds a vertex.
00161         void addVertex (void);
00162 
00163         /// Adds an edge starting at the last vertex.
00164         /// Note: The range of the target vertex number is not verified.
00165         void addEdge (int_t target);
00166 
00167         /// Adds an edge from the last vertex to the given one
00168         /// and sets the weight of this edge.
00169         void addEdge (int_t target, const wType &weight);
00170 
00171         /// Sets the weight of the given edge.
00172         void setWeight (int_t vertex, int_t i, const wType &weight);
00173 
00174         /// Sets the weight of the given edge.
00175         void setWeight (int_t edge, const wType &weight);
00176         
00177         /// Sets the weights of all the edges at a time.
00178         /// The weights are pulled from the table with the [] operator.
00179         template <class Table>
00180         void setWeights (const Table &tab);
00181 
00182         /// Removes the last vertex and all the edges going out from it.
00183         /// This is done in the time O(1).
00184         void removeVertex (void);
00185 
00186         /// Removes the given vertex and all the edges going out from it,
00187         /// as well as the edges going towards it.
00188         /// If requested, the weights in the graph are also updated.
00189         /// This function might be slow - it is done in the time O(V+E).
00190         void removeVertex (int_t vertex, bool updateweights = false);
00191 
00192         /// Returns the number of vertices.
00193         int_t countVertices (void) const;
00194 
00195         /// Returns the number of edges.
00196         int_t countEdges (void) const;
00197 
00198         /// Counts the number of edges leaving the given vertex.
00199         int_t countEdges (int_t vertex) const;
00200         
00201         /// Retrieves the given edge that leaves the given vertex.
00202         int_t getEdge (int_t vertex, int_t i) const;
00203 
00204         /// Retrieves the weight of the given edge.
00205         const wType &getWeight (int_t vertex, int_t i) const;
00206 
00207         /// Retrieves the weight of the given edge.
00208         const wType &getWeight (int_t edge) const;
00209 
00210         /// Gets the weights of all the edges at a time.
00211         /// The weights are put into the table with the [] operator.
00212         template <class Table>
00213         void getWeights (Table &tab) const;
00214 
00215         /// Fills out a table that represents all the edges of the graph.
00216         /// The indices of a starting and ending vertex of the n-th edge
00217         /// are written to "tab [n] [0]" and "tab [n] [1]", respectively.
00218         template <class Table>
00219         void writeEdges (Table &tab) const;
00220 
00221         /// Creates a transposed graph.
00222         template <class wType1>
00223         void transpose (diGraph<wType1> &result,
00224                 bool copyweights = false) const;
00225 
00226         /// Computes a restriction of the graph to its subgraph. The subgraph
00227         /// vertices are defined by nonzero entries in the supplied table.
00228         /// The result must be initially empty.
00229         template <class Table, class wType1>
00230         void subgraph (diGraph<wType1> &result, const Table &tab,
00231                 bool copyweights = false) const;
00232 
00233         /// Marks each vertex visited by DFS with the given color,
00234         /// starting with the given vertex. Runs for one component only.
00235         /// The initial color in 'tab' must be different than the given one.
00236         template <class Table, class Color>
00237         void DFScolor (Table &tab, const Color &color,
00238                 int_t vertex = 0) const;
00239 
00240         /// The recurrent procedure for DFScolor.
00241         template <class Table, class Color>
00242         void DFScolorRecurrent (Table &tab, const Color &color,
00243                 int_t vertex = 0) const;
00244 
00245         /// A stack version of the recurrent procedure for DFScolor.
00246         template <class Table, class Color>
00247         void DFScolorStack (Table &tab, const Color &color,
00248                 int_t vertex = 0) const;
00249 
00250         /// Computes the DFS finishing time for each vertex.
00251         /// Note: The time begins with 1, not with 0.
00252         template <class Table>
00253         void DFSfinishTime (Table &tab) const;
00254 
00255         /// Computes the DFS forest.
00256         /// Considers the vertices in the given order.
00257         /// Saves the numbers of vertices of each tree in 'compVertices',
00258         /// and keeps the one-beyond-the-end offsets of the trees in the
00259         /// table 'compEnds'. Records the connections between the trees
00260         /// in 'scc' (which must be empty when this function is called).
00261         /// If requested, only those single-vertex trees are counted
00262         /// which have an edge that loops back to themselves.
00263         /// Returns the number of trees in the computed forest.
00264         template <class Table1, class Table2, class Table3>
00265         int_t DFSforest (const Table1 &ordered, Table2 &compVertices,
00266                 Table3 &compEnds, bool nontrivial = false,
00267                 diGraph<wType> *sccGraph = 0) const;
00268 
00269         /// Computes the length of the shortest nontrivial path
00270         /// from the given vertex to another one.
00271         /// The length is measured by counting edges.
00272         /// Uses a stack version of the BFS algorithm.
00273         /// Returns the length of the path or 0 if none.
00274         int_t shortestPath (int_t source, int_t destination) const;
00275 
00276         /// Computes the length of the shortest loop from the given vertex
00277         /// to itself. The length is measured by counting edges on the way.
00278         /// Uses a stack version of the BFS algorithm.
00279         /// Returns the length of the loop or 0 if none.
00280         int_t shortestLoop (int_t origin) const;
00281 
00282         /// Dijkstra's algorithm for solving the single-source shortest
00283         /// paths problem if all the edge weights are nonnegative.
00284         /// The table 'len' is used to store the path lengths during the
00285         /// computations and contains the final result. A negative value
00286         /// stands for the infinity (no path to the given vertex).
00287         /// This is a special version that uses the given edge weights
00288         /// instead of the weights contained in the definition of the graph.
00289         template <class lenTable, class weightsType, class roundType>
00290         void Dijkstra (const roundType &rounding, int_t source,
00291                 lenTable &len, weightsType &edgeWeights) const;
00292 
00293         /// Dijkstra's algorithm running on the graph's own weights.
00294         template <class lenTable, class roundType>
00295         void Dijkstra (const roundType &rounding, int_t source,
00296                 lenTable &len) const;
00297 
00298         /// The above algorithm without rounding control.
00299         template <class lenTable>
00300         void Dijkstra (int_t source, lenTable &len) const;
00301 
00302         /// Runs the Bellman-Ford algorithm which computes the single-source
00303         /// shortest paths in a weighted directed graph, where some edge
00304         /// weights may be negative. Runs in O(V*E).
00305         /// The table 'len' is used to store the path lengths during the
00306         /// computations and contains the final result. The number for
00307         /// infinity is set to indicate unreachable vertices.
00308         /// The table with predecessors allows to retrieve shortest paths;
00309         /// this must be a pointer to an array-like structure, e.g., int **.
00310         /// To ignore this data, use an object of the 'dummyArray' class.
00311         /// Returns true if successful, false if a negative cycle is found.
00312         template <class lenTable, class predTable, class roundType>
00313         bool BellmanFord (const roundType &rounding, int_t source,
00314                 lenTable &len, wType *infinity, predTable pred) const;
00315 
00316         /// The above algorithm without rounding control.
00317         template <class lenTable, class predTable>
00318         bool BellmanFord (int_t source, lenTable &len, wType *infinity,
00319                 predTable pred) const;
00320 
00321         /// Runs the Bellman-Ford algorithm (see above) without storing the
00322         /// distances, only returns the information about the existence
00323         /// of a negative-weight cycle.
00324         template <class roundType>
00325         bool BellmanFord (const roundType &rounding, int_t source) const;
00326 
00327         /// The above algorithm without rounding control.
00328         bool BellmanFord (int_t source) const;
00329 
00330         /// Runs the Edmonds algorithm to compute the shortest path
00331         /// that runs through all the vertices of the graph.
00332         /// Computation time: O (n log n). The length of the path
00333         /// is measured as the sum of the weights of the edges.
00334         /// The path does not contain any loops.
00335         /// The graph should be strongly connected.
00336         wType Edmonds () const;
00337 
00338         /// An old implementation of the Edmonds algorithm (less efficient).
00339         wType EdmondsOld () const;
00340 
00341         /// Runs the Floyd-Warshall algorithm to compute the shortest
00342         /// paths between all pairs of vertices in the graph.
00343         /// The position [i] [j] of the array contains
00344         /// the length of the shortest path from vertex i to vertex j.
00345         /// Provides a rigorous lower bound in interval arithmetic,
00346         /// provided that intervals are compared with "<" and "<="
00347         /// by comparing their lower ends only.
00348         /// If "setInfinity" is "true", then computes a value that serves
00349         /// as the infinity, fills in the corresponding entries in "arr",
00350         /// and returns this value.
00351         /// Otherwise, returns the length of the shortest path. In this
00352         /// case, arr [i] [j] is undefined if there is no path i -> j.
00353         /// Throws an error message if a negative loop is found,
00354         /// unless "ignoreNegLoop" is set to "true".
00355         template <class arrayType, class roundType>
00356         wType FloydWarshall (const roundType &rounding, arrayType &arr,
00357                 bool setInfinity = true, bool ignoreNegLoop = false) const;
00358 
00359         /// The above algorithm without rounding control.
00360         template <class arrayType>
00361         wType FloydWarshall (arrayType &arr, bool setInfinity = true,
00362                 bool ignoreNegLoop = false) const;
00363 
00364         /// Runs Johnson's algorithm to compute the minimum path weight
00365         /// between any vertices in the graph. The time complexity of this
00366         /// algorithm is essentially O (V^2 log V + VE log V),
00367         /// which is better than the complexity of the Warshall-Floyd
00368         /// algorithm for sparse graphs, that is, graphs in which the number
00369         /// of edges E is of order smaller than V^2.
00370         /// The meaning of the arguments and the returned value is the same
00371         /// as in 'FloydWarshall'.
00372         template <class arrayType, class roundType>
00373         wType Johnson (const roundType &rounding, arrayType &arr,
00374                 bool setInfinity = true, bool ignoreNegLoop = false) const;
00375 
00376         /// The above algorithm without rounding control.
00377         template <class arrayType>
00378         wType Johnson (arrayType &arr, bool setInfinity = true,
00379                 bool ignoreNegLoop = false) const;
00380 
00381         /// Uses the Floyd-Warshall algorithm or Johnson's algorithm,
00382         /// depending on the number of edges, to compute the minimum
00383         /// path weight between any vertices in the graph.
00384         /// Throws an error message if a negative loop exists in the graph,
00385         /// unless "ignoreNegLoop" is set to "true".
00386         /// To force the use of Johnson's algorithm, set "sparseGraph" to 1,
00387         /// to force the use of Warshall-Floyd, set "sparseGraph" to 0,
00388         /// otherwise it will be determined heuristically which algorithm
00389         /// should be used.
00390         template <class roundType>
00391         wType minPathWeight (const roundType &rounding,
00392                 bool ignoreNegLoop = false, int sparseGraph = -1) const;
00393 
00394         /// The above algorithm without rounding control.
00395         wType minPathWeight (bool ignoreNegLoop = false,
00396                 int sparseGraph = -1) const;
00397 
00398         /// Runs the Karp algorithm for each strongly connected component
00399         /// of the graph and returns the minimum mean cycle weight,
00400         /// which can be negative.
00401         /// As a byproduct, saves the transposed graph, if requested to.
00402         wType minMeanCycleWeight (diGraph<wType> *transposed = 0) const;
00403 
00404         /// A version of the above function modified for the purpose
00405         /// of interval arithmetic to provide the correct lower bound
00406         /// for the minimum mean cycle weight in a graph.
00407         /// This specialization is necessary, because of the subtraction
00408         /// operation used in Karp's algorithm. Therefore, upper and lower
00409         /// bounds for the minimum path progression weights must be computed
00410         /// independently.
00411         /// The intervals are compared by comparing their lower bounds only.
00412         template <class roundType>
00413         wType minMeanCycleWeight (const roundType &rounding,
00414                 diGraph<wType> *transposed) const;
00415 
00416         /// Runs an algorithm based on Karp's idea to compute the minimum
00417         /// mean path weight for paths starting at any of the given
00418         /// n vertices and of length not exceeding the number of vertices
00419         /// in the graph.
00420         /// Returns 0 if no path starts at any of the vertices.
00421         template <class arrayType, class roundType>
00422         wType minMeanPathWeight (const roundType &rounding,
00423                 const arrayType &starting, int_t n) const;
00424 
00425         /// The above algorithm without rounding control.
00426         template <class arrayType>
00427         wType minMeanPathWeight (const arrayType &starting, int_t n) const;
00428 
00429         /// Operator == for comparing digraphs. No isomorphism check, just
00430         /// comparing with the same order of vertices. Ignores weights.
00431         template <class wType1, class wType2>
00432         friend bool operator == (const diGraph<wType1> &g1,
00433                 const diGraph<wType2> &g2);
00434 
00435         /// Outputs the graph to a text stream in a human-readable format.
00436         template <class outType>
00437         outType &show (outType &out, bool showWeights = false) const;
00438 
00439 protected:
00440         /// The number of vertices.
00441         int_t nVertices;
00442 
00443         /// A table with the offsets of the one-after-the-last edge
00444         /// of each vertex.
00445         multitable<int_t> edgeEnds;
00446 
00447         /// A table with edge target numbers.
00448         multitable<int_t> edges;
00449 
00450         /// A table with edge weights.
00451         multitable<wType> weights;
00452 
00453 private:
00454         /// The recurrent procedure for DFSfinishTime.
00455         template <class Table>
00456         void DFSfinishTimeRecurrent (Table &tab, int_t vertex,
00457                 int_t &counter) const;
00458 
00459         /// A stack version of the recurrent procedure for DFSfinishTime.
00460         template <class Table>
00461         void DFSfinishTimeStack (Table &tab, int_t vertex,
00462                 int_t &counter) const;
00463 
00464         /// The recurrent procedure for DFSforest.
00465         /// Returns true iff there is a loop within the tree found.
00466         template <class Table1, class Table2>
00467         bool DFSforestRecurrent (Table1 &tab, Table1 &ntab, int_t vertex,
00468                 int_t treeNumber, int_t countTrees, Table2 &compVertices,
00469                 int_t &curVertex, diGraph *sccGraph, int_t *sccEdgeAdded)
00470                 const;
00471 
00472         /// A stack version of the recurrent procedure for DFSforest.
00473         template <class Table1, class Table2>
00474         bool DFSforestStack (Table1 &tab, Table1 &ntab, int_t vertex,
00475                 int_t treeNumber, int_t countTrees, Table2 &compVertices,
00476                 int_t &curVertex, diGraph *sccGraph, int_t *sccEdgeAdded)
00477                 const;
00478 
00479         /// An edge with a weight (used by the Edmonds algorithm).
00480         struct edgeTriple
00481         {
00482                 /// The starting and ending vertices of the edge.
00483                 int_t vert1;
00484                 int_t vert2;
00485                 
00486                 /// The weight of the edge.
00487                 wType weight;
00488 
00489                 /// The < operator for comparing the weights of edges.
00490                 friend bool operator < (const edgeTriple &x,
00491                         const edgeTriple &y)
00492                 {
00493                         return (x. weight < y. weight);
00494                 }
00495         };
00496 
00497         /// A class for representing a positive number with negative values
00498         /// serving as the infinity (used in the Dijkstra algorithm).
00499         class posWeight
00500         {
00501         public:
00502                 /// The default constructor.
00503                 posWeight ()
00504                 {
00505                         value = 0;
00506                         return;
00507                 }
00508 
00509                 /// An optional constructor.
00510                 explicit posWeight (const wType &_value)
00511                 {
00512                         if (_value < 0)
00513                                 value = 0;
00514                         else
00515                                 value = _value;
00516                         return;
00517                 }
00518 
00519                 /// Sets the value to the infinity.
00520                 void setInfinity ()
00521                 {
00522                         value = -1;
00523                         return;
00524                 }
00525 
00526                 /// Returns true iff the value corresponds to the infinity.
00527                 bool isInfinity () const
00528                 {
00529                         return (value < 0);
00530                 }
00531 
00532                 /// Returns the value stored in this object.
00533                 const wType &getValue () const
00534                 {
00535                         return value;
00536                 }
00537 
00538                 /// The < operator for comparing the numbers.
00539                 friend bool operator < (const posWeight &x,
00540                         const posWeight &y)
00541                 {
00542                         if (y. isInfinity ())
00543                                 return !(x. isInfinity ());
00544                         else if (x. isInfinity ())
00545                                 return false;
00546                         return (x. value < y. value);
00547                 }
00548 
00549                 /// The operator for showing the number to the output stream.
00550                 friend std::ostream &operator << (std::ostream &out,
00551                         const posWeight &x)
00552                 {
00553                         if (x. isInfinity ())
00554                                 out << "+oo";
00555                         else
00556                                 out << x. getValue ();
00557                         return out;
00558                 }
00559 
00560         private:
00561                 /// The actual number.
00562                 wType value;
00563         };
00564 
00565 }; /* class diGraph */
00566 
00567 // --------------------------------------------------
00568 
00569 template <class wType>
00570 inline diGraph<wType>::diGraph (): nVertices (0),
00571         edgeEnds (1024), edges (4096), weights (4096)
00572 {
00573         return;
00574 } /* diGraph::diGraph */
00575 
00576 template <class wType>
00577 inline diGraph<wType>::~diGraph ()
00578 {
00579         return;
00580 } /* diGraph::~diGraph */
00581 
00582 // --------------------------------------------------
00583 
00584 template <class wType1, class wType2>
00585 inline bool operator == (const diGraph<wType1> &g1,
00586         const diGraph<wType2> &g2)
00587 {
00588         if (g1. nVertices != g2. nVertices)
00589                 return false;
00590         if (!g1. nVertices)
00591                 return true;
00592         for (int_t i = 0; i < g1. nVertices; ++ i)
00593         {
00594                 if (g1. edgeEnds [i] != g2. edgeEnds [i])
00595                         return false;
00596         }
00597         int_t nEdges = g1. edgeEnds [g1. nVertices - 1];
00598         for (int_t i = 0; i < nEdges; ++ i)
00599         {
00600                 if (g1. edges [i] != g2. edges [i])
00601                         return false;
00602         }
00603         return true;
00604 } /* operator == */
00605 
00606 template <class wType1, class wType2>
00607 inline bool operator != (const diGraph<wType1> &g1,
00608         const diGraph<wType2> &g2)
00609 {
00610         return !(g1 == g2);
00611 } /* operator != */
00612 
00613 // --------------------------------------------------
00614 
00615 template <class wType>
00616 inline void diGraph<wType>::swap (diGraph<wType> &g)
00617 {
00618         std::swap (nVertices, g. nVertices);
00619         edgeEnds. swap (g. edgeEnds);
00620         edges. swap (g. edges);
00621         weights. swap (g. weights);
00622         return;
00623 } /* diGraph::swap */
00624 
00625 template <class wType>
00626 inline void diGraph<wType>::addVertex (void)
00627 {
00628         edgeEnds [nVertices] = nVertices ? edgeEnds [nVertices - 1] :
00629                 static_cast<int_t> (0);
00630         ++ nVertices;
00631         return;
00632 } /* diGraph::addVertex */
00633 
00634 template <class wType>
00635 inline void diGraph<wType>::addEdge (int_t target)
00636 {
00637         if (!nVertices)
00638                 throw "Trying to add an edge to an empty graph.";
00639 //      if (target >= nVertices)
00640 //              throw "Trying to add an edge to a nonexistent vertex.";
00641         if (target < 0)
00642                 throw "Trying to add an edge to a negative vertex.";
00643         int_t &offset = edgeEnds [nVertices - 1];
00644         if (offset + 1 <= 0)
00645                 throw "Too many edges in a diGraph (limit = 2,147,483,647).";
00646         edges [offset] = target;
00647         ++ offset;
00648         return;
00649 } /* diGraph::addEdge */
00650 
00651 template <class wType>
00652 inline void diGraph<wType>::addEdge (int_t target, const wType &weight)
00653 {
00654         if (!nVertices)
00655                 throw "Trying to add an edge to an empty graph.";
00656 //      if (target >= nVertices)
00657 //              throw "Trying to add an edge to a nonexistent vertex.";
00658         if (target < 0)
00659                 throw "Trying to add an edge to a negative vertex.";
00660         int_t &offset = edgeEnds [nVertices - 1];
00661         if (offset + 1 <= 0)
00662                 throw "Too many edges in a diGraph (limit = 2,147,483,647).";
00663         edges [offset] = target;
00664         weights [offset] = weight;
00665         ++ offset;
00666         return;
00667 } /* diGraph::addEdge */
00668 
00669 template <class wType>
00670 inline void diGraph<wType>::setWeight (int_t vertex, int_t i,
00671         const wType &weight)
00672 {
00673         if (!vertex)
00674                 weights [i] = weight;
00675         else
00676                 weights [edgeEnds [vertex - 1] + i] = weight;
00677         return;
00678 } /* diGraph::setWeight */
00679 
00680 template <class wType>
00681 inline void diGraph<wType>::setWeight (int_t edge, const wType &weight)
00682 {
00683         weights [edge] = weight;
00684         return;
00685 } /* diGraph::setWeight */
00686 
00687 template <class wType> template <class Table>
00688 inline void diGraph<wType>::setWeights (const Table &tab)
00689 {
00690         if (!nVertices)
00691                 return;
00692         int_t nEdges = edgeEnds [nVertices - 1];
00693         for (int_t i = 0; i < nEdges; ++ i)
00694                 weights [i] = tab [i];
00695         return;
00696 } /* diGraph::setWeights */
00697 
00698 template <class wType>
00699 inline void diGraph<wType>::removeVertex (void)
00700 {
00701         if (!nVertices)
00702                 throw "Trying to remove a vertex from the empty graph.";
00703         -- nVertices;
00704         return;
00705 } /* diGraph::removeVertex */
00706 
00707 template <class wType>
00708 inline void diGraph<wType>::removeVertex (int_t vertex, bool updateweights)
00709 {
00710         // make sure that the vertex number is within the scope
00711         if ((vertex < 0) || (vertex >= nVertices))
00712                 throw "Trying to remove a vertex that is not in the graph.";
00713 
00714         // remove edges that begin or end at the given vertex
00715         int_t curEdge = 0;
00716         int_t newEdge = 0;
00717         int_t skipped = 0;
00718         for (int_t v = 0; v < nVertices; ++ v)
00719         {
00720                 // skip the edges that begin at the given vertex
00721                 if (!skipped && (v == vertex))
00722                 {
00723                         curEdge = edgeEnds [v];
00724                         ++ skipped;
00725                         continue;
00726                 }
00727 
00728                 // skip the edges that point to the given vertex
00729                 int_t maxEdge = edgeEnds [v];
00730                 for (; curEdge < maxEdge; ++ curEdge)
00731                 {
00732                         if (edges [curEdge] == vertex)
00733                                 continue;
00734                         int_t target = edges [curEdge];
00735                         edges [newEdge] = (target < vertex) ? target :
00736                                 (target - 1);
00737                         if (updateweights)
00738                                 weights [newEdge] = weights [curEdge];
00739                         ++ newEdge;
00740                 }
00741                 edgeEnds [v - skipped] = newEdge;
00742         }
00743 
00744         // decrease the number of vertices
00745         nVertices -= skipped;
00746 
00747         return;
00748 } /* diGraph::removeVertex */
00749 
00750 template <class wType>
00751 inline int_t diGraph<wType>::countVertices (void) const
00752 {
00753         return nVertices;
00754 } /* diGraph::countVertices */
00755 
00756 template <class wType>
00757 inline int_t diGraph<wType>::countEdges (void) const
00758 {
00759         if (!nVertices)
00760                 return 0;
00761         else
00762                 return edgeEnds [nVertices - 1];
00763 } /* diGraph::countEdges */
00764 
00765 template <class wType>
00766 inline int_t diGraph<wType>::countEdges (int_t vertex) const
00767 {
00768         if (!vertex)
00769                 return edgeEnds [0];
00770         else
00771                 return edgeEnds [vertex] - edgeEnds [vertex - 1];
00772 } /* diGraph::countEdges */
00773 
00774 template <class wType>
00775 inline int_t diGraph<wType>::getEdge (int_t vertex, int_t i) const
00776 {
00777         if (!vertex)
00778                 return edges [i];
00779         else
00780                 return edges [edgeEnds [vertex - 1] + i];
00781 } /* diGraph::getEdge */
00782 
00783 template <class wType>
00784 inline const wType &diGraph<wType>::getWeight (int_t vertex, int_t i) const
00785 {
00786         if (!vertex)
00787                 return weights [i];
00788         else
00789                 return weights [edgeEnds [vertex - 1] + i];
00790 } /* diGraph::getWeight */
00791 
00792 template <class wType>
00793 inline const wType &diGraph<wType>::getWeight (int_t edge) const
00794 {
00795         return weights [edge];
00796 } /* diGraph::getWeight */
00797 
00798 template <class wType> template <class Table>
00799 inline void diGraph<wType>::getWeights (Table &tab) const
00800 {
00801         if (!nVertices)
00802                 return;
00803         int_t nEdges = edgeEnds [nVertices - 1];
00804         for (int_t i = 0; i < nEdges; ++ i)
00805                 tab [i] = weights [i];
00806         return;
00807 } /* diGraph::getWeights */
00808 
00809 template <class wType> template <class Table>
00810 inline void diGraph<wType>::writeEdges (Table &tab) const
00811 {
00812         int_t curEdge = 0;
00813         for (int_t curVertex = 0; curVertex < nVertices; ++ curVertex)
00814         {
00815                 int_t maxEdge = edgeEnds [curVertex];
00816                 while (curEdge < maxEdge)
00817                 {
00818                         tab [curEdge] [0] = curVertex;
00819                         tab [curEdge] [1] = edges [curEdge];
00820                         ++ curEdge;
00821                 }
00822         }
00823         return;
00824 } /* diGraph::writeEdges */
00825 
00826 template <class wType> template <class wType1>
00827 inline void diGraph<wType>::transpose (diGraph<wType1> &result,
00828         bool copyweights) const
00829 {
00830         // prepare the resulting graph
00831         result. nVertices = nVertices;
00832         if (!nVertices)
00833                 return;
00834 
00835         // compute the initial offsets for the edge numbers
00836         for (int_t i = 0; i < nVertices; ++ i)
00837                 result. edgeEnds [i] = 0;
00838         int_t nEdges = edgeEnds [nVertices - 1];
00839         for (int_t i = 0; i < nEdges; ++ i)
00840         {
00841                 if (edges [i] < nVertices - 1)
00842                         ++ result. edgeEnds [edges [i] + 1];
00843         }
00844         for (int_t i = 2; i < nVertices; ++ i)
00845                 result. edgeEnds [i] += result. edgeEnds [i - 1];
00846 
00847         // compute the reversed edges
00848         int_t curEdge = 0;
00849         for (int_t i = 0; i < nVertices; ++ i)
00850         {
00851                 for (; curEdge < edgeEnds [i]; ++ curEdge)
00852                 {
00853                         int_t j = edges [curEdge];
00854                         int_t &offset = result. edgeEnds [j];
00855                         result. edges [offset] = i;
00856                         if (copyweights)
00857                                 result. weights [offset] = weights [curEdge];
00858                         ++ offset;
00859                 }
00860         }
00861         return;
00862 } /* diGraph::transpose */
00863 
00864 template <class wType> template <class Table, class wType1>
00865 inline void diGraph<wType>::subgraph (diGraph<wType1> &g,
00866         const Table &tab, bool copyweights) const
00867 {
00868         // compute the new numbers of vertices that remain in the graph
00869         int_t *numbers = new int_t [nVertices];
00870         int_t curNumber = 0;
00871         for (int_t i = 0; i < nVertices; ++ i)
00872         {
00873                 if (tab [i])
00874                         numbers [i] = curNumber ++;
00875                 else
00876                         numbers [i] = -1;
00877         }
00878 
00879         // copy the edges from the old graph to the new one
00880         for (int_t i = 0; i < nVertices; ++ i)
00881         {
00882                 if (numbers [i] < 0)
00883                         continue;
00884                 g. addVertex ();
00885                 int_t firstEdge = i ? edgeEnds [i - 1] :
00886                         static_cast<int_t> (0);
00887                 int_t endEdge = edgeEnds [i];
00888                 for (int_t j = firstEdge; j < endEdge; ++ j)
00889                 {
00890                         int_t edgeEnd = edges [j];
00891                         if (numbers [edgeEnd] >= 0)
00892                         {
00893                                 if (copyweights)
00894                                         g. addEdge (numbers [edgeEnd],
00895                                                 weights [j]);
00896                                 else
00897                                         g. addEdge (numbers [edgeEnd]);
00898                         }
00899                 }
00900         }
00901 
00902         // clean up memory and exit
00903         delete [] numbers;
00904         return;
00905 } /* diGraph::subgraph */
00906 
00907 // --------------------------------------------------
00908 
00909 template <class wType> template <class Table, class Color>
00910 inline void diGraph<wType>::DFScolorRecurrent (Table &tab,
00911         const Color &color, int_t vertex) const
00912 {
00913         tab [vertex] = color;
00914         int_t maxEdge = edgeEnds [vertex];
00915         for (int_t i = vertex ? edgeEnds [vertex - 1] :
00916                 static_cast<int_t> (0); i < maxEdge; ++ i)
00917         {
00918                 int_t next = edges [i];
00919                 if (tab [next] != color)
00920                         DFScolorRecurrent (tab, color, next);
00921         }
00922         return;
00923 } /* diGraph::DFScolorRecurrent */
00924 
00925 template <class wType> template <class Table, class Color>
00926 inline void diGraph<wType>::DFScolorStack (Table &tab, const Color &color,
00927         int_t vertex) const
00928 {
00929         // prepare stacks for the recursion
00930         std::stack<int_t> s_vertex;
00931         std::stack<int_t> s_edge;
00932         std::stack<int_t> s_maxedge;
00933 
00934         // mark the current vertex as visited
00935         tab [vertex] = color;
00936 
00937         // determine the edges to be visited
00938         int_t edge = vertex ? edgeEnds [vertex - 1] :
00939                 static_cast<int_t> (0);
00940         int_t maxedge = edgeEnds [vertex];
00941 
00942         while (1)
00943         {
00944                 // return to the previous recursion level
00945                 // if all the edges have been followed
00946                 if (edge >= maxedge)
00947                 {
00948                         // return if this is the initial recursion level
00949                         if (s_vertex. empty ())
00950                                 return;
00951 
00952                         // restore the variables from the previous level
00953                         vertex = s_vertex. top ();
00954                         s_vertex. pop ();
00955                         edge = s_edge. top ();
00956                         s_edge. pop ();
00957                         maxedge = s_maxedge. top ();
00958                         s_maxedge. pop ();
00959                         continue;
00960                 }
00961 
00962                 // go to the deeper recursion level if possible
00963                 int_t next = edges [edge ++];
00964                 if (tab [next] != color)
00965                 {
00966                         // store the previous variables at the stacks
00967                         s_vertex. push (vertex);
00968                         s_edge. push (edge);
00969                         s_maxedge. push (maxedge);
00970 
00971                         // set the new vertex
00972                         vertex = next;
00973                         
00974                         // mark the new vertex as visited
00975                         tab [vertex] = color;
00976                         
00977                         // determine the edges to be visited
00978                         edge = vertex ? edgeEnds [vertex - 1] :
00979                                 static_cast<int_t> (0);
00980                         maxedge = edgeEnds [vertex];
00981                 }
00982         }
00983 } /* diGraph::DFScolorStack */
00984 
00985 template <class wType> template <class Table, class Color>
00986 inline void diGraph<wType>::DFScolor (Table &tab, const Color &color,
00987         int_t vertex) const
00988 {
00989         DFScolorStack (tab, color, vertex);
00990         return;
00991 } /* diGraph::DFScolor */
00992 
00993 // --------------------------------------------------
00994 
00995 template <class wType> template <class Table>
00996 inline void diGraph<wType>::DFSfinishTimeRecurrent (Table &tab,
00997         int_t vertex, int_t &counter) const
00998 {
00999         // mark the current vertex as visited
01000         tab [vertex] = -1;
01001 
01002         // call DFS for the other vertices
01003         for (int_t edge = vertex ? edgeEnds [vertex - 1] :
01004                 static_cast<int_t> (0);
01005                 edge < edgeEnds [vertex]; ++ edge)
01006         {
01007                 int_t next = edges [edge];
01008                 if (!tab [next])
01009                         DFSfinishTimeRecurrent (tab, next, counter);
01010         }
01011 
01012         // record the finishing time for the current vertex and return
01013         tab [vertex] = ++ counter;
01014         return;
01015 } /* diGraph::DFSfinishTimeRecurrent */
01016 
01017 template <class wType> template <class Table>
01018 inline void diGraph<wType>::DFSfinishTimeStack (Table &tab, int_t vertex,
01019         int_t &counter) const
01020 {
01021         // prepare stacks for the recursion
01022         std::stack<int_t> s_vertex;
01023         std::stack<int_t> s_edge;
01024         std::stack<int_t> s_maxedge;
01025 
01026         // mark the current vertex as visited
01027         tab [vertex] = -1;
01028 
01029         // determine the edges to be visited
01030         int_t edge = vertex ? edgeEnds [vertex - 1] :
01031                 static_cast<int_t> (0);
01032         int_t maxedge = edgeEnds [vertex];
01033 
01034         while (1)
01035         {
01036                 // return to the previous recursion level
01037                 // if all the edges have been followed
01038                 if (edge >= maxedge)
01039                 {
01040                         // record the finishing time
01041                         // for the current vertex
01042                         tab [vertex] = ++ counter;
01043 
01044                         // return if this is the initial recursion level
01045                         if (s_vertex. empty ())
01046                                 return;
01047 
01048                         // restore the variables from the previous level
01049                         vertex = s_vertex. top ();
01050                         s_vertex. pop ();
01051                         edge = s_edge. top ();
01052                         s_edge. pop ();
01053                         maxedge = s_maxedge. top ();
01054                         s_maxedge. pop ();
01055                         continue;
01056                 }
01057 
01058                 // go to the deeper recursion level if possible
01059                 int_t next = edges [edge ++];
01060                 if (!tab [next])
01061                 {
01062                         // store the previous variables at the stacks
01063                         s_vertex. push (vertex);
01064                         s_edge. push (edge);
01065                         s_maxedge. push (maxedge);
01066 
01067                         // set the new vertex
01068                         vertex = next;
01069                         
01070                         // mark the new vertex as visited
01071                         tab [vertex] = -1;
01072                         
01073                         // determine the edges to be visited
01074                         edge = vertex ? edgeEnds [vertex - 1] :
01075                                 static_cast<int_t> (0);
01076                         maxedge = edgeEnds [vertex];
01077                 }
01078         }
01079 
01080         return;
01081 } /* diGraph::DFSfinishTimeStack */
01082 
01083 template <class wType> template <class Table>
01084 inline void diGraph<wType>::DFSfinishTime (Table &tab) const
01085 {
01086         // initialize the table and the counter
01087         for (int_t i = 0; i < nVertices; ++ i)
01088                 tab [i] = 0;
01089         int_t counter = 0;
01090 
01091         // compute the finishing time for each tree in the DFS forest
01092         for (int_t i = 0; i < nVertices; ++ i)
01093         {
01094                 if (!tab [i])
01095                         DFSfinishTimeStack (tab, i, counter);
01096         }
01097 
01098         #ifdef DIGRAPH_DEBUG
01099         int_t *tabdebug = new int_t [nVertices];
01100         for (int_t i = 0; i < nVertices; ++ i)
01101                 tabdebug [i] = 0;
01102         int_t counterdebug = 0;
01103         for (int_t i = 0; i < nVertices; ++ i)
01104                 if (!tabdebug [i])
01105                         DFSfinishTimeRecurrent (tabdebug, i, counterdebug);
01106         if (counter != counterdebug)
01107                 throw "DFSfinishTime: Wrong counter.";
01108         for (int_t i = 0; i < nVertices; ++ i)
01109         {
01110                 if (tab [i] != tabdebug [i])
01111                 {
01112                         sbug << "\nDFSfinishTime error: tabRec [" << i <<
01113                                 "] = " << tab [i] << ", tabStack [" << i <<
01114                                 "] = " << tabdebug [i] << ".\n";
01115                         throw "DFSfinishTime: Wrong numbers.";
01116                 }
01117         }
01118         sbug << "DEBUG: DFSfinishTime OK. ";
01119         #endif // DIGRAPH_DEBUG
01120         return;
01121 } /* diGraph::DFSfinishTime */
01122 
01123 // --------------------------------------------------
01124 
01125 template <class wType> template <class Table1, class Table2>
01126 inline bool diGraph<wType>::DFSforestRecurrent (Table1 &tab, Table1 &ntab,
01127         int_t vertex, int_t treeNumber, int_t countTrees,
01128         Table2 &compVertices, int_t &curVertex,
01129         diGraph<wType> *sccGraph, int_t *sccEdgeAdded) const
01130 {
01131         // add the vertex to the tree
01132         compVertices [curVertex ++] = vertex;
01133 
01134         // mark the vertex as belonging to the current tree
01135         tab [vertex] = treeNumber;
01136 //      if (sccGraph)
01137 //              ntab [treeNumber - 1] = countTrees;
01138 
01139         // build the tree recursively or record connections
01140         bool loop = false;
01141         for (int_t edge = vertex ? edgeEnds [vertex - 1] :
01142                 static_cast<int_t> (0);
01143                 edge < edgeEnds [vertex]; ++ edge)
01144         {
01145                 int_t next = edges [edge];
01146                 if (!tab [next])
01147                         loop |= DFSforestRecurrent (tab, ntab, next,
01148                                 treeNumber, countTrees, compVertices,
01149                                 curVertex, sccGraph, sccEdgeAdded);
01150                 else if (tab [next] == treeNumber)
01151                 {
01152                         if (sccGraph)
01153                         {
01154                                 int_t target = ntab [treeNumber - 1];
01155                                 if (sccEdgeAdded [target] != treeNumber)
01156                                 {
01157                                         sccGraph -> addEdge (target);
01158                                         sccEdgeAdded [target] = treeNumber;
01159                                 }
01160                         }
01161                         loop = true;
01162                 }
01163                 else if (sccGraph)
01164                 {
01165                         int_t target = ntab [tab [next] - 1];
01166                         if ((target >= 0) &&
01167                                 (sccEdgeAdded [target] != treeNumber))
01168                         {
01169                                 sccGraph -> addEdge (target);
01170                                 sccEdgeAdded [target] = treeNumber;
01171                         }
01172                 }
01173         }
01174 
01175         return loop;
01176 } /* diGraph::DFSforestRecurrent */
01177 
01178 template <class wType> template <class Table1, class Table2>
01179 inline bool diGraph<wType>::DFSforestStack (Table1 &tab, Table1 &ntab,
01180         int_t vertex, int_t treeNumber, int_t countTrees,
01181         Table2 &compVertices, int_t &curVertex,
01182         diGraph<wType> *sccGraph, int_t *sccEdgeAdded) const
01183 {
01184         // prepare stacks for the recursion
01185         std::stack<int_t> s_vertex;
01186         std::stack<int_t> s_edge;
01187         std::stack<int_t> s_maxedge;
01188         std::stack<bool> s_loop;
01189 
01190         // add the vertex to the tree
01191         compVertices [curVertex ++] = vertex;
01192 
01193         // mark the vertex as belonging to the current tree
01194         tab [vertex] = treeNumber;
01195 //      if (sccGraph)
01196 //              ntab [vertex] = countTrees;
01197 
01198         // build the tree recursively or record connections
01199         bool loop = false;
01200         int_t edge = vertex ? edgeEnds [vertex - 1] :
01201                 static_cast<int_t> (0);
01202         int_t maxedge = edgeEnds [vertex];
01203         while (1)
01204         {
01205                 // return to the previous recursion level
01206                 // if all the edges have been followed
01207                 if (edge >= maxedge)
01208                 {
01209                         // return if this is the initial recursion level
01210                         if (s_vertex. empty ())
01211                                 return loop;
01212 
01213                         // restore the variables from the previous level
01214                         vertex = s_vertex. top ();
01215                         s_vertex. pop ();
01216                         edge = s_edge. top ();
01217                         s_edge. pop ();
01218                         maxedge = s_maxedge. top ();
01219                         s_maxedge. pop ();
01220                         loop |= s_loop. top ();
01221                         s_loop. pop ();
01222                         continue;
01223                 }
01224 
01225                 // go to the deeper recursion level if possible
01226                 int_t next = edges [edge ++];
01227                 if (!tab [next])
01228                 {
01229                         // store the previous variables at the stacks
01230                         s_vertex. push (vertex);
01231                         s_edge. push (edge);
01232                         s_maxedge. push (maxedge);
01233                         s_loop. push (loop);
01234 
01235                         // set the new vertex
01236                         vertex = next;
01237                         
01238                         // add the vertex to the tree
01239                         compVertices [curVertex ++] = vertex;
01240 
01241                         // mark the vertex as belonging to the current tree
01242                         tab [vertex] = treeNumber;
01243 
01244                         // determine the edges to be visited
01245                         loop = false;
01246                         edge = vertex ? edgeEnds [vertex - 1] :
01247                                 static_cast<int_t> (0);
01248                         maxedge = edgeEnds [vertex];
01249                 }
01250                 else if (tab [next] == treeNumber)
01251                 {
01252                         if (sccGraph)
01253                         {
01254                                 int_t target = ntab [treeNumber - 1];
01255                                 if (sccEdgeAdded [target] != treeNumber)
01256                                 {
01257                                         sccGraph -> addEdge (target);
01258                                         sccEdgeAdded [target] = treeNumber;
01259                                 }
01260                         }
01261                         loop = true;
01262                 }
01263                 else if (sccGraph)
01264                 {
01265                         int_t target = ntab [tab [next] - 1];
01266                         if ((target >= 0) &&
01267                                 (sccEdgeAdded [target] != treeNumber))
01268                         {
01269                                 sccGraph -> addEdge (target);
01270                                 sccEdgeAdded [target] = treeNumber;
01271                         }
01272                 }
01273         }
01274 
01275         return loop;
01276 } /* diGraph::DFSforestStack */
01277 
01278 template <class wType> template <class Table1, class Table2, class Table3>
01279 inline int_t diGraph<wType>::DFSforest (const Table1 &ordered,
01280         Table2 &compVertices, Table3 &compEnds, bool nontrivial,
01281         diGraph<wType> *sccGraph) const
01282 {
01283         // prepare a table to record the numbers of DFS trees
01284         // to which the vertices belong (the tree numbers begin with 1)
01285         int_t *tab = new int_t [nVertices];
01286         for (int_t i = 0; i < nVertices; ++ i)
01287                 tab [i] = 0;
01288 
01289         // prepare a table to record the numbers of nontrivial trees
01290         // that correspond to the trees in 'tab' (these numbers begin with 0)
01291         int_t *ntab = new int_t [nVertices];
01292 
01293         // prepare a table to record the numbers of edges already in the
01294         // scc graph; "sccEdgeAdded [n] = m" indicates that the edge
01295         // m -> n has been added to the scc graph
01296         int_t *sccEdgeAdded = sccGraph ? new int_t [nVertices] :
01297                 static_cast<int_t *> (0);
01298         if (sccGraph)
01299         {
01300                 for (int_t n = 0; n < nVertices; ++ n)
01301                         sccEdgeAdded [n] = -1;
01302         }
01303 
01304         // prepare the official DFS tree number
01305         int_t treeNumber = 0;
01306 
01307         // prepare the data for keeping the nontrivial trees information
01308         int_t countTrees = 0;
01309         int_t curVertex = 0;
01310 
01311         // compute the DFS trees and connections between them
01312         for (int_t i = 0; i < nVertices; ++ i)
01313         {
01314                 // take the next vertex
01315                 int_t vertex = ordered [i];
01316 
01317                 // if the vertex already belongs to some tree, skip it
01318                 if (tab [vertex])
01319                         continue;
01320 
01321                 // add a vertex corresponding to the component
01322                 if (sccGraph)
01323                         sccGraph -> addVertex ();
01324 
01325                 // remember the previous vertex number
01326                 int_t prevVertex = curVertex;
01327 
01328                 // mark the entire component and record connections graph
01329                 if (sccGraph)
01330                         ntab [treeNumber] = countTrees;
01331                 ++ treeNumber;
01332                 bool loop = DFSforestStack (tab, ntab, vertex,
01333                         treeNumber, countTrees, compVertices,
01334                         curVertex, sccGraph, sccEdgeAdded);
01335 
01336                 // update the index bound for the vertex list
01337                 compEnds [countTrees ++] = curVertex;
01338 
01339                 // remove the component if it is trivial
01340                 if (nontrivial && !loop)
01341                 {
01342                         -- countTrees;
01343                         curVertex = prevVertex;
01344                         if (sccGraph)
01345                         {
01346                                 ntab [treeNumber - 1] = -1;
01347                                 sccGraph -> removeVertex ();
01348                         }
01349                 }
01350         }
01351 
01352         #ifdef DIGRAPH_DEBUG
01353         diGraph<wType> *sccGraphdebug = 0;
01354         if (sccGraph)
01355                 sccGraphdebug = new diGraph<wType>;
01356         // prepare a table to record the numbers of DFS trees
01357         // to which the vertices belong (the tree numbers begin with 1)
01358         int_t *tabdebug = new int_t [nVertices];
01359         for (int_t i = 0; i < nVertices; ++ i)
01360                 tabdebug [i] = 0;
01361 
01362         // prepare a table to record the numbers of nontrivial trees
01363         // to which the vertices belong (the tree numbers begin with 0)
01364         int_t *ntabdebug = new int_t [nVertices];
01365 
01366         // prepare a table to record the numbers of vertices from which
01367         // edges were added to the scc graph
01368         int_t *sccEdgeAddeddebug = sccGraph ? new int_t [nVertices] :
01369                 static_cast<int_t> (0);
01370         if (sccGraph)
01371         {
01372                 for (int_t n = 0; n < nVertices; ++ n)
01373                         sccEdgeAddeddebug [n] = -1;
01374         }
01375         // prepare the official DFS tree number
01376         int_t treeNumberdebug = 0;
01377 
01378         // prepare the data for keeping the nontrivial trees information
01379         int_t countTreesdebug = 0;
01380         int_t curVertexdebug = 0;
01381 
01382         int_t *compVerticesdebug = new int_t [nVertices];
01383         int_t *compEndsdebug = new int_t [nVertices];
01384         
01385         // compute the DFS trees and connections between them
01386         for (int_t i = 0; i < nVertices; ++ i)
01387         {
01388                 // take the next vertex
01389                 int_t vertex = ordered [i];
01390 
01391                 // if the vertex already belongs to some tree, skip it
01392                 if (tabdebug [vertex])
01393                         continue;
01394 
01395                 // add a vertex corresponding to the component
01396                 if (sccGraphdebug)
01397                         sccGraphdebug -> addVertex ();
01398 
01399                 // remember the previous vertex number
01400                 int_t prevVertex = curVertexdebug;
01401 
01402                 // mark the entire component and record connections graph
01403                 if (sccGraphdebug)
01404                         ntabdebug [treeNumberdebug] = countTreesdebug;
01405                 ++ treeNumberdebug;
01406                 bool loop = DFSforestRecurrent (tabdebug, ntabdebug, vertex,
01407                         treeNumberdebug, countTreesdebug, compVerticesdebug,
01408                         curVertexdebug, sccGraphdebug, sccEdgeAddeddebug);
01409 
01410                 // update the index bound for the vertex list
01411                 compEndsdebug [countTreesdebug ++] = curVertexdebug;
01412 
01413                 // remove the component if it is trivial
01414                 if (nontrivial && !loop)
01415                 {
01416                         -- countTreesdebug;
01417                         curVertexdebug = prevVertex;
01418                         if (sccGraphdebug)
01419                         {
01420                                 ntabdebug [treeNumberdebug - 1] = -1;
01421                                 sccGraphdebug -> removeVertex ();
01422                         }
01423                 }
01424         }
01425         if (countTrees != countTreesdebug)
01426                 throw "DFSforest: Wrong countTrees.";
01427         for (int_t i = 0; i < countTrees; ++ i)
01428                 if (compEnds [i] != compEndsdebug [i])
01429                         throw "DFSforest: Wrong compEnds.";
01430         for (int_t i = 0; i < compEndsdebug [countTrees - 1]; ++ i)
01431                 if (compVertices [i] != compVerticesdebug [i])
01432                         throw "DFSforest: Wrong vertices.";
01433         if (curVertex != curVertexdebug)
01434                 throw "DFSforest: Wrong curVertex.";
01435         for (int_t i = 0; i < nVertices; ++ i)
01436                 if (tab [i] != tabdebug [i])
01437                         throw "DFSforest: Wrong tab.";
01438         if (sccGraph)
01439         {
01440                 for (int_t i = 0; i < nVertices; ++ i)
01441                         if (ntab [i] != ntabdebug [i])
01442                                 throw "DFSforest: Wrong ntab.";
01443                 if (*sccGraph != *sccGraphdebug)
01444                         throw "DFSforest: Wrong graph.";
01445         }
01446         if (sccEdgeAdded)
01447         {
01448                 for (int_t i = 0; i < nVertices; ++ i)
01449                         if (sccEdgeAdded [i] != sccEdgeAddeddebug [i])
01450                                 throw "DFSforest: Wrong sccEdgeAdded.";
01451         }
01452         if (treeNumber != treeNumberdebug)
01453                 throw "DFSforest: Wrong treeNumber.";
01454         sbug << "DEBUG: DFSforest OK. ";
01455         if (!sccGraph)
01456                 sbug << "(Graphs not compared.) ";
01457         delete [] compVerticesdebug;
01458         delete [] compEndsdebug;
01459         if (sccGraphdebug)
01460                 delete sccGraphdebug;
01461         delete [] ntabdebug;
01462         delete [] tabdebug;
01463         if (sccEdgeAddeddebug)
01464                 delete [] sccEdgeAddeddebug;
01465         #endif // DIGRAPH_DEBUG
01466 
01467         if (sccEdgeAdded)
01468                 delete [] sccEdgeAdded;
01469         delete [] ntab;
01470         delete [] tab;
01471         return countTrees;
01472 } /* diGraph::DFSforest */
01473 
01474 template <class wType>
01475 inline int_t diGraph<wType>::shortestPath (int_t source, int_t destination)
01476         const
01477 {
01478         // make sure that the given vertex is present in the graph
01479         if ((source < 0) || (source >= nVertices) ||
01480                 (destination < 0) || (destination >= nVertices))
01481         {
01482                 throw "diGraph - shortest path: Wrong vertex number.";
01483         }
01484 
01485         // prepare an array of bits to store the information
01486         // on whether the given vertices have been visited or not
01487         BitField visited;
01488         visited. allocate (nVertices);
01489         visited. clearall (nVertices);
01490 
01491         // prepare queues for the BFS algorithm
01492         std::queue<int_t> q_vertex;
01493         std::queue<int_t> q_depth;
01494 
01495         // set the initial vertex
01496         int_t vertex = source;
01497         int_t depth = 0;
01498 
01499         while (1)
01500         {
01501                 // mark the current vertex as visited
01502                 visited. set (vertex);
01503 
01504                 // determine the depth of the vertices that will be analyzed
01505                 ++ depth;
01506 
01507                 // determine the edges to be checked
01508                 int_t firstedge = vertex ? edgeEnds [vertex - 1] :
01509                         static_cast<int_t> (0);
01510                 int_t maxedge = edgeEnds [vertex];
01511 
01512                 // put all the unvisited destination vertices on the stack
01513                 for (int_t edge = firstedge; edge < maxedge; ++ edge)
01514                 {
01515                         // determine the vertex pointed to by this edge
01516                         int_t next = edges [edge];
01517 
01518                         // if this is the destination vertex then return
01519                         // the shortest path length; note: this vertex
01520                         // might be visited if checking a loop, so it is
01521                         // important to check the destination first
01522                         if (next == destination)
01523                         {
01524                                 visited. free ();
01525                                 return depth;
01526                         }
01527 
01528                         // if the vertex was already visited then skip it
01529                         if (visited. test (next))
01530                                 continue;
01531 
01532                         // add the vertex to the queue
01533                         q_vertex. push (next);
01534                         q_depth. push (depth);
01535                 }
01536 
01537                 // if there are no vertices whose neighbors are to be
01538                 // analyzed and the destination vertex was not found
01539                 // then there is no path to that vertex
01540                 if (q_vertex. empty ())
01541                 {
01542                         visited. free ();
01543                         return 0;
01544                 }
01545 
01546                 // pick up a vertex stored in the queue
01547                 vertex = q_vertex. front ();
01548                 q_vertex. pop ();
01549                 depth = q_depth. front ();
01550                 q_depth. pop ();
01551         }
01552 } /* diGraph::shortestPath */
01553 
01554 template <class wType>
01555 inline int_t diGraph<wType>::shortestLoop (int_t origin) const
01556 {
01557         return shortestPath (origin, origin);
01558 } /* diGraph::shortestLoop */
01559 
01560 template <class wType>
01561 template <class lenTable, class weightsType, class roundType>
01562 inline void diGraph<wType>::Dijkstra (const roundType &rounding,
01563         int_t source, lenTable &len, weightsType &edgeWeights) const
01564 {
01565         // use the Fibonacci heap as a priority queue
01566         FibonacciHeap<posWeight> Q (nVertices);
01567 
01568         // add the vertices to the heap with the initial shortest path
01569         // lengths: 0 to the source, plus infinity to all the others
01570         for (int_t v = 0; v < nVertices; ++ v)
01571         {
01572                 posWeight w (0);
01573                 if (v != source)
01574                         w. setInfinity ();
01575                 int_t number = Q. Insert (w);
01576                 if (number != v)
01577                 {
01578                         throw "Wrong implementation of Fibonacci heap "
01579                                 "for this version of Dijkstra.";
01580                 }
01581         }
01582 
01583         // pick up vertices from the priority queue, record the length
01584         // of the shortest path to them, and modify the remaining paths
01585         for (int_t i = 0; i < nVertices; ++ i)
01586         {
01587                 // extract the minimal vertex from the queue
01588                 int_t minVertex = Q. Minimum ();
01589                 posWeight minWeight = Q. ExtractMinimum ();
01590 
01591                 if (minWeight. isInfinity ())
01592                 {
01593                         len [minVertex] = -1;
01594                         continue;
01595                 }
01596                 wType minValue = minWeight. getValue ();
01597                 len [minVertex] = minValue;
01598 
01599                 // go through all the edges emanating from this vertex
01600                 // and update the path lengths for the target vertices
01601                 int_t edge = minVertex ? edgeEnds [minVertex - 1] :
01602                         static_cast<int_t> (0);
01603                 int_t maxEdge = edgeEnds [minVertex];
01604                 for (; edge < maxEdge; ++ edge)
01605                 {
01606                         // determine the vertex at the other end of the edge
01607                         int_t nextVertex = edges [edge];
01608 
01609                         // if the path that runs through the extracted
01610                         // vertex is shorter, then make a correction
01611                         const posWeight &nextWeight = Q. Value (nextVertex);
01612                         wType newWeight = rounding. add_down (minValue,
01613                                 edgeWeights [edge]);
01614                         if (newWeight < 0)
01615                                 newWeight = 0;
01616                         if (nextWeight. isInfinity () ||
01617                                 (newWeight < nextWeight. getValue ()))
01618                         {
01619                                 Q. DecreaseKey (nextVertex,
01620                                         posWeight (newWeight));
01621                         }
01622                 }
01623         }
01624         return;
01625 } /* diGraph::Dijkstra */
01626 
01627 template <class wType>
01628 template <class lenTable, class roundType>
01629 inline void diGraph<wType>::Dijkstra (const roundType &rounding,
01630         int_t source, lenTable &len) const
01631 {
01632         this -> Dijkstra (rounding, source, len, this -> weights);
01633         return;
01634 } /* diGraph::Dijkstra */
01635 
01636 template <class wType>
01637 template <class lenTable>
01638 inline void diGraph<wType>::Dijkstra (int_t source, lenTable &len) const
01639 {
01640         const dummyRounding<wType> rounding = dummyRounding<wType> ();
01641         this -> Dijkstra (rounding, source, len);
01642         return;
01643 } /* diGraph::Dijkstra */
01644 
01645 template <class wType> template <class lenTable, class predTable,
01646         class roundType>
01647 inline bool diGraph<wType>::BellmanFord (const roundType &rounding,
01648         int_t source, lenTable &len, wType *infinity, predTable pred) const
01649 {
01650         // make sure the source vertex number is correct
01651         if ((source < 0) || (source >= nVertices))
01652                 throw "Bellman-Ford: Wrong source vertex number.";
01653 
01654         // prepare marks to indicate finite values (not "infinity")
01655         BitField finite;
01656         finite. allocate (nVertices);
01657         finite. clearall (nVertices);
01658         finite. set (source);
01659         len [source] = 0;
01660 
01661         // count the negative vertices
01662         int_t countNegative = 0;
01663 
01664         // set the initial predecessors
01665         for (int_t i = 0; i < nVertices; ++ i)
01666                 pred [i] = -1;
01667 
01668         // update the lenghts of the paths repeatedly (max nVertices times)
01669         bool noNegativeLoop = false;
01670         int_t counter = 0;
01671         for (; counter <= nVertices; ++ counter)
01672         {
01673                 bool modified = false;
01674                 int_t curEdge = 0;
01675                 for (int_t vertex = 0; vertex < nVertices; ++ vertex)
01676                 {
01677                         int_t maxEdge = edgeEnds [vertex];
01678                         if (!finite. test (vertex))
01679                         {
01680                                 curEdge = maxEdge;
01681                                 continue;
01682                         }
01683                         for (; curEdge < maxEdge; ++ curEdge)
01684                         {
01685                                 int_t next = edges [curEdge];
01686                                 wType newlen = rounding. add_down
01687                                         (len [vertex], weights [curEdge]);
01688                                 if (!finite. test (next))
01689                                 {
01690                                         finite. set (next);
01691                                         modified = true;
01692                                         len [next] = newlen;
01693                                         pred [next] = vertex;
01694                                         if (newlen < 0)
01695                                                 ++ countNegative;
01696                                 }
01697                                 else if (newlen < len [next])
01698                                 {
01699                                         modified = true;
01700                                         if (!(len [next] < 0) &&
01701                                                 (newlen < 0))
01702                                         {
01703                                                 ++ countNegative;
01704                                         }
01705                                         len [next] = newlen;
01706                                         pred [next] = vertex;
01707                                 }
01708                         }
01709                 }
01710                 if (countNegative == nVertices)
01711                 {
01712                         noNegativeLoop = false;
01713                         ++ counter;
01714                         break;
01715                 }
01716                 if (!modified)
01717                 {
01718                         noNegativeLoop = true;
01719                         ++ counter;
01720                         break;
01721                 }
01722         }
01723 
01724         // show a message on how many loops have been done
01725         if (false && chomp::homology::sbug. show)
01726         {
01727                 chomp::homology::sbug << "Bellman-Ford: " <<
01728                         counter << ((counter > 1) ? " loops (" : " loop (") <<
01729                         nVertices << " vertices, " << countNegative <<
01730                         " negative). " <<
01731                         (noNegativeLoop ? "No negative loops.\n" :
01732                         "A negative loop found.\n");
01733         }
01734 
01735         // compute the value for the infinity and set the undefined distances
01736         if (infinity && noNegativeLoop)
01737         {
01738                 wType infty (0);
01739                 bool first = true;
01740                 for (int_t i = 0; i < nVertices; ++ i)
01741                 {
01742                         if (!finite. test (i))
01743                                 continue;
01744                         if (first)
01745                         {
01746                                 infty = len [i];
01747                                 first = false;
01748                         }
01749                         else if (infty < len [i])
01750                         {
01751                                 infty = len [i];
01752                         }
01753                 }
01754                 infty = infty + 1;
01755                 for (int_t i = 0; i < nVertices; ++ i)
01756                 {
01757                         if (!finite. test (i))
01758                                 len [i] = infty;
01759                 }
01760                 *infinity = infty;
01761         }
01762 
01763         finite. free ();
01764         return noNegativeLoop;
01765 } /* diGraph::BellmanFord */
01766 
01767 template <class wType> template <class lenTable, class predTable>
01768 inline bool diGraph<wType>::BellmanFord (int_t source, lenTable &len,
01769         wType *infinity, predTable pred) const
01770 {
01771         const dummyRounding<wType> rounding = dummyRounding<wType> ();
01772         return this -> BellmanFord (rounding, source, len, infinity, pred);
01773 } /* diGraph::BellmanFord */
01774 
01775 template <class wType> template <class roundType>
01776 inline bool diGraph<wType>::BellmanFord (const roundType &rounding,
01777         int_t source) const
01778 {
01779         chomp::homology::auto_array<wType> len_ptr (new wType [nVertices]);
01780         wType *len = len_ptr. get ();
01781         wType *infinity = 0;
01782         dummyArray tab;
01783         return BellmanFord (rounding, source, len, infinity, tab);
01784 } /* diGraph::BellmanFord */
01785 
01786 template <class wType>
01787 inline bool diGraph<wType>::BellmanFord (int_t source) const
01788 {
01789         const dummyRounding<wType> rounding = dummyRounding<wType> ();
01790         return BellmanFord (rounding, source);
01791 } /* diGraph::BellmanFord */
01792 
01793 template <class wType>
01794 inline wType diGraph<wType>::Edmonds () const
01795 {
01796         // create a list of edges with weights and sort this list
01797         std::vector<edgeTriple> edgeTriples (countEdges ());
01798         int_t prevEdge = 0;
01799         int_t curEdge = 0;
01800         for (int_t vert = 0; vert < nVertices; ++ vert)
01801         {
01802                 while (curEdge < edgeEnds [vert])
01803                 {
01804                         edgeTriple &e = edgeTriples [curEdge];
01805                         e. vert1 = vert;
01806                         e. vert2 = edges [curEdge];
01807                         e. weight = weights [curEdge];
01808                         ++ curEdge;
01809                 }
01810                 prevEdge = curEdge;
01811         }
01812         std::sort (edgeTriples. begin (), edgeTriples. end ());
01813 
01814         // create a forest which initially contains single vertices
01815         chomp::homology::auto_array<int_t> root_ptr (new int_t [nVertices]);
01816         int_t *root = root_ptr. get ();
01817         for (int_t vert = 0; vert < nVertices; ++ vert)
01818         {
01819                 root [vert] = -1;
01820         }
01821 
01822         // go through the edges and join the trees, but avoid loops
01823         wType totalWeight = 0;
01824         int_t nEdges = countEdges ();
01825         for (int_t curEdge = 0; curEdge < nEdges; ++ curEdge)
01826         {
01827                 // determine the roots of both vertices of the edge
01828                 // and compress the paths
01829                 edgeTriple &e = edgeTriples [curEdge];
01830                 int_t root1 = e. vert1;
01831                 while (root [root1] >= 0)
01832                 {
01833                         root1 = root [root1];
01834                 }
01835                 int_t vert1 = e. vert1;
01836                 while (root [vert1] >= 0)
01837                 {
01838                         int_t next = root [vert1];
01839                         root [vert1] = root1;
01840                         vert1 = next;
01841                 }
01842                 int_t root2 = e. vert2;
01843                 while (root [root2] >= 0)
01844                 {
01845                         root2 = root [root2];
01846                 }
01847                 int_t vert2 = e. vert2;
01848                 while (root [vert2] >= 0)
01849                 {
01850                         int_t next = root [vert2];
01851                         root [vert2] = root2;
01852                         vert2 = next;
01853                 }
01854 
01855                 // skip the edge if it closes a loop
01856                 if (root1 == root2)
01857                         continue;
01858 
01859                 // add the weight
01860                 totalWeight += e. weight;
01861 
01862                 // join the trees
01863                 root [root1] = root2;
01864         }
01865         return totalWeight;
01866 } /* diGraph::Edmonds */
01867 
01868 template <class wType>
01869 inline wType diGraph<wType>::EdmondsOld () const
01870 {
01871         // create a list of edges with weights and sort this list
01872         std::vector<edgeTriple> edgeTriples (countEdges ());
01873         int_t prevEdge = 0;
01874         int_t curEdge = 0;
01875         for (int_t vert = 0; vert < nVertices; ++ vert)
01876         {
01877                 while (curEdge < edgeEnds [vert])
01878                 {
01879                         edgeTriple &e = edgeTriples [curEdge];
01880                         e. vert1 = vert;
01881                         e. vert2 = edges [curEdge];
01882                         e. weight = weights [curEdge];
01883                         ++ curEdge;
01884                 }
01885                 prevEdge = curEdge;
01886         }
01887         std::sort (edgeTriples. begin (), edgeTriples. end ());
01888 
01889         // create a forest which initially contains single vertices
01890         chomp::homology::auto_array<int_t> forest_ptr
01891                 (new int_t [nVertices]);
01892         int_t *forest = forest_ptr. get ();
01893         chomp::homology::auto_array<int_t> next_ptr
01894                 (new int_t [nVertices]);
01895         int_t *next = next_ptr. get ();
01896         chomp::homology::auto_array<int_t> prev_ptr
01897                 (new int_t [nVertices]);
01898         int_t *prev = prev_ptr. get ();
01899         for (int_t vert = 0; vert < nVertices; ++ vert)
01900         {
01901                 forest [vert] = vert;
01902                 next [vert] = -1;
01903                 prev [vert] = -1;
01904         }
01905 
01906         // go through the edges and join the trees, but avoid loops
01907         wType totalWeight = 0;
01908         int_t nEdges = countEdges ();
01909         for (int_t curEdge = 0; curEdge < nEdges; ++ curEdge)
01910         {
01911                 // check the edge and skip it if it closes a loop
01912                 edgeTriple &e = edgeTriples [curEdge];
01913                 if (forest [e. vert1] == forest [e. vert2])
01914                         continue;
01915 
01916                 // add the weight
01917                 totalWeight += e. weight;
01918 
01919                 // go to the end of the first tree
01920                 int_t vert1 = e. vert1;
01921                 while (next [vert1] >= 0)
01922                 {
01923                         vert1 = next [vert1];
01924                 }
01925                 
01926                 // go to the beginning of the second tree
01927                 int_t vert2 = e. vert2;
01928                 while (prev [vert2] >= 0)
01929                 {
01930                         vert2 = prev [vert2];
01931                 }
01932 
01933                 // join the trees and modify the numbers
01934                 next [vert1] = vert2;
01935                 prev [vert2] = vert1;
01936                 int_t treeNumber = forest [vert1];
01937                 while (vert2 >= 0)
01938                 {
01939                         forest [vert2] = treeNumber;
01940                         vert2 = next [vert2];
01941                 }
01942         }
01943         return totalWeight;
01944 } /* diGraph::EdmondsOld */
01945 
01946 template <class wType>
01947 template <class arrayType, class roundType>
01948 inline wType diGraph<wType>::FloydWarshall (const roundType &rounding,
01949         arrayType &arr, bool setInfinity, bool ignoreNegLoop) const
01950 {
01951         // do nothing if the graph is empty
01952         if (!nVertices)
01953                 return 0;
01954 
01955         // prepare marks to indicate finite values (not "infinity")
01956         BitField *finite = new BitField [nVertices];
01957         for (int_t i = 0; i < nVertices; ++ i)
01958         {
01959                 finite [i]. allocate (nVertices);
01960                 finite [i]. clearall (nVertices);
01961         }
01962 
01963         // create the initial values of the array based on the edge weights
01964         int_t curEdge = 0;
01965         for (int_t source = 0; source < nVertices; ++ source)
01966         {
01967                 bool diagset = false;
01968                 while (curEdge < edgeEnds [source])
01969                 {
01970                         int_t target = edges [curEdge];
01971                         const wType &weight = weights [curEdge];
01972                         if (source == target)
01973                         {
01974                                 if (weight < 0)
01975                                 {
01976                                         arr [source] [target] = weight;
01977                                         diagset = true;
01978                                 }
01979                         }
01980                         else
01981                         {
01982                                 arr [source] [target] = weight;
01983                                 finite [source]. set (target);
01984                         }
01985                         ++ curEdge;
01986                 }
01987                 if (!diagset)
01988                         arr [source] [source] = 0;
01989                 finite [source]. set (source);
01990         }
01991 
01992         // find the shortest paths between the vertices (dynamic programming)
01993         for (int_t k = 0; k < nVertices; ++ k)
01994         {
01995                 for (int_t i = 0; i < nVertices; ++ i)
01996                 {
01997                         if (!finite [i]. test (k))
01998                                 continue;
01999                         for (int_t j = 0; j < nVertices; ++ j)
02000                         {
02001                                 if (!finite [k]. test (j))
02002                                         continue;
02003                                 const wType sum = rounding. add_down
02004                                         (arr [i] [k], arr [k] [j]);
02005                                 if (finite [i]. test (j))
02006                                 {
02007                                         if (sum < arr [i] [j])
02008                                                 arr [i] [j] = sum;
02009                                 }
02010                                 else
02011                                 {
02012                                         arr [i] [j] = sum;
02013                                         finite [i]. set (j);
02014                                 }
02015                         }
02016                 }
02017         }
02018 
02019         // verify if a negative loop exists by checking for a negative
02020         // result in the diagonal
02021         if (!ignoreNegLoop)
02022         {
02023                 for (int_t i = 0; i < nVertices; ++ i)
02024                 {
02025                         if (arr [i] [i] < 0)
02026                                 throw "Negative loop in Floyd-Warshall.";
02027                 }
02028         }
02029 
02030         // prepare a variable to store the returned result
02031         wType result = 0;
02032 
02033         // compute the value for the infinity and fill in the array
02034         // if requested to do so
02035         if (setInfinity)
02036         {
02037                 wType &infinity = result;
02038                 for (int_t i = 0; i < nVertices; ++ i)
02039                 {
02040                         for (int_t j = 0; j < nVertices; ++ j)
02041                         {
02042                                 if (finite [i]. test (j) &&
02043                                         (infinity <= arr [i] [j]))
02044                                 {
02045                                         infinity = rounding. add_up
02046                                                 (arr [i] [j], 1);
02047                                 }
02048                         }
02049                 }
02050                 for (int_t i = 0; i < nVertices; ++ i)
02051                 {
02052                         for (int_t j = 0; j < nVertices; ++ j)
02053                         {
02054                                 if (!finite [i]. test (j))
02055                                         arr [i] [j] = infinity;
02056                         }
02057                 }
02058         }
02059 
02060         // otherwise, only compute the minimum path weight
02061         else
02062         {
02063                 wType &minWeight = result;
02064                 bool firstTime = true;
02065                 for (int_t i = 0; i < nVertices; ++ i)
02066                 {
02067                         for (int_t j = 0; j < nVertices; ++ j)
02068                         {
02069                                 if (finite [i]. test (j))
02070                                 {
02071                                         if (firstTime)
02072                                         {
02073                                                 minWeight = arr [i] [j];
02074                                                 firstTime = false;
02075                                         }
02076                                         else if (arr [i] [j] < minWeight)
02077                                         {
02078                                                 minWeight = arr [i] [j];
02079                                         }
02080                                 }
02081                         }
02082                 }
02083         }
02084 
02085         // release the 'finite' bitfields
02086         for (int_t i = 0; i < nVertices; ++ i)
02087                 finite [i]. free ();
02088         delete [] finite;
02089 
02090         return result;
02091 } /* diGraph::FloydWarshall */
02092 
02093 template <class wType>
02094 template <class arrayType>
02095 inline wType diGraph<wType>::FloydWarshall (arrayType &arr,
02096         bool setInfinity, bool ignoreNegLoop) const
02097 {
02098         const dummyRounding<wType> rounding = dummyRounding<wType> ();
02099         return FloydWarshall (rounding, arr, setInfinity, ignoreNegLoop);
02100 } /* diGraph::FloydWarshall */
02101 
02102 template <class wType>
02103 template <class arrayType, class roundType>
02104 inline wType diGraph<wType>::Johnson (const roundType &rounding,
02105         arrayType &arr, bool setInfinity, bool ignoreNegLoop) const
02106 {
02107         // DEBUG VERIFICATION
02108         if (false && sbug. show)
02109         {
02110                 timeused stopWatch;
02111                 wType res = FloydWarshall (rounding,
02112                         arr, setInfinity, ignoreNegLoop);
02113                 chomp::homology::sbug << "\n[Floyd-Warshall: " << res <<
02114                         ", " << (double) stopWatch << " sec]\n";
02115         }
02116         // debug time measurement (see below)
02117 //      timeused stopWatch;
02118 
02119         // do nothing if the graph is empty
02120         if (!nVertices)
02121                 return 0;
02122 
02123         // STEP 1: Compute the shortest paths to any vertex from an
02124         // artificial extra vertex connected to every other vertex in the
02125         // graph by an edge of weight zero. Use Bellman-Ford for this.
02126         wType *len = new wType [nVertices];
02127         for (int_t i = 0; i < nVertices; ++ i)
02128                 len [i] = 0;
02129 
02130         // update the lenghts of the paths repeatedly (max nVertices times)
02131         bool noNegativeLoop = false;
02132         int_t counter = 0;
02133         for (; counter <= nVertices; ++ counter)
02134         {
02135                 bool modified = false;
02136                 int_t curEdge = 0;
02137                 for (int_t vertex = 0; vertex < nVertices; ++ vertex)
02138                 {
02139                         int_t maxEdge = edgeEnds [vertex];
02140                         for (; curEdge < maxEdge; ++ curEdge)
02141                         {
02142                                 int_t next = edges [curEdge];
02143                                 wType newlen = rounding. add_down
02144                                         (len [vertex], weights [curEdge]);
02145                                 if (newlen < len [next])
02146                                 {
02147                                         // this "if" statement is necessary
02148                                         // because of a bug in GCC 3.4.2...
02149                                         if (counter > nVertices)
02150                                         {
02151                                                 std::cout << vertex;
02152                                         }
02153                                         modified = true;
02154                                         len [next] = newlen;
02155                                 }
02156                         }
02157                 }
02158                 if (!modified)
02159                 {
02160                         noNegativeLoop = true;
02161                         ++ counter;
02162                         break;
02163                 }
02164         }
02165         if (!ignoreNegLoop && !noNegativeLoop)
02166                 throw "Negative loop found in Johnson's algorithm.";
02167 
02168         // STEP 2: Re-weigh the edges using the computed path lengths.
02169         wType *newWeights = new wType [edgeEnds [nVertices - 1]];
02170         int_t edge = 0;
02171         for (int_t vertex = 0; vertex < nVertices; ++ vertex)
02172         {
02173                 int_t maxEdge = edgeEnds [vertex];
02174                 for (; edge < maxEdge; ++ edge)
02175                 {
02176                         wType w = weights [edge];
02177                         w = rounding. add_down (w, len [vertex]);
02178                         w = rounding. sub_down (w, len [edges [edge]]);
02179                         newWeights [edge] = (w < 0) ?
02180                                 static_cast<wType> (0) : w;
02181                 }
02182         }
02183 
02184         // STEP 3: Run the Dijkstra algorithm for every vertex to compute
02185         // the shortest paths to all the other vertices.
02186         // Negative entries indicate no path existence.
02187         for (int_t u = 0; u < nVertices; ++ u)
02188         {
02189                 this -> Dijkstra (rounding, u, arr [u], newWeights);
02190         }
02191         delete [] newWeights;
02192 
02193         // STEP 4: Make a correction to the computed path lengths.
02194         // Compute the value for infinity if requested to.
02195         // Otherwise compute the minimum of path lengths.
02196         wType result = 0;
02197         if (setInfinity)
02198         {
02199                 wType &infinity = result;
02200                 for (int_t u = 0; u < nVertices; ++ u)
02201                 {
02202                         for (int_t v = 0; v < nVertices; ++ v)
02203                         {
02204                                 wType w = arr [u] [v];
02205                                 if (w < 0)
02206                                         continue;
02207                                 w = rounding. add_down (w, len [v]);
02208                                 w = rounding. sub_down (w, len [u]);
02209                                 if (w < infinity)
02210                                         continue;
02211                                 infinity = rounding. add_up (w, 1);
02212                         }
02213                 }
02214                 for (int_t u = 0; u < nVertices; ++ u)
02215                 {
02216                         for (int_t v = 0; v < nVertices; ++ v)
02217                         {
02218                                 wType w = arr [u] [v];
02219                                 if (w < 0)
02220                                 {
02221                                         arr [u] [v] = infinity;
02222                                         continue;
02223                                 }
02224                                 w = rounding. add_down (w, len [v]);
02225                                 arr [u] [v] =
02226                                         rounding. sub_down (w, len [u]);
02227                         }
02228                 }
02229         }
02230         else
02231         {
02232                 wType &minWeight = result;
02233                 bool firstTime = true;
02234                 for (int_t u = 0; u < nVertices; ++ u)
02235                 {
02236                         for (int_t v = 0; v < nVertices; ++ v)
02237                         {
02238                                 wType w = arr [u] [v];
02239                                 if (w < 0)
02240                                         continue;
02241                                 w = rounding. add_down (w, len [v]);
02242                                 w = rounding. sub_down (w, len [u]);
02243                                 if (firstTime)
02244                                 {
02245                                         minWeight = w;
02246                                         firstTime = false;
02247                                 }
02248                                 else if (w < minWeight)
02249                                         minWeight = w;
02250                         }
02251                 }
02252         }
02253         delete [] len;
02254 
02255         // DEBUG VERIFICATION
02256         if (false && sbug. show)
02257         {
02258 //              chomp::homology::sbug << "[Johnson: " << result <<
02259 //                      ", " << (double) stopWatch << " sec]\n";
02260         }
02261 
02262         return result;
02263 } /* diGraph::Johnson */
02264 
02265 template <class wType>
02266 template <class arrayType>
02267 inline wType diGraph<wType>::Johnson (arrayType &arr,
02268         bool setInfinity, bool ignoreNegLoop) const
02269 {
02270         const dummyRounding<wType> rounding = dummyRounding<wType> ();
02271         return Johnson (rounding, arr, setInfinity, ignoreNegLoop);
02272 } /* diGraph::Johnson */
02273 
02274 template <class wType>
02275 template <class roundType>
02276 wType diGraph<wType>::minPathWeight (const roundType &rounding,
02277         bool ignoreNegLoop, int sparseGraph) const
02278 {
02279         // if the graph is empty, return 0 as the minimum path weight
02280         if (nVertices <= 0)
02281                 return 0;
02282 
02283         // allocate memory for an array of arrays to store min path weights
02284         wType **arr = new wType * [nVertices];
02285         for (int_t i = 0; i < nVertices; ++ i)
02286                 arr [i] = new wType [nVertices];
02287 
02288         // determine whether to run the Floyd-Warshall algorithm
02289         // or Johnson's algorithm
02290         bool sparse = false;
02291         if (sparseGraph == 1)
02292                 sparse = true;
02293         else if (sparseGraph != 0)
02294         {
02295                 double nEdgesD = this -> countEdges ();
02296                 double nVerticesD = this -> countVertices ();
02297                 if ((nVerticesD > 3000) && (nEdgesD < nVerticesD * 1000) &&
02298                         (nEdgesD < nVerticesD * nVerticesD / 50))
02299                 {
02300                         sparse = true;
02301                 }
02302         }
02303 
02304         // run the Johnson's or Floyd-Warshall algorithm
02305         wType minWeight = sparse ?
02306                 this -> Johnson (rounding, arr, false, ignoreNegLoop) :
02307                 this -> FloydWarshall (rounding, arr, false, ignoreNegLoop);
02308 
02309 /*      // compute the minimum of all the paths
02310         wType minWeight = arr [0] [0];
02311         for (int_t i = 0; i < nVertices; ++ i)
02312         {
02313                 for (int_t j = 0; j < nVertices; ++ j)
02314                 {
02315                         const wType &weight = arr [i] [j];
02316                         if (weight < minWeight)
02317                                 minWeight = weight;
02318                 }
02319         }
02320 */
02321         // release the memory
02322         for (int_t i = 0; i < nVertices; ++ i)
02323                 delete [] (arr [i]);
02324         delete [] arr;
02325 
02326         return minWeight;
02327 } /* diGraph::minPathWeight */
02328 
02329 template <class wType>
02330 wType diGraph<wType>::minPathWeight (bool ignoreNegLoop, int sparseGraph)
02331         const
02332 {
02333         const dummyRounding<wType> rounding = dummyRounding<wType> ();
02334         return this -> minPathWeight (rounding, ignoreNegLoop, sparseGraph);
02335 } /* diGraph::minPathWeight */
02336 
02337 template <class wType> template <class outType>
02338 inline outType &diGraph<wType>::show (outType &out, bool showWeights) const
02339 {
02340         out << "; Directed graph: " << nVertices << " vertices.\n";
02341         int_t curEdge = 0;
02342         for (int_t i = 0; i < nVertices; ++ i)
02343         {
02344                 for (; curEdge < edgeEnds [i]; ++ curEdge)
02345                 {
02346                         out << i << " -> " << edges [curEdge];
02347                         if (showWeights)
02348                                 out << " [" << weights [curEdge] << "]\n";
02349                         else
02350                                 out << "\n";
02351                 }
02352         }
02353         out << "; This is the end of the graph.\n";
02354         return out;
02355 } /* diGraph::show */
02356 
02357         
02358 // --------------------------------------------------
02359 
02360 /// Writes a graph in the text mode to an output stream.
02361 template <class wType>
02362 inline std::ostream &operator << (std::ostream &out, const diGraph<wType> &g)
02363 {
02364         return g. show (out, false);
02365 } /* operator << */
02366 
02367 // --------------------------------------------------
02368 
02369 /// Computes strongly connected components of the graph 'g'. Creates
02370 /// the graph 'scc' in which each vertex corresponds to one component.
02371 /// The graph 'scc' given as an argument must be initially empty.
02372 /// The table 'compVertices' is filled with the numbers of vertices in 'g'
02373 /// which form the components, and the indices that end the listing
02374 /// for each component are stored in the table 'compEnds'.
02375 /// Returns the number of strongly connected components found.
02376 template <class wType, class Table1, class Table2>
02377 inline int_t SCC (const diGraph<wType> &g, Table1 &compVertices,
02378         Table2 &compEnds, diGraph<wType> *scc = 0,
02379         diGraph<wType> *transposed = 0, bool copyweights = false)
02380 {
02381         // prepare two tables
02382         int_t nVert = g. countVertices ();
02383         int_t *ordered = new int_t [nVert];
02384         int_t *tab = new int_t [nVert];
02385 
02386         // compute the list of vertices in the descending finishing time
02387         g. DFSfinishTime (tab);
02388         for (int_t i = 0; i < nVert; ++ i)
02389                 ordered [nVert - tab [i]] = i;
02390         delete [] tab;
02391 
02392         // create the transposed graph
02393         diGraph<wType> gT;
02394         if (!transposed)
02395                 transposed = &gT;
02396         g. transpose (*transposed, copyweights);
02397 
02398         // extract the DFS forest of gT in the given order of vertices
02399         int_t n = transposed -> DFSforest (ordered, compVertices, compEnds,
02400                 true, scc);
02401 
02402         // cleanup memory and return the number of components
02403         delete [] ordered;
02404         return n;
02405 } /* SCC */
02406 
02407 // --------------------------------------------------
02408 
02409 /// Computes strongly connected components of the graph 'g'
02410 /// using Tarjan's algorithm (as described in the Wikipedia).
02411 /// Tha advantage of this approach over the one described in Cormen's
02412 /// textbook is that the transposed graph need not be computed.
02413 /// However, this algorithm might be slightly slower than the other one.
02414 /// The table 'compVertices' is filled with the numbers of vertices in 'g'
02415 /// which form the components, and the indices that end the listing
02416 /// for each component are stored in the table 'compEnds'.
02417 /// Returns the number of strongly connected components found.
02418 template <class wType, class Table1, class Table2>
02419 inline int_t SCC_Tarjan (const diGraph<wType> &g, Table1 &compVertices,
02420         Table2 &compEnds)
02421 {
02422         // return the obvious result if the graph is empty
02423         int_t nVertices = g. countVertices ();
02424         if (!nVertices)
02425                 return 0;
02426 
02427         // prepare an array of discovery times for all the vertices
02428         // (zero == not yet discovered)
02429         std::vector<int_t> dfsIndex (nVertices, 0);
02430 
02431         // prepare an array of the minimal index of a node reachable
02432         // from each of the vertices
02433         std::vector<int_t> lowLink (nVertices, 0);
02434 
02435         // prepare an empty stack of nodes
02436         std::stack<int_t> s_nodes;
02437 
02438         // prepare an array of bits indicating whether the vertices are
02439         // in the stack of nodes or not
02440         BitField inTheStack;
02441         inTheStack. allocate (nVertices);
02442         inTheStack. clearall (nVertices);
02443 
02444         // prepare the number of strongly connected components
02445         int_t nComponents = 0;
02446 
02447         // prepare the current position in the array 'compVertices'
02448         int_t posVertices = 0;
02449 
02450         // remember the next vertex number in the graph to scan
02451         // whether this vertex has already been visited or not
02452         int_t vertexToScan = 0;
02453 
02454         // prepare a variable for storing the discovery time in the DFS
02455         int_t discoveryTime = 0;
02456 
02457         // prepare stacks for the DFS recursion
02458         std::stack<int_t> s_vertex;
02459         std::stack<int_t> s_edge;
02460         std::stack<int_t> s_maxedge;
02461 
02462         // initialize the number of the currently processed vertex
02463         int_t vertex = -1;
02464 
02465         // initialize the range of edges to be visited
02466         int_t edge = 0;
02467         int_t maxedge = 0;
02468 
02469         while (1)
02470         {
02471                 // return to the previous recursion level
02472                 // if all the edges have been checked
02473                 if (edge >= maxedge)
02474                 {
02475                         // extract a strongly connected component
02476                         if ((vertex >= 0) && (lowLink [vertex] ==
02477                                 dfsIndex [vertex]))
02478                         {
02479                                 int_t v = 0;
02480                                 do
02481                                 {
02482                                         v = s_nodes. top ();
02483                                         s_nodes. pop ();
02484                                         inTheStack. clear (v);
02485                                         compVertices [posVertices ++] = v;
02486                                 } while (v != vertex);
02487                                 compEnds [nComponents ++] = posVertices;
02488                         }
02489 
02490                         // if this is the top level of the recursion
02491                         // then find another unvisited vertex or return
02492                         if (s_vertex. empty ())
02493                         {
02494                                 // find an unvisited vertex in the graph
02495                                 while ((vertexToScan < nVertices) &&
02496                                         (dfsIndex [vertexToScan] != 0))
02497                                 {
02498                                         ++ vertexToScan;
02499                                 }
02500 
02501                                 // return the result if all visited
02502                                 if (vertexToScan == nVertices)
02503                                 {
02504                                         inTheStack. free ();
02505                                         return nComponents;
02506                                 }
02507 
02508                                 // set the new vertex
02509                                 vertex = vertexToScan ++;
02510 
02511                                 // mark the current vertex as visited
02512                                 // and initialize its low link
02513                                 dfsIndex [vertex] = ++ discoveryTime;
02514                                 lowLink [vertex] = discoveryTime;
02515 
02516                                 // push this vertex on the stack
02517                                 s_nodes. push (vertex);
02518                                 inTheStack. set (vertex);
02519 
02520                                 // determine the edges to be visited
02521                                 edge = 0;
02522                                 maxedge = g. countEdges (vertex);
02523                         }
02524 
02525                         // otherwise trace back to the previous level
02526                         else
02527                         {
02528                                 // remember the current low link index
02529                                 int_t lowLink2 = lowLink [vertex];
02530 
02531                                 // restore the variables
02532                                 vertex = s_vertex. top ();
02533                                 s_vertex. pop ();
02534                                 edge = s_edge. top ();
02535                                 s_edge. pop ();
02536                                 maxedge = s_maxedge. top ();
02537                                 s_maxedge. pop ();
02538 
02539                                 // update the current low link index
02540                                 if (lowLink [vertex] > lowLink2)
02541                                         lowLink [vertex] = lowLink2;
02542                         }
02543                 }
02544 
02545                 // analyse the next edge coming out from the current vertex
02546                 else
02547                 {
02548                         // determine the next vertex
02549                         int_t next = g. getEdge (vertex, edge ++);
02550 
02551                         // go to a deeper recursion level if unvisited
02552                         if (dfsIndex [next] == 0)
02553                         {
02554                                 // store the variables at the stacks
02555                                 s_vertex. push (vertex);
02556                                 s_edge. push (edge);
02557                                 s_maxedge. push (maxedge);
02558 
02559                                 // set the new vertex
02560                                 vertex = next;
02561                         
02562                                 // mark the new vertex as visited
02563                                 dfsIndex [vertex] = ++ discoveryTime;
02564                                 lowLink [vertex] = discoveryTime;
02565                         
02566                                 // push this vertex on the stack
02567                                 s_nodes. push (vertex);
02568                                 inTheStack. set (vertex);
02569 
02570                                 // determine the edges to be visited
02571                                 edge = 0;
02572                                 maxedge = g. countEdges (vertex);
02573                         }
02574 
02575                         // update the low link index if the vertex has been
02576                         // visited and is currently in the stack of nodes
02577                         else if (inTheStack. test (next))
02578                         {
02579                                 if (lowLink [vertex] > dfsIndex [next])
02580                                         lowLink [vertex] = dfsIndex [next];
02581                         }
02582                 }
02583         }
02584 
02585         // finalize and return the number of strongly connected components
02586         inTheStack. free ();
02587         return nComponents;
02588 } /* SCC_Tarjan */
02589 
02590 // --------------------------------------------------
02591 
02592 template <class wType>
02593 wType diGraph<wType>::minMeanCycleWeight (diGraph<wType> *transposed) const
02594 {
02595         // find the strongly connected components of the graph
02596         multitable<int_t> compVertices, compEnds;
02597         bool copyweights = !!transposed;
02598         int_t countSCC = SCC (*this, compVertices, compEnds,
02599                 static_cast<diGraph<wType> *> (0), transposed, copyweights);
02600         if (countSCC <= 0)
02601                 return 0;
02602 
02603         // compute the maximum size of each strongly connected component
02604         int_t maxCompSize = compEnds [0];
02605         for (int_t comp = 1; comp < countSCC; ++ comp)
02606         {
02607                 int_t compSize = compEnds [comp] - compEnds [comp - 1];
02608                 if (maxCompSize < compSize)
02609                         maxCompSize = compSize;
02610         }
02611 
02612         // allocate arrays for weights and bit fields
02613         wType **F = new wType * [maxCompSize + 1];
02614         BitField *finite = new BitField [maxCompSize + 1];
02615         for (int_t i = 0; i <= maxCompSize; ++ i)
02616         {
02617                 F [i] = new wType [maxCompSize];
02618                 finite [i]. allocate (maxCompSize);
02619         }
02620 
02621         // compute the numbers of vertices in each component
02622         int_t *numbers = new int_t [nVertices];
02623         int_t *components = new int_t [nVertices];
02624         for (int_t i = 0; i < nVertices; ++ i)
02625                 components [i] = -1;
02626         int_t offset = 0;
02627         for (int_t comp = 0; comp < countSCC; ++ comp)
02628         {
02629                 int_t maxOffset = compEnds [comp];
02630                 int_t pos = offset;
02631                 for (; pos < maxOffset; ++ pos)
02632                 {
02633                         numbers [compVertices [pos]] = pos - offset;
02634                         components [compVertices [pos]] = comp;
02635                 }
02636                 offset = pos;
02637         }
02638 
02639         // compute the minimum mean cycle weight for each component
02640         wType minWeight (0);
02641         for (int_t comp = 0; comp < countSCC; ++ comp)
02642         {
02643                 int_t offset = comp ? compEnds [comp - 1] :
02644                         static_cast<int_t> (0);
02645                 int_t compSize = compEnds [comp] - offset;
02646                 for (int_t i = 0; i <= compSize; ++ i)
02647                         finite [i]. clearall (compSize);
02648                 F [0] [0] = 0;
02649                 finite [0]. set (0);
02650                 // compute path progressions of given length
02651                 for (int_t len = 1; len <= compSize; ++ len)
02652                 {
02653                         // process source vertices
02654                         for (int_t i = 0; i < compSize; ++ i)
02655                         {
02656                                 if (!finite [len - 1]. test (i))
02657                                         continue;
02658                                 wType prevWeight = F [len - 1] [i];
02659                                 int_t source = compVertices [offset + i];
02660 
02661                                 // process target vertices (and edges)
02662                                 int_t edgeOffset = source ?
02663                                         edgeEnds [source - 1] :
02664                                         static_cast<int_t> (0);
02665                                 int_t edgeEnd = edgeEnds [source];
02666                                 for (; edgeOffset < edgeEnd; ++ edgeOffset)
02667                                 {
02668                                         int_t target = edges [edgeOffset];
02669                                         if (components [target] != comp)
02670                                                 continue;
02671                                         int_t j = numbers [target];
02672                                         wType newWeight = prevWeight +
02673                                                 weights [edgeOffset];
02674                                         if (!finite [len]. test (j))
02675                                         {
02676                                                 finite [len]. set (j);
02677                                                 F [len] [j] = newWeight;
02678                                         }
02679                                         else if (newWeight < F [len] [j])
02680                                                 F [len] [j] = newWeight;
02681                                 }
02682                         }
02683                 }
02684 
02685                 // compute the minimum mean cycle weight for this component
02686                 wType minCompWeight (0);
02687                 bool firstMin = true;
02688                 for (int_t vert = 0; vert < compSize; ++ vert)
02689                 {
02690                         if (!finite [compSize]. test (vert))
02691                                 continue;
02692                         bool firstAverage = true;
02693                         wType maxAverage = 0;
02694                         for (int_t first = 0; first < compSize; ++ first)
02695                         {
02696                                 if (!finite [first]. test (vert))
02697                                         continue;
02698                                 wType average = (F [compSize] [vert] -
02699                                         F [first] [vert]) /
02700                                         (compSize - first);
02701                                 if (firstAverage)
02702                                 {
02703                                         maxAverage = average;
02704                                         firstAverage = false;
02705                                 }
02706                                 else if (maxAverage < average)
02707                                         maxAverage = average;
02708                         }
02709                         if (firstMin || (maxAverage < minCompWeight))
02710                         {
02711                                 if (firstAverage)
02712                                         throw "Bug in Karp's algorithm";
02713                                 minCompWeight = maxAverage;
02714                                 firstMin = false;
02715                         }
02716                 }
02717 
02718                 // make a correction to the total minimum if necessary
02719                 if (!comp || (minCompWeight < minWeight))
02720                         minWeight = minCompWeight;
02721         }
02722 
02723         // release allocated memory
02724         delete [] components;
02725         delete [] numbers;
02726         for (int_t i = 0; i < maxCompSize; ++ i)
02727         {
02728                 finite [i]. free ();
02729                 delete [] (F [i]);
02730         }
02731         delete [] finite;
02732         delete [] F;
02733 
02734         // return the computed minimum
02735         return minWeight;
02736 } /* diGraph::minMeanCycleWeight */
02737 
02738 template <class wType>
02739 template <class roundType>
02740 wType diGraph<wType>::minMeanCycleWeight (const roundType &rounding,
02741         diGraph<wType> *transposed) const
02742 {
02743         // find the strongly connected components of the graph
02744         multitable<int_t> compVertices, compEnds;
02745         bool copyweights = !!transposed;
02746         int_t countSCC = SCC (*this, compVertices, compEnds,
02747                 static_cast<diGraph<wType> *> (0), transposed, copyweights);
02748         if (countSCC <= 0)
02749                 return 0;
02750 
02751         // compute the maximum size of each strongly connected component
02752         int_t maxCompSize = compEnds [0];
02753         for (int_t comp = 1; comp < countSCC; ++ comp)
02754         {
02755                 int_t compSize = compEnds [comp] - compEnds [comp - 1];
02756                 if (maxCompSize < compSize)
02757                         maxCompSize = compSize;
02758         }
02759 
02760         // allocate arrays for weights and bit fields
02761         wType **FUpper = new wType * [maxCompSize + 1];
02762         wType **FLower = new wType * [maxCompSize + 1];
02763         BitField *finite = new BitField [maxCompSize + 1];
02764         for (int_t i = 0; i <= maxCompSize; ++ i)
02765         {
02766                 FUpper [i] = new wType [maxCompSize];
02767                 FLower [i] = new wType [maxCompSize];
02768                 finite [i]. allocate (maxCompSize);
02769         }
02770 
02771         // compute the numbers of vertices in each component
02772         int_t *numbers = new int_t [nVertices];
02773         int_t *components = new int_t [nVertices];
02774         for (int_t i = 0; i < nVertices; ++ i)
02775                 components [i] = -1;
02776         int_t offset = 0;
02777         for (int_t comp = 0; comp < countSCC; ++ comp)
02778         {
02779                 int_t maxOffset = compEnds [comp];
02780                 int_t pos = offset;
02781                 for (; pos < maxOffset; ++ pos)
02782                 {
02783                         numbers [compVertices [pos]] = pos - offset;
02784                         components [compVertices [pos]] = comp;
02785                 }
02786                 offset = pos;
02787         }
02788 
02789         // compute the minimum mean cycle weight for each component
02790         wType minWeight (0);
02791         for (int_t comp = 0; comp < countSCC; ++ comp)
02792         {
02793                 int_t offset = comp ? compEnds [comp - 1] :
02794                         static_cast<int_t> (0);
02795                 int_t compSize = compEnds [comp] - offset;
02796                 for (int_t i = 0; i <= compSize; ++ i)
02797                         finite [i]. clearall (compSize);
02798                 FUpper [0] [0] = 0;
02799                 FLower [0] [0] = 0;
02800                 finite [0]. set (0);
02801                 // compute path progressions of given length
02802                 for (int_t len = 1; len <= compSize; ++ len)
02803                 {
02804                         // process source vertices
02805                         for (int_t i = 0; i < compSize; ++ i)
02806                         {
02807                                 if (!finite [len - 1]. test (i))
02808                                         continue;
02809                                 wType prevUpper = FUpper [len - 1] [i];
02810                                 wType prevLower = FLower [len - 1] [i];
02811                                 int_t source = compVertices [offset + i];
02812 
02813                                 // process target vertices (and edges)
02814                                 int_t edgeOffset = source ?
02815                                         edgeEnds [source - 1] :
02816                                         static_cast<int_t> (0);
02817                                 int_t edgeEnd = edgeEnds [source];
02818                                 for (; edgeOffset < edgeEnd; ++ edgeOffset)
02819                                 {
02820                                         int_t target = edges [edgeOffset];
02821                                         if (components [target] != comp)
02822                                                 continue;
02823                                         int_t j = numbers [target];
02824                                         wType newUpper = rounding. add_up
02825                                                 (prevUpper,
02826                                                 weights [edgeOffset]);
02827                                         wType newLower = rounding. add_down
02828                                                 (prevLower,
02829                                                 weights [edgeOffset]);
02830                                         if (!finite [len]. test (j))
02831                                         {
02832                                                 finite [len]. set (j);
02833                                                 FUpper [len] [j] = newUpper;
02834                                                 FLower [len] [j] = newLower;
02835                                         }
02836                                         else
02837                                         {
02838                                                 wType &curUpper =
02839                                                         FUpper [len] [j];
02840                                                 if (newUpper < curUpper)
02841                                                         curUpper = newUpper;
02842                                                 wType &curLower =
02843                                                         FLower [len] [j];
02844                                                 if (newLower < curLower)
02845                                                         curLower = newLower;
02846                                         }
02847                                 }
02848                         }
02849                 }
02850 
02851                 // compute the minimum mean cycle weight for this component
02852                 wType minCompWeight (0);
02853                 bool firstMin = true;
02854                 for (int_t vert = 0; vert < compSize; ++ vert)
02855                 {
02856                         if (!finite [compSize]. test (vert))
02857                                 continue;
02858                         bool firstAverage = true;
02859                         wType maxAverage = 0;
02860                         for (int_t first = 0; first < compSize; ++ first)
02861                         {
02862                                 if (!finite [first]. test (vert))
02863                                         continue;
02864                                 const wType diff = rounding. sub_down
02865                                         (FLower [compSize] [vert],
02866                                         FUpper [first] [vert]);
02867                                 wType average = rounding. div_down
02868                                         (diff, compSize - first);
02869                                 if (firstAverage)
02870                                 {
02871                                         maxAverage = average;
02872                                         firstAverage = false;
02873                                 }
02874                                 else if (maxAverage < average)
02875                                         maxAverage = average;
02876                         }
02877                         if (firstMin || (maxAverage < minCompWeight))
02878                         {
02879                                 if (firstAverage)
02880                                         throw "Bug in Karp's algorithm";
02881                                 minCompWeight = maxAverage;
02882                                 firstMin = false;
02883                         }
02884                 }
02885 
02886                 // make a correction to the total minimum if necessary
02887                 if (!comp || (minCompWeight < minWeight))
02888                         minWeight = minCompWeight;
02889         }
02890 
02891         // release allocated memory
02892         delete [] components;
02893         delete [] numbers;
02894         for (int_t i = 0; i < maxCompSize; ++ i)
02895         {
02896                 finite [i]. free ();
02897                 delete [] (FUpper [i]);
02898                 delete [] (FLower [i]);
02899         }
02900         delete [] finite;
02901         delete [] FUpper;
02902         delete [] FLower;
02903 
02904         // return the computed minimum
02905         return minWeight;
02906 } /* diGraph::minMeanCycleWeight_intv */
02907 
02908 template <class wType>
02909 template <class arrayType, class roundType>
02910 wType diGraph<wType>::minMeanPathWeight (const roundType &rounding,
02911         const arrayType &starting, int_t n) const
02912 {
02913         // allocate arrays for weights and bit fields
02914         const int nIndices = 2;
02915         wType **F = new wType * [nIndices];
02916         BitField *finite = new BitField [nIndices];
02917         for (int i = 0; i < nIndices; ++ i)
02918         {
02919                 F [i] = new wType [nVertices];
02920                 finite [i]. allocate (nVertices);
02921         }
02922 
02923         // set the zero path lengths from the initial vertices
02924         for (int_t i = 0; i < n; ++ i)
02925         {
02926                 int_t v = starting [i];
02927                 if ((v < 0) || (v >= nVertices))
02928                         throw "Starting vertex out of range.";
02929                 F [0] [v] = 0;
02930                 finite [0]. set (v);
02931         }
02932 
02933         // compute path progressions of given length and average weights
02934         wType minWeight (0);
02935         bool firstWeight = true;
02936         for (int_t len = 1; len <= nVertices; ++ len)
02937         {
02938                 // determine the indices for previous and current paths
02939                 int_t prevIndex = (len - 1) & 1;
02940                 int_t curIndex = len & 1;
02941                 finite [curIndex]. clearall (nVertices);
02942 
02943                 // process source vertices
02944                 for (int_t source = 0; source < nVertices; ++ source)
02945                 {
02946                         if (!finite [prevIndex]. test (source))
02947                                 continue;
02948                         wType prevWeight = F [prevIndex] [source];
02949 
02950                         // process target vertices (and edges)
02951                         int_t edgeOffset = source ?
02952                                 edgeEnds [source - 1] :
02953                                 static_cast<int_t> (0);
02954                         int_t edgeEnd = edgeEnds [source];
02955                         for (; edgeOffset < edgeEnd; ++ edgeOffset)
02956                         {
02957                                 int_t target = edges [edgeOffset];
02958                                 wType newWeight = rounding. add_down
02959                                         (prevWeight, weights [edgeOffset]);
02960                                 if (!finite [curIndex]. test (target))
02961                                 {
02962                                         finite [curIndex]. set (target);
02963                                         F [curIndex] [target] = newWeight;
02964                                 }
02965                                 else if (newWeight < F [curIndex] [target])
02966                                         F [curIndex] [target] = newWeight;
02967                         }
02968                 }
02969 
02970                 // update the minimum mean path weight
02971                 for (int_t vert = 0; vert < nVertices; ++ vert)
02972                 {
02973                         if (!finite [curIndex]. test (vert))
02974                                 continue;
02975                         wType average = rounding. div_down
02976                                 (F [curIndex] [vert], len);
02977                         if (firstWeight)
02978                         {
02979                                 minWeight = average;
02980                                 firstWeight = false;
02981                         }
02982                         else if (average < minWeight)
02983                                 minWeight = average;
02984                 }
02985         }
02986 
02987         // release allocated memory
02988         for (int i = 0; i < nIndices; ++ i)
02989         {
02990                 finite [i]. free ();
02991                 delete [] (F [i]);
02992         }
02993         delete [] finite;
02994         delete [] F;
02995 
02996         // return the computed minimum
02997         return minWeight;
02998 } /* diGraph::minMeanPathWeight */
02999 
03000 template <class wType>
03001 template <class arrayType>
03002 wType diGraph<wType>::minMeanPathWeight (const arrayType &starting, int_t n)
03003         const
03004 {
03005         const dummyRounding<wType> rounding = dummyRounding<wType> ();
03006         return minMeanPathWeight (rounding, starting, n);
03007 } /* diGraph::minMeanPathWeight */
03008 
03009 
03010 // --------------------------------------------------
03011 // ------------- OTHER GRAPH ALGORITHMS -------------
03012 // --------------------------------------------------
03013 
03014 /// A helper function for "collapseVertices".
03015 template <class wType>
03016 inline int_t addRenumEdges (const diGraph<wType> &g, int_t vertex,
03017         const int_t *newNum, int_t cur, int_t *srcVert,
03018         diGraph<wType> &result)
03019 {
03020         int_t nEdges = g. countEdges (vertex);
03021         // add all the edges that start at the component
03022         for (int_t edge = 0; edge < nEdges; ++ edge)
03023         {
03024                 // determine the dest. vertex of the edge
03025                 int_t dest = newNum [g. getEdge (vertex, edge)];
03026                 // if this is an edge to self, then ignore it
03027                 if (dest == cur)
03028                         continue;
03029                 // if the edge has already been added,
03030                 // then skip it
03031                 if (srcVert [dest] == cur)
03032                         continue;
03033                 // remember that the dest. vertex has an edge
03034                 // pointing out from the current vertex
03035                 srcVert [dest] = cur;
03036                 // add the edge to the result graph
03037                 result. addEdge (dest);
03038         }
03039         return 0;
03040 } /* addRenumEdges */
03041 
03042 /// Collapses disjoint subsets of vertices to single vertices and creates
03043 /// a corresponding graph in which the first vertices come from
03044 /// the collapsed ones. The result graph must be initially empty.
03045 /// Saves translation table from old to new vertex numbers.
03046 /// The table must have sufficient length.
03047 template <class wType, class Table1, class Table2>
03048 inline int_t collapseVertices (const diGraph<wType> &g,
03049         int_t compNum, Table1 &compVertices, Table2 &compEnds,
03050         diGraph<wType> &result, int_t *newNumbers = 0)
03051 {
03052         // do nothing if the input graph is empty
03053         int_t nVert = g. countVertices ();
03054         if (!nVert)
03055                 return 0;
03056 
03057         // compute the new numbers of vertices: newNum [i] is the
03058         // number of vertex 'i' from 'g' in the resulting graph
03059         int_t *newNum = newNumbers ? newNumbers : new int_t [nVert];
03060         for (int_t i = 0; i < nVert; ++ i)
03061                 newNum [i] = -1;
03062         int_t cur = 0; // current vertex number in the new graph
03063         int_t pos = 0; // position in compVertices
03064         while (cur < compNum)
03065         {
03066                 int_t posEnd = compEnds [cur];
03067                 while (pos < posEnd)
03068                         newNum [compVertices [pos ++]] = cur;
03069                 ++ cur;
03070         }
03071         for (int_t i = 0; i < nVert; ++ i)
03072         {
03073                 if (newNum [i] < 0)
03074                         newNum [i] = cur ++;
03075         }
03076 
03077         // prepare a table to mark the most recent source vertex for edges:
03078         // srcVert [j] contains i if the edge i -> j has just been added
03079         int_t newVert = nVert - (compNum ? compEnds [compNum - 1] :
03080                 static_cast<int_t> (0)) + compNum;
03081         // debug:
03082         if (cur != newVert)
03083                 throw "DEBUG: Wrong numbers.";
03084         int_t *srcVert = new int_t [newVert];
03085         for (int_t i = 0; i < newVert; ++ i)
03086                 srcVert [i] = -1;
03087 
03088         // scan the input graph and create the resulting graph: begin with
03089         // the vertices in the collapsed groups
03090         cur = 0, pos = 0;
03091         while (cur < compNum)
03092         {
03093                 // add a new vertex to the result graph
03094                 result. addVertex ();
03095                 // for all the vertices in this component...
03096                 int_t posEnd = compEnds [cur];
03097                 while (pos < posEnd)
03098                 {
03099                         // take the next vertex from the component
03100                         int_t vertex = compVertices [pos ++];
03101                         // add the appropriate edges to the result graph
03102                         addRenumEdges (g, vertex, newNum, cur,
03103                                 srcVert, result);
03104                 }
03105                 ++ cur;
03106         }
03107         // process the remaining vertices of the graph
03108         for (int_t vertex = 0; vertex < nVert; ++ vertex)
03109         {
03110                 // debug
03111                 if (newNum [vertex] > cur)
03112                         throw "DEBUG: Wrong order.";
03113                 // if the vertex has already been processed, skip it
03114                 if (newNum [vertex] != cur)
03115                         continue;
03116                 // add the appropriate edges to the result graph
03117                 result. addVertex ();
03118                 addRenumEdges (g, vertex, newNum, cur, srcVert, result);
03119                 ++ cur;
03120         }
03121 
03122         // free memory and exit
03123         delete [] srcVert;
03124         if (!newNumbers)
03125                 delete [] newNum;
03126         return 0;
03127 } /* collapseVertices */
03128 
03129 /// A helper function for "collapseVertices2".
03130 template <class wType, class TabSets>
03131 inline int_t addRenumEdges2 (const diGraph<wType> &g, int_t vertex,
03132         const int_t *newNum, TabSets &numSets,
03133         int_t cur, int_t *srcVert, diGraph<wType> &result)
03134 {
03135         int_t nEdges = g. countEdges (vertex);
03136 
03137         // add all the edges that start at this vertex
03138         for (int_t edge = 0; edge < nEdges; ++ edge)
03139         {
03140                 // determine the dest. vertex of the edge
03141                 int_t destNumber = newNum [g. getEdge (vertex, edge)];
03142 
03143                 // consider all the destination vertices
03144                 int_t destSize = (destNumber < 0) ?
03145                         numSets [-destNumber]. size () :
03146                         static_cast<int_t> (1);
03147                 for (int_t i = 0; i < destSize; ++ i)
03148                 {
03149                         // determine the consecutive destination vertex
03150                         int_t dest = (destNumber < 0) ?
03151                                 numSets [-destNumber] [i] : destNumber;
03152 
03153                         // if this is an edge to self, then ignore it
03154                         if (dest == cur)
03155                                 continue;
03156 
03157                         // if the edge has already been added,
03158                         // then skip it
03159                         if (srcVert [dest] == cur)
03160                                 continue;
03161 
03162                         // add the edge to the result graph
03163                         result. addEdge (dest);
03164                 }
03165         }
03166         return 0;
03167 } /* addRenumEdges2 */
03168 
03169 /// Collapses (possibly non-disjoint) subsets of vertices
03170 /// to single vertices and creates a corresponding graph
03171 /// in which the first vertices come from the collapsed ones.
03172 /// The result graph must be initially empty.
03173 template <class wType, class Table1, class Table2>
03174 inline int_t collapseVertices2 (const diGraph<wType> &g,
03175         int_t compNum, Table1 &compVertices, Table2 &compEnds,
03176         diGraph<wType> &result)
03177 {
03178         // do nothing if the input graph is empty
03179         int_t nVert = g. countVertices ();
03180         if (!nVert)
03181                 return 0;
03182 
03183         // compute the new numbers of vertices: newNum [i] is the
03184         // number of vertex 'i' from 'g' in the resulting graph,
03185         // unless negative; then it points to a set of numbers,
03186         // with -1 meaning "not defined, yet"
03187         int_t *newNum = new int_t [nVert];
03188         for (int_t i = 0; i < nVert; ++ i)
03189                 newNum [i] = -1;
03190 
03191         // sets of numbers of vertices; the numbers refering to the sets
03192         // begin with 2, thus the first two sets are skipped
03193         multitable<hashedset<int_t> > numSets;
03194         // the number of the set waiting to be used next
03195         int_t numSetCur = 2;
03196 
03197         int_t cur = 0; // current vertex number in the new graph
03198         int_t pos = 0; // position in compVertices
03199         while (cur < compNum)
03200         {
03201                 int_t posEnd = compEnds [cur];
03202                 while (pos < posEnd)
03203                 {
03204                         int_t number = compVertices [pos ++];
03205                         if (newNum [number] == -1)
03206                                 newNum [number] = cur;
03207                         else if (newNum [number] < 0)
03208                                 numSets [-newNum [number]]. add (cur);
03209                         else
03210                         {
03211                                 numSets [numSetCur]. add (newNum [number]);
03212                                 numSets [numSetCur]. add (cur);
03213                                 newNum [number] = -(numSetCur ++);
03214                         }
03215                 }
03216                 ++ cur;
03217         }
03218         for (int_t i = 0; i < nVert; ++ i)
03219         {
03220                 if (newNum [i] == -1)
03221                         newNum [i] = cur ++;
03222         }
03223 
03224         // prepare a table to mark the most recent source vertex for edges:
03225         // srcVert [j] contains i if the edge i -> j has just been added
03226         int_t newVert = cur;
03227         int_t *srcVert = new int_t [newVert];
03228         for (int_t i = 0; i < newVert; ++ i)
03229                 srcVert [i] = -1;
03230 
03231         // scan the input graph and create the resulting graph: begin with
03232         // the vertices in the collapsed groups
03233         cur = 0;
03234         pos = 0;
03235         while (cur < compNum)
03236         {
03237                 // add a new vertex to the result graph
03238                 result. addVertex ();
03239 
03240                 // for all the vertices in this component...
03241                 int_t posEnd = compEnds [cur];
03242                 while (pos < posEnd)
03243                 {
03244                         // take the next vertex from the component
03245                         int_t vertex = compVertices [pos ++];
03246 
03247                         // add the appropriate edges to the result graph
03248                         addRenumEdges2 (g, vertex, newNum, numSets, cur,
03249                                 srcVert, result);
03250                 }
03251                 ++ cur;
03252         }
03253 
03254         // process the remaining vertices of the graph
03255         for (int_t vertex = 0; vertex < nVert; ++ vertex)
03256         {
03257                 // debug
03258                 if (newNum [vertex] > cur)
03259                         throw "DEBUG: Wrong order 2.";
03260 
03261                 // if the vertex has already been processed, skip it
03262                 if (newNum [vertex] != cur)
03263                         continue;
03264 
03265                 // add the appropriate edges to the result graph
03266                 result. addVertex ();
03267                 addRenumEdges2 (g, vertex, newNum, numSets, cur,
03268                         srcVert, result);
03269                 ++ cur;
03270         }
03271 
03272         // free memory and exit
03273         delete [] srcVert;
03274         delete [] newNum;
03275         return 0;
03276 } /* collapseVertices2 */
03277 
03278 /// Computes the graph that represents connections between a number
03279 /// of the first vertices in the given graph.
03280 /// The vertices of the result graph are the first "nVert" vertices
03281 /// from the source graph. An edge is added to the new graph
03282 /// iff there exists a path from one vertex to another,
03283 /// without going through any other vertex in that set.
03284 /// Runs DFS starting from each of the first "nVert" vertices,
03285 /// and thus may be a little inefficient in some cases.
03286 template <class wType>
03287 inline int_t connectionGraph (const diGraph<wType> &g, int_t nVert,
03288         diGraph<wType> &result)
03289 {
03290         // remember the number of vertices in the input graph
03291         int_t nVertG = g. countVertices ();
03292         if (!nVertG)
03293                 return 0;
03294 
03295         // prepare a bitfield in which visited vertices will be marked
03296         BitField visited;
03297         visited. allocate (nVertG);
03298         visited. clearall (nVertG);
03299 
03300         // run DFS for each of the starting vertices
03301         for (int_t startVertex = 0; startVertex < nVert; ++ startVertex)
03302         {
03303                 // add this vertex to the resulting graph
03304                 result. addVertex ();
03305 
03306                 // prepare a counter and a stack of visited vertices
03307                 int_t nVisited = 0;
03308                 std::stack<int_t> s_visited;
03309 
03310                 // prepare stacks for the DFS algorithm
03311                 std::stack<int_t> s_vertex;
03312                 std::stack<int_t> s_edge;
03313 
03314                 // mark the starting vertex as visited
03315                 visited. set (startVertex);
03316                 s_visited. push (startVertex);
03317                 ++ nVisited;
03318 
03319                 // use DFS to visit vertices reachable from that vertex
03320                 int_t vertex = startVertex;
03321                 int_t edge = 0;
03322                 while (1)
03323                 {
03324                         // go back with the recursion
03325                         // if all the edges have been processed
03326                         if (edge >= g. countEdges (vertex))
03327                         {
03328                                 // if this was the top recursion level
03329                                 // then quit
03330                                 if (s_vertex. empty ())
03331                                         break;
03332 
03333                                 // return from the recursion
03334                                 vertex = s_vertex. top ();
03335                                 s_vertex. pop ();
03336                                 edge = s_edge. top ();
03337                                 s_edge. pop ();
03338                                 continue;
03339                         }
03340 
03341                         // take the next vertex
03342                         int_t next = g. getEdge (vertex, edge ++);
03343 
03344                         // go to the deeper recursion level if not visited
03345                         if (!visited. test (next))
03346                         {
03347                                 // add an edge to the result if necessary
03348                                 if (next < nVert)
03349                                         result. addEdge (next);
03350 
03351                                 // store the previous variables at the stacks
03352                                 s_vertex. push (vertex);
03353                                 s_edge. push (edge);
03354 
03355                                 // take the new vertex
03356                                 vertex = next;
03357                                 edge = 0;
03358 
03359                                 // mark the new vertex as visited
03360                                 visited. set (vertex);
03361                                 s_visited. push (vertex);
03362                                 ++ nVisited;
03363                         }
03364                 }
03365 
03366                 // if this was the last vertex then skip the rest
03367                 if (startVertex == nVert - 1)
03368                         break;
03369 
03370                 // mark each vertex as non-visited
03371                 if (nVisited > (nVertG >> 6))
03372                 {
03373                         visited. clearall (nVertG);
03374                 }
03375                 else
03376                 {
03377                         while (!s_visited. empty ())
03378                         {
03379                                 int_t vertex = s_visited. top ();
03380                                 s_visited. pop ();
03381                                 visited. clear (vertex);
03382                         }
03383                 }
03384         }
03385 
03386         // free memory and exit
03387         visited. free ();
03388         return 0;
03389 } /* connectionGraph */
03390 
03391 /// Computes the period of a strongly connected graph.
03392 /// The period of a graph is defined as the greatest common divisor
03393 /// of the lengths of all the cycles in the graph.
03394 /// Period 1 means that the graph is aperiodic.
03395 /// The complexity of this operation is linear (one DFS).
03396 template <class GraphType>
03397 inline int_t computePeriod (const GraphType &g)
03398 {
03399         // remember the number of vertices in the input graph
03400         int_t nVertG = g. countVertices ();
03401         if (!nVertG)
03402                 return 0;
03403 
03404         // prepare an array to record the depth of each visited vertex
03405         std::vector<int_t> visited (nVertG, 0);
03406 
03407         // run DFS starting at the first vertex
03408         // prepare stacks for the DFS algorithm
03409         std::stack<int_t> s_vertex;
03410         std::stack<int_t> s_edge;
03411 
03412         // mark the starting vertex as visited
03413         visited [0] = 1;
03414 
03415         // use DFS to visit vertices reachable from that vertex
03416         int_t vertex = 0;
03417         int_t edge = 0;
03418         int_t level = 1;
03419         int_t period = 0;
03420         while (1)
03421         {
03422                 // go back with the recursion
03423                 // if all the edges have been processed
03424                 if (edge >= g. countEdges (vertex))
03425                 {
03426                         // if this was the top recursion level then quit
03427                         if (s_vertex. empty ())
03428                                 break;
03429 
03430                         // return from the recursion
03431                         vertex = s_vertex. top ();
03432                         s_vertex. pop ();
03433                         edge = s_edge. top ();
03434                         s_edge. pop ();
03435                         -- level;
03436                         continue;
03437                 }
03438 
03439                 // take the next vertex
03440                 int_t next = g. getEdge (vertex, edge ++);
03441 
03442                 // update the GCD of the cycle lengths
03443                 if (visited [next])
03444                 {
03445                         // determine the period provided by the cycle
03446                         int_t newPeriod = visited [next] - level - 1;
03447                         if (newPeriod < 0)
03448                                 newPeriod = -newPeriod;
03449 
03450                         // compute the GCD of the old and new periods
03451                         int_t a = newPeriod;
03452                         int_t b = period;
03453                         while (b)
03454                         {
03455                                 int_t t = b;
03456                                 b = a % b;
03457                                 a = t;
03458                         }
03459                         period = a;
03460 
03461                         // if the graph turns out to be aperiodic
03462                         // then immediately return the result
03463                         if (period == 1)
03464                                 return period;
03465                 }
03466 
03467                 // go to the deeper recursion level if not visited
03468                 else
03469                 {
03470                         // store the previous variables at the stacks
03471                         s_vertex. push (vertex);
03472                         s_edge. push (edge);
03473 
03474                         // take the new vertex
03475                         vertex = next;
03476                         edge = 0;
03477                         ++ level;
03478 
03479                         // mark the new vertex as visited
03480                         visited [vertex] = level;
03481                 }
03482         }
03483 
03484         // return the computed peirod and exit
03485         return period;
03486 } /* computePeriod */
03487 
03488 /// Computes the graph that represents flow-induced relations on Morse sets.
03489 /// The vertices of the result graph are the first "nVert" vertices
03490 /// from the source graph. An edge is added to the new graph
03491 /// iff there exists a path from one vertex to another.
03492 /// Edges that come from the transitive closure are not added.
03493 template <class wType>
03494 inline int_t inclusionGraph (const diGraph<wType> &g, int_t nVert,
03495         diGraph<wType> &result)
03496 {
03497         // remember the number of vertices in the input graph
03498         int_t nVertG = g. countVertices ();
03499         if (!nVertG)
03500                 return 0;
03501 
03502         // mark each vertex as non-visited
03503         BitField visited;
03504         visited. allocate (nVertG);
03505         visited. clearall (nVertG);
03506 
03507         // create the sets of vertices reachable from each vertex
03508         BitSets lists (nVertG, nVert);
03509 
03510         // prepare stacks for the DFS algorithm
03511         std::stack<int_t> s_vertex;
03512         std::stack<int_t> s_edge;
03513 
03514         // use DFS to propagate the reachable sets
03515         int_t startVertex = 0;
03516         int_t vertex = 0;
03517         int_t edge = 0;
03518         visited. set (vertex);
03519         while (1)
03520         {
03521                 // go back with the recursion
03522                 // if all the edges have been processed
03523                 if (edge >= g. countEdges (vertex))
03524                 {
03525                         // if this was the top recursion level,
03526                         // then find another starting point or quit
03527                         if (s_vertex. empty ())
03528                         {
03529                                 while ((startVertex < nVertG) &&
03530                                         (visited. test (startVertex)))
03531                                         ++ startVertex;
03532                                 if (startVertex >= nVertG)
03533                                         break;
03534                                 vertex = startVertex;
03535                                 edge = 0;
03536                                 visited. set (vertex);
03537                                 continue;
03538                         }
03539 
03540                         // determine the previous vertex (to trace back to)
03541                         int_t prev = s_vertex. top ();
03542 
03543                         // add the list of the current vertex
03544                         // to the list of the previous vertex
03545                         // unless that was one of the first 'nVert' vertices
03546                         if (vertex >= nVert)
03547                                 lists. addset (prev, vertex);
03548 
03549                         // return from the recursion
03550                         vertex = prev;
03551                         s_vertex. pop ();
03552                         edge = s_edge. top ();
03553                         s_edge. pop ();
03554                         continue;
03555                 }
03556 
03557                 // take the next vertex
03558                 int_t next = g. getEdge (vertex, edge ++);
03559 
03560                 // add the number to the list if appropriate
03561                 if (next < nVert)
03562                         lists. add (vertex, next);
03563 
03564                 // go to the deeper recursion level if vertex not visited
03565                 if (!visited. test (next))
03566                 {
03567                         // store the previous variables at the stacks
03568                         s_vertex. push (vertex);
03569                         s_edge. push (edge);
03570 
03571                         // take the new vertex and mark as visited
03572                         vertex = next;
03573                         edge = 0;
03574                         visited. set (vertex);
03575                 }
03576 
03577                 // add the list of the next vertex to the current one
03578                 // if the vertex has already been visited before
03579                 // and the vertex is not one of the first 'nVert' ones
03580                 else if (next >= nVert)
03581                 {
03582                         lists. addset (vertex, next);
03583                 }
03584         }
03585 
03586         // create the result graph based on the lists of edges
03587         for (int_t vertex = 0; vertex < nVert; ++ vertex)
03588         {
03589                 result. addVertex ();
03590                 int_t edge = lists. get (vertex);
03591                 while (edge >= 0)
03592                 {
03593                         result. addEdge (edge);
03594                         edge = lists. get (vertex, edge + 1);
03595                 }
03596         }
03597 
03598         // free memory and exit
03599         visited. free ();
03600         return 0;
03601 } /* inclusionGraph */
03602 
03603 // --------------------------------------------------
03604 
03605 /// A more complicated version of the 'inclusionGraph' function.
03606 /// Computes the graph that represents flow-induced relations on Morse sets.
03607 /// The vertices of the result graph are the first "nVert" vertices
03608 /// from the source graph. An edge is added to the new graph
03609 /// iff there exists a path from one vertex to another.
03610 /// Edges that come from the transitive closure are not added.
03611 /// Records vertices along connecting orbits using operator () with the
03612 /// following arguments: source, target, vertex.
03613 template<class wType, class TConn>
03614 inline int_t inclusionGraph (const diGraph<wType> &g, int_t nVert,
03615         diGraph<wType> &result, TConn &connections)
03616 {
03617         // remember the number of vertices in the input graph
03618         int_t nVertG = g. countVertices ();
03619         if (!nVertG)
03620                 return 0;
03621 
03622         // mark each vertex as non-visited
03623         BitField visited;
03624         visited. allocate (nVertG);
03625         visited. clearall (nVertG);
03626 
03627         // create the sets of vertices reachable from each vertex
03628         BitSets forwardlists (nVertG, nVert);
03629 
03630         // prepare stacks for the DFS algorithm
03631         std::stack<int_t> s_vertex;
03632         std::stack<int_t> s_edge;
03633 
03634         // use DFS to propagate the reachable sets
03635         int_t startVertex = 0;
03636         int_t vertex = 0;
03637         int_t edge = 0;
03638         visited. set (vertex);
03639         while (1)
03640         {
03641                 // go back with the recursion
03642                 // if all the edges have been processed
03643                 if (edge >= g. countEdges (vertex))
03644                 {
03645                         // if this was the top recursion level,
03646                         // then find another starting point or quit
03647                         if (s_vertex. empty ())
03648                         {
03649                                 while ((startVertex < nVertG) &&
03650                                         (visited. test (startVertex)))
03651                                         ++ startVertex;
03652                                 if (startVertex >= nVertG)
03653                                         break;
03654                                 vertex = startVertex;
03655                                 edge = 0;
03656                                 visited. set (vertex);
03657                                 continue;
03658                         }
03659 
03660                         // determine the previous vertex (to trace back to)
03661                         int_t prev = s_vertex. top ();
03662 
03663                         // add the list of the current vertex
03664                         // to the list of the previous vertex
03665                         // unless that was one of the first 'nVert' vertices
03666                         if (vertex >= nVert)
03667                                 forwardlists. addset (prev, vertex);
03668 
03669                         // return from the recursion
03670                         vertex = prev;
03671                         s_vertex. pop ();
03672                         edge = s_edge. top ();
03673                         s_edge. pop ();
03674                         continue;
03675                 }
03676 
03677                 // take the next vertex
03678                 int_t next = g. getEdge (vertex, edge ++);
03679 
03680                 // add the number to the list if appropriate
03681                 if (next < nVert)
03682                         forwardlists. add (vertex, next);
03683 
03684                 // go to the deeper recursion level if vertex not visited
03685                 if (!visited. test (next))
03686                 {
03687                         // store the previous variables at the stacks
03688                         s_vertex. push (vertex);
03689                         s_edge. push (edge);
03690 
03691                         // take the new vertex and mark as visited
03692                         vertex = next;
03693                         edge = 0;
03694                         visited. set (vertex);
03695                 }
03696 
03697                 // add the list of the next vertex to the current one
03698                 // if the vertex has already been visited before
03699                 // and the vertex is not one of the first 'nVert' ones
03700                 else if (next >= nVert)
03701                 {
03702                         forwardlists. addset (vertex, next);
03703                 }
03704         }
03705 
03706         // ------------------------------------------
03707 
03708         // create the sets of vertices that can reach each vertex
03709         BitSets backlists (nVertG, nVert);
03710 
03711         diGraph<wType> gT;
03712         g. transpose (gT);
03713 
03714         // do a simple test for debugging purposes
03715         if (!s_vertex. empty () || !s_edge. empty ())
03716                 throw "DEBUG: Nonempty stacks in 'inclusionGraph'.";
03717 
03718         // mark all the vertices as unvisited for the backward run
03719         visited. clearall (nVertG);
03720 
03721         // use DFS to propagate the reachable sets
03722         startVertex = 0;
03723         vertex = 0;
03724         edge = 0;
03725         visited. set (vertex);
03726         while (1)
03727         {
03728                 // go back with the recursion
03729                 // if all the edges have been processed
03730                 if (edge >= gT. countEdges (vertex))
03731                 {
03732                         // if this was the top recursion level,
03733                         // then find another starting point or quit
03734                         if (s_vertex. empty ())
03735                         {
03736                                 while ((startVertex < nVertG) &&
03737                                         (visited. test (startVertex)))
03738                                         ++ startVertex;
03739                                 if (startVertex >= nVertG)
03740                                         break;
03741                                 vertex = startVertex;
03742                                 edge = 0;
03743                                 visited. set (vertex);
03744                                 continue;
03745                         }
03746 
03747                         // determine the previous vertex (to trace back to)
03748                         int_t prev = s_vertex. top ();
03749 
03750                         // add the list of the current vertex
03751                         // to the list of the previous vertex
03752                         // unless that was one of the first 'nVert' vertices
03753                         if (vertex >= nVert)
03754                                 backlists. addset (prev, vertex);
03755 
03756                         // return from the recursion
03757                         vertex = prev;
03758                         s_vertex. pop ();
03759                         edge = s_edge. top ();
03760                         s_edge. pop ();
03761                         continue;
03762                 }
03763 
03764                 // take the next vertex
03765                 int_t next = gT. getEdge (vertex, edge ++);
03766 
03767                 // add the number to the list if appropriate
03768                 if (next < nVert)
03769                         backlists. add (vertex, next);
03770 
03771                 // go to the deeper recursion level if vertex not visited
03772                 if (!visited. test (next))
03773                 {
03774                         // store the previous variables at the stacks
03775                         s_vertex. push (vertex);
03776                         s_edge. push (edge);
03777 
03778                         // take the new vertex and mark as visited
03779                         vertex = next;
03780                         edge = 0;
03781                         visited. set (vertex);
03782                 }
03783 
03784                 // add the list of the next vertex to the current one
03785                 // if the vertex has already been visited before
03786                 // and the vertex is not one of the first 'nVert' ones
03787                 else if (next >= nVert)
03788                 {
03789                         backlists. addset (vertex, next);
03790                 }
03791         }
03792 
03793         // ------------------------------------------
03794 
03795         // record the connections to the obtained data structure
03796         for (int_t v = nVert; v < nVertG; ++ v)
03797         {
03798                 int_t bvertex = backlists. get (v);
03799                 while (bvertex >= 0)
03800                 {
03801                         int_t fvertex = forwardlists. get (v);
03802                         while (fvertex >= 0)
03803                         {
03804                                 connections (bvertex, fvertex, v);
03805                                 fvertex = forwardlists. get (v, fvertex + 1);
03806                         }
03807                         bvertex = backlists. get (v, bvertex + 1);
03808                 }
03809         }
03810 
03811         // ------------------------------------------
03812         
03813         // create the result graph based on the lists of edges
03814         for (int_t vertex = 0; vertex < nVert; ++ vertex)
03815         {
03816                 result. addVertex ();
03817                 int_t edge = forwardlists. get (vertex);
03818                 while (edge >= 0)
03819                 {
03820                         result. addEdge (edge);
03821                         edge = forwardlists. get (vertex, edge + 1);
03822                 }
03823         }
03824 
03825         // free memory and exit
03826         visited. free ();
03827         return 0;
03828 } /* inclusionGraph */
03829 
03830 // --------------------------------------------------
03831 
03832 /// Creates the adjacency matrix of the given graph.
03833 /// m [i] [j] is set to 1 if the graph g contains the edge i -> j,
03834 /// using the assignment operator.
03835 template <class wType, class matrix>
03836 inline void graph2matrix (const diGraph<wType> &g, matrix &m)
03837 {
03838         int_t nVert = g. countVertices ();
03839         for (int_t v = 0; v < nVert; ++ v)
03840         {
03841                 int_t nEdges = g. countEdges (v);
03842                 for (int_t e = 0; e < nEdges; ++ e)
03843                 {
03844                         int_t w = g. getEdge (v, e);
03845                         m [v] [w] = 1;
03846                 }
03847         }
03848         return;
03849 } /* graph2matrix */
03850 
03851 /// Creates a graph based on the given adjacency matrix.
03852 /// If m [i] [j] is nonzero, then the edge i -> j is added to the graph.
03853 /// It is assumed that the graph g is initially empty.
03854 /// The size of the matrix (the number of vertices), n, must be given.
03855 template <class wType, class matrix>
03856 inline void matrix2graph (const matrix &m, int_t n, diGraph<wType> &g)
03857 {
03858         for (int_t v = 0; v < n; ++ v)
03859         {
03860                 g. addVertex ();
03861                 for (int_t w = 0; w < n; ++ w)
03862                         if (m [v] [w])
03863                                 g. addEdge (w);
03864         }
03865         return;
03866 } /* matrix2graph */
03867 
03868 // --------------------------------------------------
03869 
03870 /// Computes the transitive closure of an acyclic graph
03871 /// defined by its adjacency matrix, using the Warshall's algorithm:
03872 /// S. Warshall, A theorem on Boolean matrices, J. ACM 9 (1962) 11-12.
03873 template <class matrix>
03874 inline void transitiveClosure (matrix &m, int_t n)
03875 {
03876         for (int_t k = 0; k < n; ++ k)
03877         {
03878                 for (int_t i = 0; i < n; ++ i)
03879                 {
03880                         for (int_t j = 0; j < n; ++ j)
03881                         {
03882                                 if (m [i] [k] && m [k] [j])
03883                                         m [i] [j] = 1;
03884                         }
03885                 }
03886         }
03887         return;
03888 } /* transitiveClosure */
03889 
03890 /// Computes the transitive reduction of a CLOSED acyclic graph
03891 /// defined by its adjacency matrix, using the algorithm by D. Gries,
03892 /// A.J. Martin, J.L.A. van de Snepscheut and J.T. Udding, An algorithm
03893 /// for transitive reduction of an acyclic graph, Science of Computer
03894 /// Programming 12 (1989), 151-155.
03895 /// WARNING: The input graph MUST BE CLOSED, use the "transitiveClosure"
03896 /// algorithm first if this is not the case.
03897 template <class matrix>
03898 inline void transitiveReduction (matrix &m, int_t n)
03899 {
03900         for (int_t k = n - 1; k >= 0; -- k)
03901         {
03902                 for (int_t i = 0; i < n; ++ i)
03903                 {
03904                         for (int_t j = 0; j < n; ++ j)
03905                         {
03906                                 if (m [i] [k] && m [k] [j])
03907                                         m [i] [j] = 0;
03908                         }
03909                 }
03910         }
03911         return;
03912 } /* transitiveReduction */
03913 
03914 /// Computes the transitive reduction of an arbitrary acyclic graph.
03915 /// The output graph must be initially empty.
03916 template <class wType>
03917 inline void transitiveReduction (const diGraph<wType> &g,
03918         diGraph<wType> &gRed)
03919 {
03920         int_t nVert = g. countVertices ();
03921         if (nVert <= 0)
03922                 return;
03923         flatMatrix<char> m (nVert);
03924         m. clear (0);
03925         graph2matrix (g, m);
03926         transitiveClosure (m, nVert);
03927         transitiveReduction (m, nVert);
03928         matrix2graph (m, nVert, gRed);
03929         return;
03930 } /* transitiveReduction */
03931 
03932 
03933 } // namespace homology
03934 } // namespace chomp
03935 
03936 #endif // _CHOMP_STRUCT_DIGRAPH_H_
03937 
03938 /// @}
03939 

Generated on Sun Feb 3 2013 12:40:30 for The Uniform Expansion Software by  doxygen 1.7.2