From 20a385c2ea4d1b87aae468b8828b52ccd6ebd34d Mon Sep 17 00:00:00 2001 From: Bradley Meyer Date: Sun, 23 Aug 2015 11:43:45 -0400 Subject: [PATCH] Add rank spanning branchings routine. The development history is in the messages below: Add proposed new files. Use std::less for default comparison. Use std::greater in example. (Sat Aug 29 17:18:27 2015 -0400) Allow for other vertex list type than vecS. Change test code to use setS for vertices. (Wed Sep 2 18:10:41 2015 -0400) Move rank routine to implementation. (Fri Sep 4 18:56:40 2015 -0400) Remove unnecessary template parameters--already defined as global types. (Mon Sep 7 08:32:53 2015 -0400) Move BranchingGraph and BranchingVertex into detail namespace. (Thu Sep 10 13:27:29 2015 -0400) Change to have functor iterate on edges. (Mon Sep 28 16:09:20 2015 -0400) Fix template typename issues. (Mon Sep 28 16:35:14 2015 -0400) Rename files. Generalize input of graph in rank-branchings2.cpp. (Wed Sep 30 17:45:04 2015 -0400) Do some clean up. (Fri Oct 2 09:59:43 2015 -0400) Simplify number of branchings to compute by computing it through the called functor. (Mon Oct 5 13:54:32 2015 -0400) Update test routines. (Mon Oct 5 14:50:48 2015 -0400) Use iterators on edges. (Thu Oct 8 07:45:25 2015 -0400) Add return value. (Thu Oct 8 13:51:29 2015 -0400) Use pair for edge iterators in user-supplied functor. (Sat Oct 10 09:29:19 2015 -0400) Use named parameters. (Sun Nov 15 10:00:42 2015 -0500) Use distance_compare since this compares both edge and branching weights. Restore named_functions. Use named parameters in test_real_ranking (for testing). (Sun Nov 15 13:40:44 2015 -0500) Explicitly call named parameters to test. Add removal for Cygwin in test. Move test for empty graph. (Mon Nov 16 07:49:44 2015 -0500) Add documentation. Alphabetize includes. (Mon Nov 16 17:00:42 2015 -0500) Update. (Wed Nov 18 18:37:29 2015 -0500) Update documentation and add references. (Sun Nov 22 10:26:58 2015 -0500) Remove unnecessary concept checks. Update documentation. (Sun Nov 22 12:40:32 2015 -0500) Update documentation. (Mon Nov 23 08:18:57 2015 -0500) Update references. Add example usage statement. (Tue Nov 24 07:50:30 2015 -0500) Add sentence on weight modification. (Thu Nov 26 08:05:46 2015 -0500) Simplify removal of parallel arcs. (Fri Dec 11 15:05:50 2015 -0500) Concept check that edge weight is less than comparable only if edge comparison function is not supplied (to allow, for example, complex edge weights if user supplies the appropriate comparison function). (Sun Feb 28 11:26:30 2016 -0500) Do some clean up. (Tue Mar 1 18:53:57 2016 -0500) Update jam file. Move tests to single code. (Fri Apr 8 15:12:34 2016 -0400) Add test of graphs with complex weight. (Sat Apr 9 09:47:51 2016 -0400) Remove unused typedef Vertex and place check on typedef Edge. (Wed Apr 20 16:11:14 2016 -0400) Simplify initialization. (Mon Aug 8 10:09:40 2016 -0400) Change to use filtered graph for branching processor. Explicitly add branching edge to in edges in condensed vertex in next to correct problem when there are many in edges with the same weight such that the branching edge might be skipped. (Thu Sep 29 17:48:39 2016 -0400) Add examples to Jamfile. Use scoped pointer for summing branching weights to allow general weight type. Add class weight to test. (Tue Oct 4 17:07:43 2016 -0400) Remove extraneous comment. (Sat Nov 5 10:41:54 2016 -0400) --- doc/bibliography.html | 30 + doc/rank_spanning_branchings.html | 381 +++++++ example/Jamfile.v2 | 2 + example/branching_input.txt | 16 + example/rank-branchings1.cpp | 110 ++ example/rank-branchings2.cpp | 233 +++++ include/boost/graph/rank_spanning_branchings.hpp | 1159 ++++++++++++++++++++++ test/Jamfile.v2 | 2 + test/rank_branchings_test.cpp | 835 ++++++++++++++++ 9 files changed, 2768 insertions(+) create mode 100644 doc/rank_spanning_branchings.html create mode 100644 example/branching_input.txt create mode 100644 example/rank-branchings1.cpp create mode 100644 example/rank-branchings2.cpp create mode 100644 include/boost/graph/rank_spanning_branchings.hpp create mode 100644 test/rank_branchings_test.cpp diff --git a/doc/bibliography.html b/doc/bibliography.html index 7b5bca55..24aac759 100644 --- a/doc/bibliography.html +++ b/doc/bibliography.html @@ -438,6 +438,36 @@
Generating random spanning trees more quickly than the cover time. ACM Symposium on the Theory of Computing, pp. 296-303, 1996. +

74 +
+P. M. Camerini, L. Fratta, and F. Maffioli +
The K Best Spanning Arborescences of a Network
+Networks 10: 91-110, 1980. + +

75 +
+Y. J. Chu and T. H. Liu +
On the Shortest Arborescence of a Directed Graph
+Science Sinica 14: 1396-1400, 1965. + +

76 +
+J. Edmonds +
Optimum Branchings
+J. Res. Nat. Bur. Standards 71B: 233-240, 1967. + +

77 +
+R. E. Tarjan +
Finding Optimum Branchings
+Networks 7: 25-35, 1977. + +

78 +
+P. M. Camerini, L. Fratta, and F. Maffioli +
A Note on Finding Optimum Branchings
+Networks 9: 309-312, 1979. +
diff --git a/doc/rank_spanning_branchings.html b/doc/rank_spanning_branchings.html new file mode 100644 index 00000000..11453a07 --- /dev/null +++ b/doc/rank_spanning_branchings.html @@ -0,0 +1,381 @@ + + + +Boost Graph Library: Rank Spanning Branchings + +C++ Boost + +
+ + +

+rank_spanning_branchings +

+ +
+template <class Graph, class BranchingProcessor, class P, class T, class R>
+void
+rank_spanning_branchings(Graph& g, BranchingProcessor bp, 
+    const bgl_named_params<P, T, R>& params = all defaults);
+
+ +

+The rank_spanning_branchings() function finds, in order by branching +weight, the spanning branchings in a directed graph with weighted edges. +If the order of a directed graph is N, a spanning branching is an +acyclic subgraph of +that graph with N-1 edges such that no vertex has indegree larger than +one. The weight of the branching is the sum of the weights of the edges +in the branching. +This function uses the Camerini et al. algorithm +[74] +to compute the +branchings in order by branching weight. Once the routine finds a branching, +the routine calls a user-supplied +BranchingProcessor +to process the branching and then +moves on to find the next branching. The +user-supplied +BranchingProcessor +must take as input +a filtered view of the parent graph +that gives the branching. +

+

+The Camerini et al. algorithm starts with a directed graph that has +a single root vertex (that is, a vertex that has indegree equal to zero +and outdegree not equal to zero). The procedure RANK +finds from the procedure BEST the spanning +branching A with the largest +branching weight, denoted d(A). +The +BranchingProcessor +is invoked on A, and if it returns true, the procedure continues. +The procedure NEXT is then invoked to find +the edge e such that when e is removed from G, +BEST will find the next largest branching. +NEXT also returns δ +such that d(A)-δ is the weight of the next largest branching. +

+ +

+RANK keeps a priority queue of tuples with elements (weight, edge, +branching, Y, Z). The element +branching is a set of edges giving a branching. Y is a set of +edges that must be included in a branching. +Z is a set of edges that must be excluded from a branching. After +each evaluation of NEXT, +the appropriate tuple is inserted into Q. +The priority queue ranks tuples by the weight entry. +Once the maximum weight branching is found, RANK finds subsequent +branchings by popping a tuple off the top of Q and running +BEST. Possible descendants are then evaluated with NEXT. +

+ +
+RANK(G, BranchingProcessor())
+  A := BEST(G, Ø, Ø)
+  if (!BranchingProcessor(A)) return
+  (e,δ) := NEXT(G, A, Ø, Ø)
+  INSERT( Q, (d(A) - δ, e, A, Ø, Ø) )
+  while (Q != Ø)
+    (w, e, A, Y, Z) := POP(Q)
+    Y' := Y U e
+    Z' := Z U e
+    A' := BEST(G, Y, Z')
+    if (!BranchingProcessor(A')) return
+    (e',δ) := NEXT(G, A, Y', Z)
+    INSERT( Q, (d(A) - δ, e', A, Y', Z))
+    (e',δ) := NEXT(G, A', Y, Z')
+    INSERT( Q, (w - δ, e', A', Y, Z'))
+  end while 
+  return
+
+ +

+The procedure BEST is the Chu-Liu/Edmond's algorithm +[75, +76] +in the form adopted by +Tarjan [77] and +Camerini et al. +[78]. +The procedure begins by +constructing a subgraph H of G. To construct H, +the procedure loops over each edge in Y and removes from the graph +every other edge with the same target. This ensures that the edges in +Y will be part of the branching found in this invocation of +BEST. The procedure then removes the edges in Z from the +graph; thus, the edges in Z will not be part of the branching. +BEST finds the largest edge into each vertex. When a cycle forms, +the cycle is condensed to a new vertex (such that the graph becomes +HC) and the weight of +edges into the cycle (that is, into the new vertex) +are modified. The modification is such that, if v is a vertex +in the cycle and is the invertex for edge e whose source is outside +the cycle, the weight of e is modified +by adding the weight of the least costly edge within the cycle and +subtracting the weight of the cycle edge into v. +This proceeds until all +vertices have been examined. The condensed graph is then expanded to give the +best branching given the constraints Y and Z. +

+ +
+BEST(G, Y, Z)
+  H := GYZ
+  B := Ø
+  Γ := (V,Ø)
+  while there exists an exposed vertex v ≠ root in H with respect to B
+    let wH be the weighting function of H
+    let Q be a priority queue of in edges of v ordered by wH
+    b := TOP(Q)
+    B := B U b
+    β(v) := b
+    if B contains a cycle C of H then
+      H := HC
+      let u be the vertex of HC corresponding to C
+      add u to Γ
+      for each vertex v' in C
+        make v' a child of u in Γ
+      end for
+      B := B - C
+    end if
+  end while
+  while Γ contains non-isolated vertices
+    find path in Γ with vertices v1, v2, ... vk such that v1 is any non-isolated root of Γ and vk = t(β(v1))
+    for h = 1, k - 1
+      β(vh+1) := β(vh)
+      remove from Γ vertex vh and all edges directed out of vh
+    end for
+  end while
+  A := {β(v):v ≠ root is a vertex of G}
+  return A
+
+ +

+The procedure NEXT constructs the constrained graph H as +in BEST. It examines each edge in the branching to determine +whether swapping an edge not in the input branching for an edge in the +branching results in a smaller +weight difference. Only those alternative +edges for which the target of the +branching edge (to be replaced) is not an ancestor of the source of the +alternative edge in the branching are allowed as possible +swapped edges (as determined by procedure SEEK). The implementation +of rank_spanning_branchings() does not in fact use a separate +SEEK routine but rather simply an ordered iteration over the +priority queue containing in edges to the vertex under +examination. Cycles are handled as in BEST. NEXT +returns the edge e that when removed from the input graph to +BEST will +give the next best branching and the weight difference between the current +best branching and the next best descendent. +

+ +
+NEXT(G, A, Y, Z)
+  H := GYZ
+  B := Ø
+  δ := ∞
+  while there exists an exposed vertex v ≠ root in H with respect to B
+    let wH be the weighting function of H
+    let Q be a priority queue of in edges of v ordered by wH
+    b := TOP(Q) (ties among edges must be solved in favor of edges in A)
+    B := B U b
+    if b ∈ A - Y 
+       w' := SEEK(b, Q, A, H)
+       if wH(b) - w' < δ
+         e := b
+         δ := wH(b) - w'
+       end if
+    end if
+    if B contains a cycle C of H
+      H := HC
+      B := B - C
+    end if
+  end while
+  return (e, δ)
+
+ +
+SEEK(b, Q, A, H )
+  for f in Q (ordered iteration by wH)
+    if f != b and target of b is not an ancestor of source of f in A
+      return wH(f)
+    end if
+  end for
+  return -∞
+
+ +

Where Defined

+ +

+boost/graph/rank_spanning_branchings.hpp + +

+ +

+BranchingProcessor +

+ +

+BranchingProcessor is used in the rank_spanning_branchings() +function to process the current branching. It is a functor that is called with +a filtered view of the parent graph. +The processor must return a boolean. If the return value is +true, rank_spanning_branchings() +continues and seeks the next branching. +If the return value is false, rank_spanning_branchings() stops. +

+ +

+The following functor is an example +BranchingProcessor that prints out the edges in a branching. +

+struct print_branching
+{
+
+  print_branching() {} 
+
+  template<class BranchingGraph>
+  bool operator()( const BranchingGraph& bg )
+  {
+    std::cout << "Branching:";
+    BGL_FORALL_EDGES_T( e, bg, BranchingGraph )
+    {
+      std::cout << " (" << source( e, bg ) << "," << target( e, bg ) << ")";
+    }
+    std::cout << std::endl;
+    return true;
+  }
+};
+
+It returns true, so rank_spanning_branchings() would always continue +on to find the next branching. +

+ +

+This example +BranchingProcessor could be called on a Graph g as +

+rank_spanning_branchings( g, print_branching() );
+
+This would print out all spanning branchings in the graph in descending order +by branching weight. +

+ +

Parameters

+ +IN: const Graph& g +
+ A directed graph. The graph type must be a model of + VertexAndEdgeListGraph. + The graph should + have a single root vertex, that is, a single vertex that has indegree + equal to zero and an edge directed to at least one other vertex. + No other vertex should have indegree equal to zero. + A graph with at least one spanning branching + can always be constructed from a general + directed graph by adding a vertex (the root vertex) and then adding edges + with weight zero directed from the root vertex to every other vertex.
+
+ + + IN: BranchingProcessor bp +
+ A functor that models BranchingProcessor to process the edges in a branching.
+
+ +

Named Parameters

+ +IN: weight_map(WeightMap w_map) +
+The weight or ``length'' of + each edge in the graph. + The WeightMap type must be a model + of Readable + Property Map. If no edge weight comparison function is supplied, + its value type must be Less Than + Comparable. The key type of this map needs to be the graph's + edge descriptor type.
+ Default: get(edge_weight, g)
+
+ +IN: vertex_index_map(VertexIndexMap i_map) +
+ This maps each vertex to an integer in the range [0, + num_vertices(g)). + The type VertexIndexMap + must be a model of Readable Property + Map. The value type of the map must be an integer type. The + vertex descriptor type of the graph needs to be usable as the key + type of the map.
+ Default: get(vertex_index, g) + Note: if you use this default, make sure your graph has + an internal vertex_index property. For example, + adjacenty_list with VertexList=listS does + not have an internal vertex_index property. +
+
+ +IN: distance_compare(CompareFunction cmp) +
+ This function is used to compare edge weights to determine which + edges to include + in the next best branching. It is also used to compare the weights of + branchings. The weight of a branching is the sum of weights of its edges. + The CompareFunction type must be a model of Binary + Predicate and have argument types that match the value type of + the WeightMap property map.
+ + Default: + std::less<W> with W=typename + property_traits<WeightMap>::value_type
+
+ +

Complexity

+ +

+The time complexity is O(kE log V), where k is the number +of spanning branchings found. + +

Example

+ +

+The files example/rank-branching1.cpp +and +example/rank-branching2.cpp +contain examples of using the Camerini et al. algorithm. + +

Acknowledgments

+ +

+The development of rank_spanning_branchings() was inspired by +the C++ implementation of +Edmond's algorithm by Ali Tofigh and Erik Sjölund. +

+ +
+
+ + +
Copyright © 2015 +Clemson University, Bradley S. Meyer (mbradle@clemson.edu) +
+ + + diff --git a/example/Jamfile.v2 b/example/Jamfile.v2 index 13cdf75f..fa8aa346 100644 --- a/example/Jamfile.v2 +++ b/example/Jamfile.v2 @@ -54,4 +54,6 @@ exe hawick_circuits : hawick_circuits.cpp ; exe edge_coloring : edge_coloring.cpp ; exe successive_shortest_path_nonnegative_weights_example : successive_shortest_path_nonnegative_weights_example.cpp ; exe cycle_canceling_example : cycle_canceling_example.cpp ; +exe rank-branchings1 : rank-branchings1.cpp ; +exe rank-branchings2 : rank-branchings2.cpp ; diff --git a/example/branching_input.txt b/example/branching_input.txt new file mode 100644 index 00000000..2d6f336f --- /dev/null +++ b/example/branching_input.txt @@ -0,0 +1,16 @@ +v 1, v 2, 1 +v 1, v 3, 2 +v 1, v 4, -3 +v 2, v 1, -2 +v 2, v 3, 6 +v 2, v 4, 4 +v 3, v 1, 1 +v 3, v 2, 5 +v 3, v 4, 3 +v 4, v 1, -2 +v 4, v 2, 7 +v 4, v 3, 1 +root, v 1, 0 +root, v 2, 0 +root, v 3, 0 +root, v 4, 0 diff --git a/example/rank-branchings1.cpp b/example/rank-branchings1.cpp new file mode 100644 index 00000000..4281af55 --- /dev/null +++ b/example/rank-branchings1.cpp @@ -0,0 +1,110 @@ +// Copyright 2015 Clemson University +// Authors: Bradley S. Meyer +// +// Distributed under the Boost Software License, Version 1.0. (See +// accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) +//======================================================================= +#include +#include +#include +#include + +struct print_branching +{ + + print_branching(){} + + template + bool operator()( BranchingGraph& bg ) + { + + typedef + typename + boost::property_map::const_type + WeightMap; + + WeightMap w; + typename boost::property_traits::value_type weight; + + std::cout << "Branching:"; + + weight = 0.; + + BGL_FORALL_EDGES_T( e, bg, BranchingGraph ) + { + + std::cout << " (" << boost::source( e, bg ) << "," << + boost::target( e, bg ) << ")"; + + weight += get( w, e ); + + } + + std::cout << std::endl << " Weight: " << weight << std::endl << std::endl; + + return true; + + } + +}; + +int +main() +{ + typedef boost::adjacency_list < boost::vecS, boost::vecS, boost::directedS, + boost::no_property, boost::property < boost::edge_weight_t, int > > Graph; +#if defined(BOOST_MSVC) && BOOST_MSVC <= 1300 + typedef boost::graph_traits < Graph >::edge_descriptor Edge; +#endif + typedef std::pair E; + + const int num_nodes = 4; + E edge_array[] = { E(0, 1), E(0, 2), E(0, 3), E(1, 2), E(2, 1), E(2, 3), + E(3, 2), E(1, 3), E(3, 1) + }; + int weights[] = { 5, 1, 1, 11, 10, 5, 8, 4, 9 }; + + std::size_t num_edges = sizeof(edge_array) / sizeof(E); +#if defined(BOOST_MSVC) && BOOST_MSVC <= 1300 + Graph g(num_nodes); + boost::property_map::type weightmap = + get(edge_weight, g); + for (std::size_t j = 0; j < num_edges; ++j) { + Edge e; bool inserted; + boost::tie(e, inserted) = add_edge(edge_array[j].first, edge_array[j].second, g); + weightmap[e] = weights[j]; + } +#else + Graph g(edge_array, edge_array + num_edges, weights, num_nodes); +#endif + + // Rank branchings in descending order of weight (use default + // weight comparison). + + std::cout << std::endl; + + std::cout << "Spanning branchings in descending order of weight:" + << std::endl << std::endl; + + boost::rank_spanning_branchings( + g, + print_branching() + ); + + std::cout << std::endl; + + // Rank branchings in ascending order of weight (use supplied + // weight comparison). + + std::cout << "Spanning branchings in ascending order of weight:" + << std::endl << std::endl; + + boost::rank_spanning_branchings( + g, + print_branching(), + boost::distance_compare( std::greater() ) + ); + + return EXIT_SUCCESS; +} diff --git a/example/rank-branchings2.cpp b/example/rank-branchings2.cpp new file mode 100644 index 00000000..6efa163f --- /dev/null +++ b/example/rank-branchings2.cpp @@ -0,0 +1,233 @@ +//======================================================================= +// Copyright 2015 Clemson University +// Authors: Bradley S. Meyer +// +// Distributed under the Boost Software License, Version 1.0. (See +// accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) +//======================================================================= + +/* + Example execution: + + ./rank-branchings2 branching_input.txt -10 + +*/ + +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +// When rank_spanning_branchings finds a branching, it calls the user-defined +// functor (in this case, my_function). The functor takes as input +// a filtered view of the parent graph that gives the branching. If the +// functor returns false, rank_spanning_branchings will stop ranking the +// branchings and will return. If the functor returns true, +// rank_spanning_branchings will continue to find the next branching. + +using namespace boost; + +template +class my_function +{ + + public: + + typedef + typename property_map::const_type WeightMap; + + typedef + typename property_traits::value_type weight_t; + + my_function( weight_t& _max_weight, weight_t _cut ) : + max_weight( _max_weight ), cut( _cut ) + {} + + template + bool operator()( BranchingGraph& bg ) + { + + WeightMap w; + + typename property_map::const_type name_map; + + weight_t weight = 0; + + b_string = ""; + + BGL_FORALL_EDGES_T( e, bg, BranchingGraph ) + { + + weight += get( w, e ); + + std::stringstream ss; + ss << "(" << name_map[source( e, bg )] << ", " << + name_map[target( e, bg )] << ") "; + + b_string += ss.str(); + } + + if( max_weight == -std::numeric_limits::infinity() ) + { + max_weight = weight; + } + + d_diff = weight - max_weight; + + if( d_diff < cut ) + { + std::cout << std::endl; + return false; + } // Stop before output. + + std::cout << "Branching: " << b_string << std::endl; + + std::cout << " Weight = " << weight << std::endl; + + std::cout << " Weight - Max Weight = " << d_diff << std::endl; + + std::cout << std::endl; + + return true; + + } + + private: + weight_t& max_weight; + weight_t cut; + std::string b_string; + weight_t d_diff; + +}; + +// Check for vertex from input file and add if not already present. + +template +void +check_vertex( Graph& g, std::map& s_map, std::string str ) +{ + + typedef typename property_map < Graph, vertex_index_t >::type index_map_t; + typedef typename property_map < Graph, vertex_name_t >::type name_map_t; + + index_map_t index_map = get( vertex_index, g ); + name_map_t name_map = get( vertex_name, g ); + + if( s_map.find( str ) == s_map.end() ) + { + Vertex u = add_vertex( g ); + index_map[u] = s_map.size(); + put( name_map, u, str ); + s_map[str] = u; + } + +} + +// Read in the graph from the file. Each line in the input file should be +// a comma-delimited triplet (source, target, weight). + +template +void +read_graph_file( + std::istream & graph_in, + Graph & g +) +{ + + typedef typename property_map::const_type WeightMap; + + typedef typename property_traits::value_type weight_t; + + typedef typename graph_traits::vertex_descriptor Vertex; + std::map s_map; + std::string text; + + char_separator sep(","); + + while( std::getline( graph_in, text ) ) + { + + std::vector str; + tokenizer > tokens( text, sep ); + + BOOST_FOREACH( const std::string& name, tokens ) + { + std::string s = name; + algorithm::trim( s ); + str.push_back( s ); + } + + check_vertex( g, s_map, str[0] ); // Check source; + check_vertex( g, s_map, str[1] ); // Check target; + + add_edge( + s_map[str[0]], s_map[str[1]], lexical_cast( str[2] ), g + ); + + } + +} + +int main( int argc, char **argv ) +{ + + typedef int my_type; + + typedef adjacency_list < + listS, listS, directedS, + property < + vertex_index_t, size_t, + property< vertex_name_t, std::string > + >, + property < edge_weight_t, my_type> + > Graph; + std::ifstream input_file; + my_type max_weight = -std::numeric_limits::infinity(); + + if( argc != 3 ) + { + std::cerr << std::endl; + std::cerr << " Usage: " << argv[0] << + " file cut" << std::endl << std::endl; + std::cerr << " file = input file" << std::endl; + std::cerr << " cut = branching weight cut" << std::endl << std::endl; + std::cerr << " Example usage: " << argv[0] << + " branching_input.txt -10" << std::endl << std::endl;; + return EXIT_FAILURE; + } + + Graph g; + + input_file.open( argv[1] ); + + if( !input_file.is_open() ) + { + std::cerr << "Invalid input file." << std::endl; + return EXIT_FAILURE; + } + + std::ifstream file_in( argv[1] ); + + read_graph_file( file_in, g ); + + std::cout << std::endl; + + rank_spanning_branchings( + g, + my_function( + max_weight, + lexical_cast( argv[2] ) + ) + ); + + return EXIT_SUCCESS; + +} diff --git a/include/boost/graph/rank_spanning_branchings.hpp b/include/boost/graph/rank_spanning_branchings.hpp new file mode 100644 index 00000000..a4640597 --- /dev/null +++ b/include/boost/graph/rank_spanning_branchings.hpp @@ -0,0 +1,1159 @@ +//======================================================================= +// Copyright 2015 Clemson University +// Authors: Bradley S. Meyer +// +// Distributed under the Boost Software License, Version 1.0. (See +// accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) +//======================================================================= +#ifndef BOOST_GRAPH_RANK_SPANNING_BRANCHINGS_HPP +#define BOOST_GRAPH_RANK_SPANNING_BRANCHINGS_HPP + +/* + * Rank spanning branchings + * Camerini et al Algorithm + * + * Requirement: + * directed graph with single root vertex + */ + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace boost { + + namespace detail { + + typedef adjacency_list < vecS, vecS, bidirectionalS, + no_property, no_property > BranchingGraph; + + typedef graph_traits::vertex_descriptor BranchingVertex; + + template + struct EdgeNode + { + Edge edge; + typename property_traits::value_type weight; + Compare compare; + EdgeNode(){} + EdgeNode( + const Edge& e, + const typename property_traits::value_type & w, + Compare c ) : + edge( e ), weight( w ), compare( c ) {} + bool operator<( EdgeNode const & rhs ) const + { return compare( weight, rhs.weight ); } + }; + + // Insert edges from graph into heap, taking into account constraints. + + template + bool + insert_edges + ( + const Graph& g, + WeightMap& w, + Compare& comp, + std::map< + Vertex, + heap::fibonacci_heap > + >& in_edges, + const unordered_set& include_edges, + const unordered_set& exclude_edges + ) + { + unordered_set include_vertices, in_set, out_set; + typename unordered_set::iterator it; + + // Insert vertices in sets to check for spanning branching. + + BGL_FORALL_VERTICES_T( v, g, Graph ) + { + in_set.insert( v ); + out_set.insert( v ); + } + + // Insert edges that must be present and note the invertex. Remove + // vertices from relevant sets. + + BOOST_FOREACH( const Edge& e, include_edges ) + { + include_vertices.insert( target( e, g ) ); + in_edges[target( e, g )].push( + EdgeNode( e, get( w, e ), comp ) + ); + it = in_set.find( target( e, g ) ); + if( it != in_set.end() ) { in_set.erase( it ); } + it = out_set.find( source( e, g ) ); + if( it != out_set.end() ) { out_set.erase( it ); } + } + + // Insert edges, but not edges to be excluded or into vertices that + // already have in edges. Remove vertices from relevant sets. + + BOOST_FOREACH( const Edge &e, edges(g) ) + { + if( + include_vertices.find( target(e, g) ) == include_vertices.end() + && + exclude_edges.find( e ) == exclude_edges.end() + ) + { + if( source(e, g) != target(e, g) ) + { + in_edges[target(e,g)].push( + EdgeNode( e, get( w, e ), comp ) + ); + it = in_set.find( target( e, g ) ); + if( it != in_set.end() ) { in_set.erase( it ); } + it = out_set.find( source( e, g ) ); + if( it != out_set.end() ) { out_set.erase( it ); } + } + } + } + + // Check for correct number of roots. + // + if( in_set.size() != 1 ) + return false; // Zero roots or more than one root. + else if( out_set.find( *in_set.begin() ) != out_set.end() ) + return false; // Root is isolated. + else + return true; + + } + + // Retrieve the path back from leaf to root in cycle branching. + + void + find_back_path( + BranchingGraph& C, + std::vector& bv + ) + { + + BGL_FORALL_INEDGES( bv[0], e, C, BranchingGraph ) + { + bv.insert( bv.begin(), source( e, C ) ); + find_back_path( C, bv ); + } + + } + + // Expand cycles. + + template + void + expand( + const Graph& g, + IndexMap& v_id, + BranchingGraph& C, + unordered_set& root_set, + std::vector >& beta + ) + { + + BOOST_FOREACH( BranchingVertex v, root_set ) + { + if( in_degree( v, C ) == 0 && out_degree( v, C ) != 0 ) + { + std::vector bv; + bv.push_back( v_id[target( beta[v].edge, g )] ); + find_back_path( C, bv ); + for( std::size_t i = 0; i < bv.size() - 1; i++ ) + { + beta[bv[i+1]] = beta[bv[i]]; + clear_vertex( bv[i], C ); + } + } + } + + // Remove isolated vertices. + + std::vector vertices_to_remove_from_set; + + BOOST_FOREACH( BranchingVertex v, root_set ) + { + if( in_degree( v, C ) == 0 && out_degree( v, C ) == 0 ) + { + vertices_to_remove_from_set.push_back( v ); + } + } + + BOOST_FOREACH( BranchingVertex v, vertices_to_remove_from_set ) + { + root_set.erase( v ); + } + + } + + // Camerini et al. BEST routine. + + template + void + best_spanning_branching( const Graph& g, + unordered_set& branching, + IndexMap& v_id, + WeightMap& w, + Compare& comp, + Rank rank, + Pred pred1, + Pred pred2, + unordered_set& include_edges, + unordered_set& exclude_edges + ) + { + + typedef typename graph_traits::vertex_descriptor Vertex; + typedef typename graph_traits::vertices_size_type vertex_idx_t; + + typedef EdgeNode edge_node_t; + + typedef heap::fibonacci_heap f_heap_t; + typedef std::map exit_map_t; + + unordered_set unvisited_vertex_set; + + Vertex root_vertex; + + BranchingGraph C; + + vertex_idx_t n = num_vertices(g); + + disjoint_sets W( rank, pred1 ), S( rank, pred2 ); + + std::map in_edges; + + if( + !insert_edges + ( + g, + w, + comp, + in_edges, + include_edges, + exclude_edges + ) + ) return; + + std::vector beta( 2 * n ); + + std::map parent; + + BGL_FORALL_VERTICES_T( v, g, Graph ) + { + W.make_set( v ); + S.make_set( v ); + parent[v] = v_id[v]; + add_vertex( C ); + unvisited_vertex_set.insert( v ); + } + + while( !unvisited_vertex_set.empty() ) + { + + typename unordered_set::iterator it = + unvisited_vertex_set.begin(); + + if( in_edges[*it].empty() ) + { + root_vertex = *it; + unvisited_vertex_set.erase( it ); + } + else + { + + edge_node_t critical_edge_node = in_edges[*it].top(); + + beta[parent[*it]] = critical_edge_node; + + unvisited_vertex_set.erase( it ); + + if( + W.find_set( source( critical_edge_node.edge, g ) ) != + W.find_set( target( critical_edge_node.edge, g ) ) + ) + { + W.union_set( + source( critical_edge_node.edge, g ), + target( critical_edge_node.edge, g ) + ); + } + else + { + BranchingVertex v_new = add_vertex( C ); + + edge_node_t least_costly_edge_node = critical_edge_node; + + boost::unordered_set cycle_vertex_set; + + for( + Vertex v = source( critical_edge_node.edge, g ); + cycle_vertex_set.find( S.find_set( v ) ) == + cycle_vertex_set.end(); + v = source( beta[parent[S.find_set( v )]].edge, g ) + ) + { + Vertex u = S.find_set( v ); + cycle_vertex_set.insert( u ); + add_edge( v_new, parent[u], C ); + if( + comp( beta[parent[u]].weight, least_costly_edge_node.weight ) + ) + { + least_costly_edge_node = beta[parent[u]]; + } + } + + Vertex new_repr = *cycle_vertex_set.begin(); + + BOOST_FOREACH( Vertex u, cycle_vertex_set ) + { + S.link( u, new_repr ); + new_repr = S.find_set( new_repr ); + } + + BOOST_FOREACH( Vertex v, cycle_vertex_set ) + { + exit_map_t v_exit; + BOOST_FOREACH( edge_node_t en, in_edges[v] ) + { + if( S.find_set( source( en.edge, g ) ) != new_repr ) + { + en.weight += least_costly_edge_node.weight - + beta[parent[v]].weight; + Vertex u = S.find_set( source( en.edge, g ) ); + if( v_exit.find(u) != v_exit.end() ) + { + if( comp( v_exit[u].weight, en.weight ) ) + { + v_exit[u] = en; + } + } + else + { + v_exit[u] = en; + } + } + } + f_heap_t tmp_heap; + BOOST_FOREACH( typename exit_map_t::value_type& t, v_exit ) + { + tmp_heap.push( t.second ); + } + in_edges[v].swap( tmp_heap ); + } + + BOOST_FOREACH( Vertex v, cycle_vertex_set ) + { + if( v != new_repr ) in_edges[new_repr].merge( in_edges[v] ); + } + + unvisited_vertex_set.insert( new_repr ); + + parent[new_repr] = v_new; + + } + + } + + } + + // Create a set containing possible roots of the cycle branching. + // In each cycle expansion, remove isolated vertices of C to avoid + // considering them in subsequent cycle expansions. + + unordered_set root_set; + + BGL_FORALL_VERTICES( u, C, BranchingGraph ) + { + root_set.insert( u ); + } + + while( !root_set.empty() ) + { + expand( g, v_id, C, root_set, beta ); + } + + BGL_FORALL_VERTICES_T( v, g, Graph ) + { + if( v != root_vertex ) + { + branching.insert( beta[v_id[v]].edge ); + } + } + + } + + // Depth-first visitor to set up pre-order and post-order maps. + + template + class dfs_order_visitor:public default_dfs_visitor + { + + public: + dfs_order_visitor( + OrderMap& pr, + OrderMap& po + ) : m_pr( pr ), m_po( po ) { td = 0; tf = 0; } + + template + void discover_vertex(const Vertex u, const Graph & g) + { + m_pr[u] = td++; + } + + template + void finish_vertex(const Vertex u, const Graph & g) + { + m_po[u] = tf++; + } + + OrderMap& m_pr; + OrderMap& m_po; + + private: + std::size_t td; + std::size_t tf; + + }; + + // The ancestor checker. A vertex u is a proper ancestor of + // vertex v (in the parent branching) if pr[u] < pr[v] and po[u] > po[v]. + + template + struct ancestor_checker + { + OrderMap& m_pr; + OrderMap& m_po; + ancestor_checker( OrderMap& pr, OrderMap& po ) : + m_pr( pr ), m_po( po ) {} + + template + bool operator()( const Vertex& v1, const Vertex& v2 ) const + { return ( m_pr[v1] < m_pr[v2] && m_po[v1] > m_po[v2] ); } + }; + + // Camerini et al. NEXT routine. + + template + shared_ptr< + std::pair::value_type> + > + next_spanning_branching( const Graph& g, + const unordered_set& branching, + IndexMap& v_id, + WeightMap& w, + Compare& comp, + Rank rank, + Pred pred1, + Pred pred2, + unordered_set& include_edges, + unordered_set& exclude_edges + ) + { + + typedef typename graph_traits::vertex_descriptor Vertex; + + Edge return_edge; + + typedef typename property_traits::value_type weight_t; + + scoped_ptr delta; + + typedef EdgeNode edge_node_t; + + edge_node_t b; + + // Initialize delta. + + delta.reset(); + + // Create a branching graph to check whether a vertex is an ancestor of + // another vertex in the branching. + + BranchingGraph C; + + BGL_FORALL_VERTICES_T( v, g, Graph ) + { + add_vertex( C ); + } + + BOOST_FOREACH( const Edge& e, branching ) + { + add_edge( v_id[source( e, g)], v_id[target( e, g )], C ); + } + + // Create other types and data. + + typedef heap::fibonacci_heap f_heap_t; + typedef std::map exit_map_t; + + shared_ptr< + std::pair::value_type> + > p( + new std::pair::value_type> + ); + + unordered_set unvisited_vertex_set; + + disjoint_sets W( rank, pred1 ), S( rank, pred2 ); + + std::map in_edges; + + std::map max_e; + + if( + !insert_edges + ( + g, + w, + comp, + in_edges, + include_edges, + exclude_edges + ) + ) + { + p.reset(); + return p; + } + + // Initialize data structures and find start vertex for depth + // first search. + + BranchingVertex start_vertex; + + BGL_FORALL_VERTICES_T( v, g, Graph ) + { + W.make_set( v ); + S.make_set( v ); + unvisited_vertex_set.insert( v ); + if( in_edges[v].empty() ) { start_vertex = v_id[v]; } + } + + // Create the ancestor checker from C. + + typedef std::map OrderMap; + + OrderMap pr; + OrderMap po; + + dfs_order_visitor vis(pr, po); + + depth_first_search(C, visitor(vis).root_vertex( start_vertex ) ); + + ancestor_checker is_ancestor(pr, po); + + // Main loop. + + while( !unvisited_vertex_set.empty() ) + { + + typename unordered_set::iterator it = + unvisited_vertex_set.begin(); + + if( in_edges[*it].empty() ) + { + unvisited_vertex_set.erase( it ); + } + else + { + + // Get largest in edge with ties solved in favor of edges in + // input branching. + + for( + typename f_heap_t::ordered_iterator ei = + in_edges[*it].ordered_begin(); + ei != in_edges[*it].ordered_end(); + ei++ + ) + { + if( comp( (*ei).weight, in_edges[*it].top().weight ) ) + break; + b = *ei; + if( branching.find( b.edge ) != branching.end() ) + break; + } + + max_e[*it] = b; + + // Camerini et al. SEEK routine. + + if( branching.find( b.edge ) != branching.end() ) + { + for( + typename f_heap_t::ordered_iterator ei = + in_edges[*it].ordered_begin(); + ei != in_edges[*it].ordered_end(); + ei++ + ) + { + if( (*ei).edge != b.edge ) + { + if( + !is_ancestor( + v_id[target( b.edge, g )], + v_id[source( (*ei).edge, g )] + ) + ) + { + if( !delta ) + { + delta.reset( new weight_t ); + *delta = b.weight - (*ei).weight; + return_edge = b.edge; + break; + } + else if( + comp( b.weight - (*ei).weight, *delta ) + ) + { + *delta = b.weight - (*ei).weight; + return_edge = b.edge; + break; + } + } + } + } + } + + unvisited_vertex_set.erase( it ); + + if( + W.find_set( source( b.edge, g ) ) != + W.find_set( target( b.edge, g ) ) + ) + { + W.union_set( source( b.edge, g ), target( b.edge, g ) ); + } + else + { + + edge_node_t least_costly_edge_node = b; + + boost::unordered_set cycle_vertex_set; + + for( + Vertex v = source( b.edge, g ); + cycle_vertex_set.find( S.find_set( v ) ) == + cycle_vertex_set.end(); + v = source( max_e[S.find_set( v )].edge, g ) + ) + { + Vertex u = S.find_set( v ); + cycle_vertex_set.insert( u ); + if( comp( max_e[u].weight, least_costly_edge_node.weight ) ) + { + least_costly_edge_node = max_e[u]; + } + } + + Vertex new_repr = *cycle_vertex_set.begin(); + + BOOST_FOREACH( Vertex v, cycle_vertex_set ) + { + S.link( v, new_repr ); + new_repr = S.find_set( new_repr ); + } + + // Adjust arc weights and remove parallel arcs. Keep the + // the largest weight in arc and the largest viable alternative + // arc from each source outside the cycle. Make sure that + // an arc from the branching is among the added edges, if present. + + BOOST_FOREACH( Vertex v, cycle_vertex_set ) + { + exit_map_t b_exit, v_exit1, v_exit2; + BOOST_FOREACH( edge_node_t en, in_edges[v] ) + { + if( S.find_set( source( en.edge, g ) ) != new_repr ) + { + en.weight += least_costly_edge_node.weight - max_e[v].weight; + Vertex u = S.find_set( source( en.edge, g ) ); + if( branching.find( en.edge ) != branching.end() ) + { + b_exit[u] = en; + } + else + { + if( v_exit1.find(u) != v_exit1.end() ) + { + if( comp( v_exit1[u].weight, en.weight ) ) + { + if( + !is_ancestor( + v_id[target( v_exit1[u].edge, g )], + v_id[source( v_exit1[u].edge, g )] + ) + ) + { + v_exit2[u] = v_exit1[u]; + } + v_exit1[u] = en; + } + else if( + v_exit2.find(u) == v_exit2.end() || + comp( v_exit2[u].weight, en.weight ) + ) + { + if( + !is_ancestor( + v_id[target( en.edge, g )], + v_id[source( en.edge, g )] + ) + ) + { + v_exit2[u] = en; + } + } + } + else + { + v_exit1[u] = en; + } + } + } + } + f_heap_t tmp_heap; + BOOST_FOREACH( typename exit_map_t::value_type& t, b_exit ) + { + tmp_heap.push( t.second ); + } + BOOST_FOREACH( typename exit_map_t::value_type& t, v_exit1 ) + { + tmp_heap.push( t.second ); + } + BOOST_FOREACH( typename exit_map_t::value_type& t, v_exit2 ) + { + tmp_heap.push( t.second ); + } + in_edges[v].swap( tmp_heap ); + } + + BOOST_FOREACH( Vertex v, cycle_vertex_set ) + { + if( v != new_repr ) in_edges[new_repr].merge( in_edges[v] ); + } + + unvisited_vertex_set.insert( new_repr ); + + } + + } + + } + + if( delta ) + { + *p = std::make_pair( return_edge, *delta ); + return p; + } + else + { + p.reset(); + return p; + } + + } + + template + class branching_filter + { + + public: + branching_filter(){} + branching_filter( const EdgeSet * _es ) : p_es( _es ){} + + template + bool operator()( const Edge& e ) const + { + if( p_es->find( e ) != p_es->end() ) + return true; + else + return false; + } + + private: + const EdgeSet * p_es; + + }; + + template + struct BranchingEntry + { + Edge edge; + typename property_traits::value_type weight; + Compare compare; + unordered_set branching; + unordered_set include_edges; + unordered_set exclude_edges; + BranchingEntry( + const typename property_traits::value_type& w, + const Edge& e, + Compare& comp, + const unordered_set& b, + const unordered_set& include, + const unordered_set& exclude + ) : + edge( e ), weight( w ), compare( comp ), + branching( b ), include_edges( include ), + exclude_edges( exclude ){} + bool operator<( BranchingEntry const & rhs ) const + { return compare( weight, rhs.weight ); } + }; + + template + typename property_traits::value_type + compute_branching_weight( + WeightMap& w, + const unordered_set& branching + ) + { + + typedef typename property_traits::value_type weight_t; + + scoped_ptr weight; + weight.reset(); + + BOOST_FOREACH( const Edge& e, branching ) + { + + if( !weight ) + { + weight.reset( new weight_t ); + *weight = get( w, e ); + } + else + { + *weight += get( w, e ); + } + } + + return *weight; + + } + + template + void + rank_spanning_branchings_impl( const Graph& g, + BranchingProcessor bp, + IndexMap v_id, + WeightMap w, + Compare comp, + Rank rank, + Parent pred1, + Parent pred2 + ) + { + + typedef typename graph_traits::vertex_descriptor Vertex; + typedef typename graph_traits::edge_descriptor Edge; + + typedef typename property_traits::value_type weight_t; + + BOOST_CONCEPT_ASSERT(( VertexAndEdgeListGraphConcept )); + + BOOST_CONCEPT_ASSERT(( ReadablePropertyMapConcept )); + BOOST_CONCEPT_ASSERT(( ReadablePropertyMapConcept )); + + unordered_set best_branching, empty_set; + + shared_ptr< std::pair > p; + + Edge e; + + heap::fibonacci_heap > Q; + + best_spanning_branching( g, + best_branching, + v_id, + w, + comp, + rank, + pred1, + pred2, + empty_set, + empty_set + ); + + branching_filter > filter1( &best_branching ); + filtered_graph > > fg1( g, filter1 ); + + if( !bp( fg1 ) ) return; + + p = next_spanning_branching( g, + best_branching, + v_id, + w, + comp, + rank, + pred1, + pred2, + empty_set, + empty_set + ); + + if( p ) + { + Q.push( + BranchingEntry( + compute_branching_weight( w, best_branching ) - (*p.get()).second, + (*p.get()).first, + comp, + best_branching, + empty_set, + empty_set + ) + ); + } + else + { + return; + } + + while( !Q.empty() ) + { + + unordered_set branching; + + BranchingEntry P = Q.top(); + + Q.pop(); + + unordered_set include_edges = P.include_edges; + + unordered_set exclude_edges = P.exclude_edges; + + include_edges.insert( P.edge ); + + exclude_edges.insert( P.edge ); + + best_spanning_branching( g, + branching, + v_id, + w, + comp, + rank, + pred1, + pred2, + P.include_edges, + exclude_edges + ); + + branching_filter > filter( &branching ); + filtered_graph > > fg( g, filter ); + + if( !bp( fg ) ) return; + + p = + next_spanning_branching( g, + P.branching, + v_id, + w, + comp, + rank, + pred1, + pred2, + include_edges, + P.exclude_edges + ); + + if( p ) + { + Q.push( + BranchingEntry( + compute_branching_weight( w, P.branching ) - (*p.get()).second, + (*p.get()).first, + comp, + P.branching, + include_edges, + P.exclude_edges + ) + ); + } + + p = + next_spanning_branching( g, + branching, + v_id, + w, + comp, + rank, + pred1, + pred2, + P.include_edges, + exclude_edges + ); + + if( p ) + { + Q.push( + BranchingEntry( + P.weight - (*p.get()).second, + (*p.get()).first, + comp, + branching, + P.include_edges, + exclude_edges + ) + ); + } + + } + + } + + template + void + rank_spanning_branchings_dispatch2( const Graph& g, + BranchingProcessor bp, + IndexMap id_map, + WeightMap w_map, + Compare compare + ) + { + + typename graph_traits::vertices_size_type n = num_vertices(g); + + if( num_vertices( g ) == 0 ) return; // Nothing to do. + + typedef typename graph_traits::vertices_size_type vertex_idx_t; + typedef typename graph_traits::vertex_descriptor Vertex; + + // Set up rank and parent for disjoint sets. + + std::vector rank( n ); + + std::vector pred1( n ), pred2( n ); + + rank_spanning_branchings_impl( + g, + bp, + id_map, + w_map, + compare, + make_iterator_property_map( rank.begin(), id_map, rank[0] ), + make_iterator_property_map( pred1.begin(), id_map, pred1[0] ), + make_iterator_property_map( pred2.begin(), id_map, pred2[0]) + ); + + } + + template + void rank_spanning_branchings_dispatch1( + const Graph& g, + BranchingProcessor bp, + Compare compare, + const bgl_named_params& params + ) + { + + detail::rank_spanning_branchings_dispatch2( + g, + bp, + choose_param( + get_param( params, vertex_index_t()), get( vertex_index, g ) + ), + choose_param( + get_param( params, edge_weight_t()), get( edge_weight, g ) + ), + compare + ); + + } + + template + void rank_spanning_branchings_dispatch1( + const Graph& g, + BranchingProcessor bp, + param_not_found, + const bgl_named_params& params + ) + { + + typedef + typename + property_traits< + typename property_map::const_type + >::value_type weight_t; + + BOOST_CONCEPT_ASSERT(( ComparableConcept )); + + detail::rank_spanning_branchings_dispatch2( + g, + bp, + choose_param( + get_param( params, vertex_index_t()), get( vertex_index, g ) + ), + choose_param( + get_param( params, edge_weight_t()), get( edge_weight, g ) + ), + std::less() + ); + + } + + } // namespace detail + + template + inline void + rank_spanning_branchings( + const Graph& g, + BranchingProcessor bp, + const bgl_named_params& params + ) + { + + detail::rank_spanning_branchings_dispatch1( + g, + bp, + get_param( params, distance_compare_t() ), + params + ); + + } + + template + inline void + rank_spanning_branchings( const Graph& g, + Function bp + ) + { + + bgl_named_params params(0); + rank_spanning_branchings( g, bp, params ); + + } + +} // namespace boost + +#endif // BOOST_GRAPH_RANK_SPANNING_BRANCHINGS_HPP diff --git a/test/Jamfile.v2 b/test/Jamfile.v2 index 0523f100..9ff88699 100644 --- a/test/Jamfile.v2 +++ b/test/Jamfile.v2 @@ -133,6 +133,8 @@ test-suite graph_test : [ run hawick_circuits.cpp ] [ run successive_shortest_path_nonnegative_weights_test.cpp ../../test/build//boost_unit_test_framework/static ] [ run cycle_canceling_test.cpp ../../test/build//boost_unit_test_framework/static ] + [ run strong_components_test.cpp ] + [ run rank_branchings_test.cpp : 4 10 ] ; # Run SDB tests only when -sSDB= is set. diff --git a/test/rank_branchings_test.cpp b/test/rank_branchings_test.cpp new file mode 100644 index 00000000..36edbff8 --- /dev/null +++ b/test/rank_branchings_test.cpp @@ -0,0 +1,835 @@ +//======================================================================= +// Copyright 2015-2016 Clemson University +// Authors: Bradley S. Meyer +// +// Distributed under the Boost Software License, Version 1.0. (See +// accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) +//======================================================================= + +// This test first creates a complete graph or multigraph with n vertices. +// Each arc is assigned a random integer, double, complex, or defined class +// weight. A root vertex is added and arcs with zero weight +// are added from the root vertex to each of the other vertices. The test +// then constructs all potential branchings and discards those +// that have indegree for any vertex greater than one or contain +// a cycle. The valid branchings are then stored in a priority +// queue. The test then +// uses the rank_spanning_branchings.hpp routine to construct all +// branchings, in order, and places them in a vector. The vector +// and the priority queue are then compared, allowing for the possibility +// of equal weight branchings. + +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +using namespace boost; + +typedef random::mt19937 base_generator_type; + +template +class my_type +{ + public: + + my_type(){} + my_type( T _t, U _u ) : t( _t ), u( _u ) {} + + my_type& operator+=( const my_type& rhs ) + { + this->t += rhs.t; + this->u += rhs.u; + return *this; + } + + friend my_type operator+( my_type lhs, const my_type& rhs ) + { + lhs += rhs; + return lhs; + } + + my_type& operator-=( const my_type& rhs ) + { + this->t -= rhs.t; + this->u -= rhs.u; + return *this; + } + + friend my_type operator-( my_type lhs, const my_type& rhs ) + { + lhs -= rhs; + return lhs; + } + + bool operator!=( const my_type& rhs ) + { + return !( ( this->t == rhs.t ) && ( this->u == rhs.u ) ); + } + + bool operator<( const my_type& rhs ) const + { + if( this->t == rhs.t ) + return this->u < rhs.u; + else + return this->t < rhs.t; + } + + friend std::ostream& operator<<( std::ostream& os, const my_type& p ) + { + os << "(" << p.t << ", " << p.u << ")"; + return os; + } + + private: + T t; + U u; + +}; + +template +struct my_compare +{ + bool operator()( const std::complex& a, const std::complex& b) const + { + if( a.real() == b.real() ) + return a.imag() < b.imag( ); + else + return a.real() < b.real( ); + } +}; + +template +struct Branching +{ + typedef + typename + property_traits< + typename property_map::const_type + >::value_type weight_t; + weight_t weight; + std::set v_edges; + Branching( weight_t & w, const std::set& v_e) : + weight( w ), v_edges( v_e ){} +}; + +template +struct compare_branchings +{ + Compare comp; + compare_branchings( Compare _comp ) : comp( _comp ){} + bool operator()( + const Branching& b1, + const Branching& b2 ) const + { + return comp( b1.weight, b2.weight ); + } +}; + + +template +struct set_rank_vector +{ + + std::vector >& rank_vector; + set_rank_vector( + std::vector >& rv + ) : rank_vector( rv ) {} + + template + bool operator()( BranchingGraph& bg ) + { + + std::set branching; + + typedef + typename + property_map::const_type weight_map_t; + weight_map_t w; + typedef typename property_traits::value_type weight_t; + + scoped_ptr weight; + weight.reset(); + + BGL_FORALL_EDGES_T( e, bg, BranchingGraph ) + { + if( !weight ) + { + weight.reset( new weight_t ); + *weight = get( w, e ); + } + else + { + *weight += get( w, e ); + } + branching.insert( e ); + } + + rank_vector.push_back(Branching( *weight, branching ) ); + + return true; + + } + +}; + +template +void +set_branching( + Graph& g, + const std::vector& v_edges, + const std::vector& combination, + Rank rank, + Parent parent, + std::vector >& branchings +) +{ + disjoint_sets W( rank, parent ); + typedef typename graph_traits::vertex_descriptor Vertex; + std::set in_vertex_set; + typedef typename property_map::const_type weight_map_t; + typedef typename property_traits::value_type weight_t; + std::set branching_edges; + + BGL_FORALL_VERTICES_T( v, g, Graph ) + { + W.make_set( v ); + } + + BOOST_FOREACH( size_t i, combination ) + { + Edge e = v_edges[i]; + if( in_vertex_set.find( target( e, g ) ) != in_vertex_set.end() ) return; + if( W.find_set( source( e, g ) ) == W.find_set( target( e, g ) ) ) return; + W.union_set( source( e, g ), target( e, g ) ); + in_vertex_set.insert( target( e, g ) ); + } + + scoped_ptr weight; + weight.reset(); + + BOOST_FOREACH( size_t i, combination ) + { + Edge e = v_edges[i]; + if( !weight ) + { + weight.reset( new weight_t ); + *weight = get( edge_weight, g, e ); + } + else + { + *weight += get( edge_weight, g, e ); + } + branching_edges.insert( e ); + } + branchings.push_back( Branching( *weight, branching_edges ) ); +} + +template +void +loop( + size_t offset, + size_t k, + Graph& g, + const std::vector& v_edges, + const std::vector& indices, + std::vector& combination, + std::vector >& branchings +) +{ + + typedef typename graph_traits::vertices_size_type size_type; + typedef typename graph_traits::vertex_descriptor vertex_t; + + size_type n = num_vertices( g ); + std::vector rank_map( n ); + std::vector pred_map( n ); + + if( k == 0 ) + { + set_branching( + g, + v_edges, + combination, + make_iterator_property_map( + rank_map.begin(), + get(vertex_index, g), + rank_map[0] + ), + make_iterator_property_map( + pred_map.begin(), + get(vertex_index, g), + pred_map[0] + ), + branchings + ); + return; + } + for (size_t i = offset; i <= indices.size() - k; ++i) { + combination.push_back(indices[i]); + loop( + i+1, + k-1, + g, + v_edges, + indices, + combination, + branchings + ); + combination.pop_back(); + } +} + +template +bool +compare_heap_and_vector( + Graph& g, + heap::fibonacci_heap< + Branching, + heap::compare + >& branching_heap, + std::vector >& rank_vector +) +{ + + // Compare branchings in heap and ranked vector. Allow for the possibility + // that two branchings might have the same weight (hence the second loop). + + typename property_map::type index_map = + get( vertex_index, g ); + + bool found; + + std::set check; + + for( size_t i = 0; i < rank_vector.size(); i++ ) + { + check.insert( i ); + } + + while( !branching_heap.empty() ) + { + + for( + std::set::iterator it = check.begin(); + it != check.end(); + it++ + ) + { + + found = true; + + if( rank_vector[*it].weight != branching_heap.top().weight ) + { + found = false; + break; + } + + BOOST_FOREACH( Edge e, branching_heap.top().v_edges ) + { + if( + rank_vector[*it].v_edges.find( e ) == rank_vector[*it].v_edges.end() + ) + { + found = false; + } + } + + if( found ) + { + check.erase( it ); + break; + } + + } + + if( !found ) + { + std::cerr << "Error: " << std::endl; + BOOST_FOREACH( Edge e, branching_heap.top().v_edges ) + { + std::cerr << + " (" << index_map[source( e, g )] << + "," << index_map[target( e, g )] << ") "; + } + std::cerr << " Weight: " << branching_heap.top().weight << std::endl; + std::cerr << std::endl; + return false; + } + branching_heap.pop(); + } + + return true; + +} + +template +void test_int_graph( base_generator_type& gen, size_t n, Compare comp ) +{ + + typedef adjacency_list < listS, setS, directedS, + property < vertex_index_t, size_t >, property < edge_weight_t, int > > + Graph; + typedef property_map < Graph, vertex_index_t >::type index_map_t; + typedef graph_traits ::edge_descriptor Edge; + + std::vector indices; + std::vector combination; + std::vector v_edges; + std::vector > branchings, rank_vector; + + typedef + heap::fibonacci_heap< + Branching, + heap::compare > + > branching_heap_t; + + branching_heap_t branching_heap( comp ); + + random::uniform_int_distribution<> tdist( -10, 10 ); + + variate_generator< + base_generator_type&, random::uniform_int_distribution<> + > rvt( gen, tdist ); + + Graph g; + + for( size_t i = 0; i <= n; i++ ) + { + add_vertex( g ); + } + + // Set index map. + + index_map_t index_map = get( vertex_index, g ); + + size_t j = 0; + BGL_FORALL_VERTICES( v, g, Graph ) + { + put( index_map, v, j++ ); + } + + // Add edges. + + BGL_FORALL_VERTICES( v, g, Graph ) + { + BGL_FORALL_VERTICES( u, g, Graph ) + { + if( u != v ) + { + if( index_map[v] != 0 ) + { + if( index_map[u] != 0 ) + { + add_edge( v, u, rvt(), g ); + } + } + else + { + add_edge( v, u, 0, g ); + } + } + } + } + + // Create vectors. + + BGL_FORALL_EDGES( e, g, Graph ) + { + v_edges.push_back( e ); + } + + for( size_t i = 0; i < v_edges.size(); ++i ) + { + indices.push_back(i); + } + + loop( + 0, + num_vertices(g) - 1, + g, + v_edges, + indices, + combination, + branchings + ); + + for( size_t i = 0; i < branchings.size(); i++ ) + { + branching_heap.push( branchings[i] ); + } + + rank_spanning_branchings( + g, + set_rank_vector( rank_vector ), + distance_compare( comp ) + ); + + // Check that number of branchings found is correct. + + assert( + pow( + static_cast( n + 1 ), + static_cast( n - 1 ) + ) + == + rank_vector.size() + ); + + assert( branching_heap.size() == rank_vector.size() ); + + // Compare branchings found by both methods. + + if( !compare_heap_and_vector( g, branching_heap, rank_vector ) ) + { + exit( EXIT_FAILURE ); + } + +} + +template +void test_real_graph( base_generator_type& gen, size_t n, Compare comp ) +{ + + typedef adjacency_list < vecS, vecS, directedS, + no_property, property < edge_weight_t, double > > Graph; + typedef graph_traits ::edge_descriptor Edge; + + std::vector indices; + std::vector combination; + std::vector v_edges; + std::vector > branchings, rank_vector; + + typedef + heap::fibonacci_heap< + Branching, + heap::compare > + > branching_heap_t; + + branching_heap_t branching_heap( comp ); + + boost::random::uniform_real_distribution<> tdist( -10, 10 ); + + boost::variate_generator< + base_generator_type&, boost::random::uniform_real_distribution<> + > rvt( gen, tdist ); + + Graph g; + + for( size_t i = 0; i < n; i++ ) + { + for( size_t j = 0; j < n; j++ ) + { + if( i != j ) + { + add_edge( i, j, rvt(), g ); + add_edge( i, j, rvt(), g ); // A parallel edge. + } + } + add_edge( n, i, 0, g ); + } + + BGL_FORALL_EDGES( e, g, Graph ) + { + v_edges.push_back( e ); + } + + for( size_t i = 0; i < v_edges.size(); ++i ) + { + indices.push_back(i); + } + + loop( + 0, + num_vertices(g) - 1, + g, + v_edges, + indices, + combination, + branchings + ); + + for( size_t i = 0; i < branchings.size(); i++ ) + { + branching_heap.push( branchings[i] ); + } + + // Call explicitly with named parameters to test. + + rank_spanning_branchings( + g, + set_rank_vector( rank_vector ), + distance_compare( comp ). + weight_map( get( edge_weight, g ) ). + vertex_index_map( get( vertex_index, g ) ) + ); + + // Check that number of branchings found is correct. The formula + // accounts for the parallel edges. + + assert( + pow( + static_cast( 2*n + 1 ), + static_cast( n - 1 ) + ) + == + rank_vector.size() + ); + + assert( branchings.size() == rank_vector.size() ); + + // Compare branchings found by both methods. + + if( !compare_heap_and_vector( g, branching_heap, rank_vector ) ) + { + exit( EXIT_FAILURE ); + } + +} + +template +void test_complex_graph( base_generator_type& gen, size_t n, Compare comp ) +{ + + typedef adjacency_list < vecS, vecS, directedS, + no_property, property < edge_weight_t, std::complex > > Graph; + typedef graph_traits ::edge_descriptor Edge; + + std::vector indices; + std::vector combination; + std::vector v_edges; + std::vector > branchings, rank_vector; + + typedef + heap::fibonacci_heap< + Branching, + heap::compare > + > branching_heap_t; + + branching_heap_t branching_heap( comp ); + + boost::random::uniform_real_distribution<> tdist( -10, 10 ); + + boost::variate_generator< + base_generator_type&, boost::random::uniform_real_distribution<> + > rvt( gen, tdist ); + + Graph g; + + for( size_t i = 0; i < n; i++ ) + { + for( size_t j = 0; j < n; j++ ) + { + if( i != j ) + { + add_edge( i, j, std::complex( rvt(), rvt() ), g ); + } + } + add_edge( n, i, std::complex( 0, 0 ), g ); + } + + BGL_FORALL_EDGES( e, g, Graph ) + { + v_edges.push_back( e ); + } + + for( size_t i = 0; i < v_edges.size(); ++i ) + { + indices.push_back(i); + } + + loop( + 0, + num_vertices(g) - 1, + g, + v_edges, + indices, + combination, + branchings + ); + + for( size_t i = 0; i < branchings.size(); i++ ) + { + branching_heap.push( branchings[i] ); + } + + rank_spanning_branchings( + g, + set_rank_vector( rank_vector ), + distance_compare( comp ) + ); + + // Check that number of branchings found is correct. + + assert( + pow( + static_cast( n + 1 ), + static_cast( n - 1 ) + ) + == + rank_vector.size() + ); + + assert( branchings.size() == rank_vector.size() ); + + // Compare branchings found by both methods. + + if( !compare_heap_and_vector( g, branching_heap, rank_vector ) ) + { + exit( EXIT_FAILURE ); + } + +} + +template +void test_class_graph( base_generator_type& gen, size_t n, Compare comp ) +{ + + typedef adjacency_list < vecS, vecS, directedS, + no_property, property < edge_weight_t, my_type > + > Graph; + typedef graph_traits ::edge_descriptor Edge; + + std::vector indices; + std::vector combination; + std::vector v_edges; + std::vector > branchings, rank_vector; + + typedef + heap::fibonacci_heap< + Branching, + heap::compare > + > branching_heap_t; + + branching_heap_t branching_heap( comp ); + + boost::random::uniform_real_distribution<> tdist( -10, 10 ); + + boost::variate_generator< + base_generator_type&, boost::random::uniform_real_distribution<> + > rvt( gen, tdist ); + + Graph g; + + for( size_t i = 0; i < n; i++ ) + { + for( size_t j = 0; j < n; j++ ) + { + if( i != j ) + { + add_edge( i, j, my_type( rvt(), rvt() ), g ); + } + } + add_edge( n, i, my_type( 0, 0 ), g ); + } + + BGL_FORALL_EDGES( e, g, Graph ) + { + v_edges.push_back( e ); + } + + for( size_t i = 0; i < v_edges.size(); ++i ) + { + indices.push_back(i); + } + + loop( + 0, + num_vertices(g) - 1, + g, + v_edges, + indices, + combination, + branchings + ); + + for( size_t i = 0; i < branchings.size(); i++ ) + { + branching_heap.push( branchings[i] ); + } + + // Call explicitly with named parameters to test. + + rank_spanning_branchings( + g, + set_rank_vector( rank_vector ), + distance_compare( comp ). + weight_map( get( edge_weight, g ) ). + vertex_index_map( get( vertex_index, g ) ) + ); + + // Check that number of branchings found is correct. + + assert( + pow( + static_cast( n + 1 ), + static_cast( n - 1 ) + ) + == + rank_vector.size() + ); + + assert( branchings.size() == rank_vector.size() ); + + // Compare branchings found by both methods. + + if( !compare_heap_and_vector( g, branching_heap, rank_vector ) ) + { + exit( EXIT_FAILURE ); + } + +} + +int main( int argc, char **argv ) +{ + + if( argc != 3 ) + { + std::cerr << std::endl; + std::cerr << " Usage: " << argv[0] << " n" << std::endl << std::endl; + std::cerr << " n = number of vertices in complete graph" << + std::endl << std::endl; + std::cerr << " k = number of trials" << std::endl << std::endl; + return EXIT_FAILURE; + } + + base_generator_type gen; + gen.seed( (unsigned int) time( NULL ) ); + + for( size_t i = 0; i < lexical_cast( argv[2] ); i++ ) + { + test_int_graph( + gen, lexical_cast( argv[1] ) + 1, std::greater() + ); + test_real_graph( + gen, lexical_cast( argv[1] ), std::less() + ); + test_complex_graph( + gen, lexical_cast( argv[1] ) + 1, my_compare() + ); + test_class_graph( + gen, + lexical_cast( argv[1] ) + 1, + std::less >() + ); + } + + std::cout << "tests passed" << std::endl; + + return EXIT_SUCCESS; + +} +