source: NonGTP/Boost/boost/graph/minimum_degree_ordering.hpp @ 857

Revision 857, 23.1 KB checked in by igarcia, 18 years ago (diff)
Line 
1//-*-c++-*-
2//=======================================================================
3// Copyright 1997-2001 University of Notre Dame.
4// Authors: Lie-Quan Lee, Jeremy Siek
5//
6// Distributed under the Boost Software License, Version 1.0. (See
7// accompanying file LICENSE_1_0.txt or copy at
8// http://www.boost.org/LICENSE_1_0.txt)
9//=======================================================================
10//
11#ifndef MINIMUM_DEGREE_ORDERING_HPP
12#define MINIMUM_DEGREE_ORDERING_HPP
13
14#include <vector>
15#include <cassert>
16#include <boost/config.hpp>
17#include <boost/pending/bucket_sorter.hpp>
18#include <boost/detail/numeric_traits.hpp> // for integer_traits
19#include <boost/graph/graph_traits.hpp>
20#include <boost/property_map.hpp>
21
22namespace boost {
23
24  namespace detail {
25
26    //
27    // Given a set of n integers (where the integer values range from
28    // zero to n-1), we want to keep track of a collection of stacks
29    // of integers. It so happens that an integer will appear in at
30    // most one stack at a time, so the stacks form disjoint sets.
31    // Because of these restrictions, we can use one big array to
32    // store all the stacks, intertwined with one another.
33    // No allocation/deallocation happens in the push()/pop() methods
34    // so this is faster than using std::stack's.
35    //
36    template <class SignedInteger>
37    class Stacks {
38      typedef SignedInteger value_type;
39      typedef typename std::vector<value_type>::size_type size_type;
40    public:
41      Stacks(size_type n) : data(n) {}
42     
43      //: stack
44      class stack {
45        typedef typename std::vector<value_type>::iterator Iterator;
46      public:
47        stack(Iterator _data, const value_type& head)
48          :  data(_data), current(head) {}
49
50        // did not use default argument here to avoid internal compiler error
51        // in g++.
52        stack(Iterator _data)
53          : data(_data), current(-(std::numeric_limits<value_type>::max)()) {}
54       
55        void pop() {
56          assert(! empty());
57          current = data[current];
58        }
59        void push(value_type v) {
60          data[v] = current;
61          current = v;
62        }
63        bool empty() {
64          return current == -(std::numeric_limits<value_type>::max)();
65        }
66        value_type& top() { return current; }
67      private:
68        Iterator data;
69        value_type current;
70      };
71
72      // To return a stack object
73      stack make_stack()
74        { return stack(data.begin()); }
75    protected:
76      std::vector<value_type> data;
77    };
78
79
80    // marker class, a generalization of coloring.
81    //
82    // This class is to provide a generalization of coloring which has
83    // complexity of amortized constant time to set all vertices' color
84    // back to be untagged. It implemented by increasing a tag.
85    //
86    // The colors are:
87    //   not tagged
88    //   tagged
89    //   multiple_tagged
90    //   done
91    //
92    template <class SignedInteger, class Vertex, class VertexIndexMap>
93    class Marker {
94      typedef SignedInteger value_type;
95      typedef typename std::vector<value_type>::size_type size_type;
96     
97      static value_type done()
98      { return (std::numeric_limits<value_type>::max)()/2; }
99    public:
100      Marker(size_type _num, VertexIndexMap index_map)
101        : tag(1 - (std::numeric_limits<value_type>::max)()),
102          data(_num, - (std::numeric_limits<value_type>::max)()),
103          id(index_map) {}
104     
105      void mark_done(Vertex node) { data[get(id, node)] = done(); }
106     
107      bool is_done(Vertex node) { return data[get(id, node)] == done(); }
108     
109      void mark_tagged(Vertex node) { data[get(id, node)] = tag; }
110     
111      void mark_multiple_tagged(Vertex node) { data[get(id, node)] = multiple_tag; }
112 
113      bool is_tagged(Vertex node) const { return data[get(id, node)] >= tag; }
114
115      bool is_not_tagged(Vertex node) const { return data[get(id, node)] < tag; }
116
117      bool is_multiple_tagged(Vertex node) const
118        { return data[get(id, node)] >= multiple_tag; }
119
120      void increment_tag() {
121        const size_type num = data.size();
122        ++tag;
123        if ( tag >= done() ) {
124          tag = 1 - (std::numeric_limits<value_type>::max)();
125          for (size_type i = 0; i < num; ++i)
126            if ( data[i] < done() )
127              data[i] = - (std::numeric_limits<value_type>::max)();
128        }
129      }
130     
131      void set_multiple_tag(value_type mdeg0)
132      {
133        const size_type num = data.size();
134        multiple_tag = tag + mdeg0;
135       
136        if ( multiple_tag >= done() ) {
137          tag = 1-(std::numeric_limits<value_type>::max)();
138         
139          for (size_type i=0; i<num; i++)
140            if ( data[i] < done() )
141              data[i] = -(std::numeric_limits<value_type>::max)();
142         
143          multiple_tag = tag + mdeg0;
144        }
145      }
146     
147      void set_tag_as_multiple_tag() { tag = multiple_tag; }
148     
149    protected:
150      value_type tag;
151      value_type multiple_tag;
152      std::vector<value_type> data;
153      VertexIndexMap id;
154    };
155   
156    template< class Iterator, class SignedInteger,
157       class Vertex, class VertexIndexMap, int offset = 1 >
158    class Numbering {
159      typedef SignedInteger number_type;
160      number_type num; //start from 1 instead of zero
161      Iterator   data;
162      number_type max_num;
163      VertexIndexMap id;
164    public:
165      Numbering(Iterator _data, number_type _max_num, VertexIndexMap id)
166        : num(1), data(_data), max_num(_max_num), id(id) {}
167      void operator()(Vertex node) { data[get(id, node)] = -num; }
168      bool all_done(number_type i = 0) const { return num + i > max_num; }
169      void increment(number_type i = 1) { num += i; }
170      bool is_numbered(Vertex node) const {
171        return data[get(id, node)] < 0;
172      }
173      void indistinguishable(Vertex i, Vertex j) {
174        data[get(id, i)] = - (get(id, j) + offset);
175      }
176    };
177
178    template <class SignedInteger, class Vertex, class VertexIndexMap>
179    class degreelists_marker {
180    public:
181      typedef SignedInteger value_type;
182      typedef typename std::vector<value_type>::size_type size_type;
183      degreelists_marker(size_type n, VertexIndexMap id)
184        : marks(n, 0), id(id) {}
185      void mark_need_update(Vertex i) { marks[get(id, i)] = 1;  }
186      bool need_update(Vertex i) { return marks[get(id, i)] == 1; }
187      bool outmatched_or_done (Vertex i) { return marks[get(id, i)] == -1; }
188      void mark(Vertex i) { marks[get(id, i)] = -1; }
189      void unmark(Vertex i) { marks[get(id, i)] = 0; }
190    private:
191      std::vector<value_type> marks;
192      VertexIndexMap id;
193    };
194
195    // Helper function object for edge removal
196    template <class Graph, class MarkerP, class NumberD, class Stack,
197      class VertexIndexMap>
198    class predicateRemoveEdge1 {
199      typedef typename graph_traits<Graph>::vertex_descriptor vertex_t;
200      typedef typename graph_traits<Graph>::edge_descriptor edge_t;
201    public:
202      predicateRemoveEdge1(Graph& _g, MarkerP& _marker,
203                           NumberD _numbering, Stack& n_e, VertexIndexMap id)
204        : g(&_g), marker(&_marker), numbering(_numbering),
205          neighbor_elements(&n_e), id(id) {}
206
207      bool operator()(edge_t e) {
208        vertex_t dist = target(e, *g);
209        if ( marker->is_tagged(dist) )
210          return true;
211        marker->mark_tagged(dist);
212        if (numbering.is_numbered(dist)) {
213          neighbor_elements->push(get(id, dist));
214          return true;
215        }
216        return false;
217      }
218    private:
219      Graph*   g;
220      MarkerP* marker;
221      NumberD  numbering;
222      Stack*   neighbor_elements;
223      VertexIndexMap id;
224    };
225
226    // Helper function object for edge removal
227    template <class Graph, class MarkerP>
228    class predicate_remove_tagged_edges
229    {
230      typedef typename graph_traits<Graph>::vertex_descriptor vertex_t;
231      typedef typename graph_traits<Graph>::edge_descriptor edge_t;
232    public:
233      predicate_remove_tagged_edges(Graph& _g, MarkerP& _marker)
234        : g(&_g), marker(&_marker) {}
235
236      bool operator()(edge_t e) {
237        vertex_t dist = target(e, *g);
238        if ( marker->is_tagged(dist) )
239          return true;
240        return false;
241      }
242    private:
243      Graph*   g;
244      MarkerP* marker;
245    };
246
247    template<class Graph, class DegreeMap,
248             class InversePermutationMap,
249             class PermutationMap,
250             class SuperNodeMap,
251             class VertexIndexMap>
252    class mmd_impl
253    {
254      // Typedefs
255      typedef graph_traits<Graph> Traits;
256      typedef typename Traits::vertices_size_type size_type;
257      typedef typename detail::integer_traits<size_type>::difference_type
258        diff_t;
259      typedef typename Traits::vertex_descriptor vertex_t;
260      typedef typename Traits::adjacency_iterator adj_iter;
261      typedef iterator_property_map<vertex_t*,
262        identity_property_map, vertex_t, vertex_t&> IndexVertexMap;
263      typedef detail::Stacks<diff_t> Workspace;
264      typedef bucket_sorter<size_type, vertex_t, DegreeMap, VertexIndexMap>
265        DegreeLists;
266      typedef Numbering<InversePermutationMap, diff_t, vertex_t,VertexIndexMap>
267        NumberingD;
268      typedef degreelists_marker<diff_t, vertex_t, VertexIndexMap>
269        DegreeListsMarker;
270      typedef Marker<diff_t, vertex_t, VertexIndexMap> MarkerP;
271
272      // Data Members
273
274      // input parameters
275      Graph& G;
276      int delta;
277      DegreeMap degree;
278      InversePermutationMap inverse_perm;
279      PermutationMap perm;
280      SuperNodeMap supernode_size;
281      VertexIndexMap vertex_index_map;
282
283      // internal data-structures
284      std::vector<vertex_t> index_vertex_vec;
285      size_type n;
286      IndexVertexMap index_vertex_map;
287      DegreeLists degreelists;
288      NumberingD numbering;
289      DegreeListsMarker degree_lists_marker;
290      MarkerP marker;
291      Workspace work_space;
292    public:
293      mmd_impl(Graph& g, size_type n_, int delta, DegreeMap degree,
294               InversePermutationMap inverse_perm,
295               PermutationMap perm,
296               SuperNodeMap supernode_size,
297               VertexIndexMap id)
298        : G(g), delta(delta), degree(degree),
299        inverse_perm(inverse_perm),
300        perm(perm),
301        supernode_size(supernode_size),
302        vertex_index_map(id),
303        index_vertex_vec(n_),
304        n(n_),
305        degreelists(n_ + 1, n_, degree, id),
306        numbering(inverse_perm, n_, vertex_index_map),
307        degree_lists_marker(n_, vertex_index_map),
308        marker(n_, vertex_index_map),
309        work_space(n_)
310      {
311        typename graph_traits<Graph>::vertex_iterator v, vend;
312        size_type vid = 0;
313        for (tie(v, vend) = vertices(G); v != vend; ++v, ++vid)
314          index_vertex_vec[vid] = *v;
315        index_vertex_map = IndexVertexMap(&index_vertex_vec[0]);
316
317        // Initialize degreelists.  Degreelists organizes the nodes
318        // according to their degree.
319        for (tie(v, vend) = vertices(G); v != vend; ++v) {
320          put(degree, *v, out_degree(*v, G));
321          degreelists.push(*v);
322        }
323      }
324
325      void do_mmd()
326      {
327        // Eliminate the isolated nodes -- these are simply the nodes
328        // with no neighbors, which are accessible as a list (really, a
329        // stack) at location 0.  Since these don't affect any other
330        // nodes, we can eliminate them without doing degree updates.
331        typename DegreeLists::stack list_isolated = degreelists[0];
332        while (!list_isolated.empty()) {
333          vertex_t node = list_isolated.top();
334          marker.mark_done(node);
335          numbering(node);
336          numbering.increment();
337          list_isolated.pop();
338        }
339        size_type min_degree = 1;
340        typename DegreeLists::stack list_min_degree = degreelists[min_degree];
341
342        while (list_min_degree.empty()) {
343          ++min_degree;
344          list_min_degree = degreelists[min_degree];
345        }
346
347        // check if the whole eliminating process is done
348        while (!numbering.all_done()) {
349
350          size_type min_degree_limit = min_degree + delta; // WARNING
351          typename Workspace::stack llist = work_space.make_stack();
352
353          // multiple elimination
354          while (delta >= 0) {
355
356            // Find the next non-empty degree
357            for (list_min_degree = degreelists[min_degree];
358                 list_min_degree.empty() && min_degree <= min_degree_limit;
359                 ++min_degree, list_min_degree = degreelists[min_degree])
360              ;
361            if (min_degree > min_degree_limit)
362              break;
363
364            const vertex_t node = list_min_degree.top();
365            const size_type node_id = get(vertex_index_map, node);
366            list_min_degree.pop();
367            numbering(node);
368
369            // check if node is the last one
370            if (numbering.all_done(supernode_size[node])) {
371              numbering.increment(supernode_size[node]);
372              break;
373            }
374            marker.increment_tag();
375            marker.mark_tagged(node);
376
377            this->eliminate(node);
378
379            numbering.increment(supernode_size[node]);
380            llist.push(node_id);
381          } // multiple elimination
382
383          if (numbering.all_done())
384            break;
385
386          this->update( llist, min_degree);
387        }
388
389      } // do_mmd()
390
391      void eliminate(vertex_t node)
392      {
393        typename Workspace::stack element_neighbor = work_space.make_stack();
394
395        // Create two function objects for edge removal
396        typedef typename Workspace::stack WorkStack;
397        predicateRemoveEdge1<Graph, MarkerP, NumberingD,
398                             WorkStack, VertexIndexMap>
399          p(G, marker, numbering, element_neighbor, vertex_index_map);
400
401        predicate_remove_tagged_edges<Graph, MarkerP> p2(G, marker);
402
403        // Reconstruct the adjacent node list, push element neighbor in a List.
404        remove_out_edge_if(node, p, G);
405        //during removal element neighbors are collected.
406
407        while (!element_neighbor.empty()) {
408          // element absorb
409          size_type e_id = element_neighbor.top();
410          vertex_t element = get(index_vertex_map, e_id);
411          adj_iter i, i_end;
412          for (tie(i, i_end) = adjacent_vertices(element, G); i != i_end; ++i){
413            vertex_t i_node = *i;
414            if (!marker.is_tagged(i_node) && !numbering.is_numbered(i_node)) {
415              marker.mark_tagged(i_node);
416              add_edge(node, i_node, G);
417            }
418          }
419          element_neighbor.pop();
420        }
421        adj_iter v, ve;
422        for (tie(v, ve) = adjacent_vertices(node, G); v != ve; ++v) {
423          vertex_t v_node = *v;
424          if (!degree_lists_marker.need_update(v_node)
425              && !degree_lists_marker.outmatched_or_done(v_node)) {
426            degreelists.remove(v_node);
427          }
428          //update out edges of v_node
429          remove_out_edge_if(v_node, p2, G);
430
431          if ( out_degree(v_node, G) == 0 ) { // indistinguishable nodes
432            supernode_size[node] += supernode_size[v_node];
433            supernode_size[v_node] = 0;
434            numbering.indistinguishable(v_node, node);
435            marker.mark_done(v_node);
436            degree_lists_marker.mark(v_node);
437          } else {                           // not indistinguishable nodes
438            add_edge(v_node, node, G);
439            degree_lists_marker.mark_need_update(v_node);
440          }
441        }
442      } // eliminate()
443
444
445      template <class Stack>
446      void update(Stack llist, size_type& min_degree)
447      {
448        size_type min_degree0 = min_degree + delta + 1;
449
450        while (! llist.empty()) {
451          size_type deg, deg0 = 0;
452
453          marker.set_multiple_tag(min_degree0);
454          typename Workspace::stack q2list = work_space.make_stack();
455          typename Workspace::stack qxlist = work_space.make_stack();
456
457          vertex_t current = get(index_vertex_map, llist.top());
458          adj_iter i, ie;
459          for (tie(i,ie) = adjacent_vertices(current, G); i != ie; ++i) {
460            vertex_t i_node = *i;
461            const size_type i_id = get(vertex_index_map, i_node);
462            if (supernode_size[i_node] != 0) {
463              deg0 += supernode_size[i_node];
464              marker.mark_multiple_tagged(i_node);
465              if (degree_lists_marker.need_update(i_node)) {
466                if (out_degree(i_node, G) == 2)
467                  q2list.push(i_id);
468                else
469                  qxlist.push(i_id);
470              }
471            }
472          }
473
474          while (!q2list.empty()) {
475            const size_type u_id = q2list.top();
476            vertex_t u_node = get(index_vertex_map, u_id);
477            // if u_id is outmatched by others, no need to update degree
478            if (degree_lists_marker.outmatched_or_done(u_node)) {
479              q2list.pop();
480              continue;
481            }
482            marker.increment_tag();
483            deg = deg0;
484
485            adj_iter nu = adjacent_vertices(u_node, G).first;
486            vertex_t neighbor = *nu;
487            if (neighbor == u_node) {
488              ++nu;
489              neighbor = *nu;
490            }
491            if (numbering.is_numbered(neighbor)) {
492              adj_iter i, ie;
493              for (tie(i,ie) = adjacent_vertices(neighbor, G);
494                   i != ie; ++i) {
495                const vertex_t i_node = *i;
496                if (i_node == u_node || supernode_size[i_node] == 0)
497                  continue;
498                if (marker.is_tagged(i_node)) {
499                  if (degree_lists_marker.need_update(i_node)) {
500                    if ( out_degree(i_node, G) == 2 ) { // is indistinguishable
501                      supernode_size[u_node] += supernode_size[i_node];
502                      supernode_size[i_node] = 0;
503                      numbering.indistinguishable(i_node, u_node);
504                      marker.mark_done(i_node);
505                      degree_lists_marker.mark(i_node);
506                    } else                        // is outmatched
507                      degree_lists_marker.mark(i_node);
508                  }
509                } else {
510                  marker.mark_tagged(i_node);
511                  deg += supernode_size[i_node];
512                }
513              }
514            } else
515              deg += supernode_size[neighbor];
516
517            deg -= supernode_size[u_node];
518            degree[u_node] = deg; //update degree
519            degreelists[deg].push(u_node);
520            //u_id has been pushed back into degreelists
521            degree_lists_marker.unmark(u_node);
522            if (min_degree > deg)
523              min_degree = deg;
524            q2list.pop();
525          } // while (!q2list.empty())
526
527          while (!qxlist.empty()) {
528            const size_type u_id = qxlist.top();
529            const vertex_t u_node = get(index_vertex_map, u_id);
530
531            // if u_id is outmatched by others, no need to update degree
532            if (degree_lists_marker.outmatched_or_done(u_node)) {
533              qxlist.pop();
534              continue;
535            }
536            marker.increment_tag();
537            deg = deg0;
538            adj_iter i, ie;
539            for (tie(i, ie) = adjacent_vertices(u_node, G); i != ie; ++i) {
540              vertex_t i_node = *i;
541              if (marker.is_tagged(i_node))
542                continue;
543              marker.mark_tagged(i_node);
544
545              if (numbering.is_numbered(i_node)) {
546                adj_iter j, je;
547                for (tie(j, je) = adjacent_vertices(i_node, G); j != je; ++j) {
548                  const vertex_t j_node = *j;
549                  if (marker.is_not_tagged(j_node)) {
550                    marker.mark_tagged(j_node);
551                    deg += supernode_size[j_node];
552                  }
553                }
554              } else
555                deg += supernode_size[i_node];
556            } // for adjacent vertices of u_node
557            deg -= supernode_size[u_node];
558            degree[u_node] = deg;
559            degreelists[deg].push(u_node);
560            // u_id has been pushed back into degreelists
561            degree_lists_marker.unmark(u_node);
562            if (min_degree > deg)
563              min_degree = deg;
564            qxlist.pop();
565          } // while (!qxlist.empty()) {
566
567          marker.set_tag_as_multiple_tag();
568          llist.pop();
569        } // while (! llist.empty())
570
571      } // update()
572
573
574      void build_permutation(InversePermutationMap next,
575                             PermutationMap prev)
576      {
577        // collect the permutation info
578        size_type i;
579        for (i = 0; i < n; ++i) {
580          diff_t size = supernode_size[get(index_vertex_map, i)];
581          if ( size <= 0 ) {
582            prev[i] = next[i];
583            supernode_size[get(index_vertex_map, i)]
584              = next[i] + 1;  // record the supernode info
585          } else
586            prev[i] = - next[i];
587        }
588        for (i = 1; i < n + 1; ++i) {
589          if ( prev[i-1] > 0 )
590            continue;
591          diff_t parent = i;
592          while ( prev[parent - 1] < 0 ) {
593            parent = - prev[parent - 1];
594          }
595
596          diff_t root = parent;
597          diff_t num = prev[root - 1] + 1;
598          next[i-1] = - num;
599          prev[root-1] = num;
600
601          parent = i;
602          diff_t next_node = - prev[parent - 1];
603          while (next_node > 0) {
604            prev[parent-1] = - root;
605            parent = next_node;
606            next_node = - prev[parent - 1];
607          }
608        }
609        for (i = 0; i < n; i++) {
610          diff_t num = - next[i] - 1;
611          next[i] = num;
612          prev[num] = i;
613        }
614      } // build_permutation()
615    };
616
617  } //namespace detail
618
619
620  // MMD algorithm
621  //
622  //The implementation presently includes the enhancements for mass
623  //elimination, incomplete degree update, multiple elimination, and
624  //external degree.
625  //
626  //Important Note: This implementation requires the BGL graph to be
627  //directed.  Therefore, nonzero entry (i, j) in a symmetrical matrix
628  //A coresponds to two directed edges (i->j and j->i).
629  //
630  //see Alan George and Joseph W. H. Liu, The Evolution of the Minimum
631  //Degree Ordering Algorithm, SIAM Review, 31, 1989, Page 1-19
632  template<class Graph, class DegreeMap,
633           class InversePermutationMap,
634           class PermutationMap,
635           class SuperNodeMap, class VertexIndexMap>
636  void minimum_degree_ordering
637    (Graph& G,
638     DegreeMap degree,
639     InversePermutationMap inverse_perm,
640     PermutationMap perm,
641     SuperNodeMap supernode_size,
642     int delta,
643     VertexIndexMap vertex_index_map)
644  {
645    detail::mmd_impl<Graph,DegreeMap,InversePermutationMap,
646      PermutationMap, SuperNodeMap, VertexIndexMap>
647      impl(G, num_vertices(G), delta, degree, inverse_perm,
648           perm, supernode_size, vertex_index_map);
649    impl.do_mmd();
650    impl.build_permutation(inverse_perm, perm);
651  }
652
653} // namespace boost
654
655#endif // MINIMUM_DEGREE_ORDERING_HPP
Note: See TracBrowser for help on using the repository browser.