+P. M. Camerini, L. Fratta, and F. Maffioli
+A Note on Finding Optimum Branchings Boost Graph Library: 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 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 example/rank-branching2.cpp
Acknowledgments
+
+
+The development of rank_spanning_branchings() was inspired by
+the C++ implementation of
+Edmond's algorithm by Ali Tofigh and Erik Sjölund.
+
+
+
+#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/
+#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 
+
+