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

Generated on Wed Nov 21 11:08:41 2007 for The Uniform Expansion Software by  doxygen 1.5.3