lemon/nagamochi_ibaraki.h
author alpar
Tue, 24 Apr 2007 09:39:01 +0000
changeset 2436 0c941c524b47
parent 2386 81b47fc5c444
child 2506 216c6bd5c18c
permissions -rw-r--r--
Integer parameters also convert to double
     1 /* -*- C++ -*-
     2  *
     3  * This file is a part of LEMON, a generic C++ optimization library
     4  *
     5  * Copyright (C) 2003-2007
     6  * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
     7  * (Egervary Research Group on Combinatorial Optimization, EGRES).
     8  *
     9  * Permission to use, modify and distribute this software is granted
    10  * provided that this copyright notice appears in all copies. For
    11  * precise terms see the accompanying LICENSE file.
    12  *
    13  * This software is provided "AS IS" with no warranty of any kind,
    14  * express or implied, and with no claim as to its suitability for any
    15  * purpose.
    16  *
    17  */
    18 
    19 #ifndef LEMON_NAGAMOCHI_IBARAKI_H
    20 #define LEMON_NAGAMOCHI_IBARAKI_H
    21 
    22 
    23 /// \ingroup min_cut
    24 /// \file 
    25 /// \brief Maximum cardinality search and minimum cut in undirected
    26 /// graphs.
    27 
    28 #include <lemon/list_graph.h>
    29 #include <lemon/bin_heap.h>
    30 #include <lemon/bucket_heap.h>
    31 
    32 #include <lemon/unionfind.h>
    33 #include <lemon/topology.h>
    34 
    35 #include <lemon/bits/invalid.h>
    36 #include <lemon/error.h>
    37 #include <lemon/maps.h>
    38 
    39 #include <functional>
    40 
    41 #include <lemon/graph_writer.h>
    42 #include <lemon/time_measure.h>
    43 
    44 namespace lemon {
    45 
    46   namespace _min_cut_bits {
    47 
    48     template <typename CapacityMap>
    49     struct HeapSelector {
    50       template <typename Value, typename Ref>
    51       struct Selector {
    52         typedef BinHeap<Value, Ref, std::greater<Value> > Heap;
    53       };
    54     };
    55 
    56     template <typename CapacityKey>
    57     struct HeapSelector<ConstMap<CapacityKey, Const<int, 1> > > {
    58       template <typename Value, typename Ref>
    59       struct Selector {
    60         typedef BucketHeap<Ref, false > Heap;
    61       };
    62     };
    63 
    64   }
    65 
    66   /// \brief Default traits class of MaxCardinalitySearch class.
    67   ///
    68   /// Default traits class of MaxCardinalitySearch class.
    69   /// \param Graph Graph type.
    70   /// \param CapacityMap Type of length map.
    71   template <typename _Graph, typename _CapacityMap>
    72   struct MaxCardinalitySearchDefaultTraits {
    73     /// The graph type the algorithm runs on. 
    74     typedef _Graph Graph;
    75 
    76     /// \brief The type of the map that stores the edge capacities.
    77     ///
    78     /// The type of the map that stores the edge capacities.
    79     /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
    80     typedef _CapacityMap CapacityMap;
    81 
    82     /// \brief The type of the capacity of the edges.
    83     typedef typename CapacityMap::Value Value;
    84 
    85     /// \brief The cross reference type used by heap.
    86     ///
    87     /// The cross reference type used by heap.
    88     /// Usually it is \c Graph::NodeMap<int>.
    89     typedef typename Graph::template NodeMap<int> HeapCrossRef;
    90 
    91     /// \brief Instantiates a HeapCrossRef.
    92     ///
    93     /// This function instantiates a \ref HeapCrossRef. 
    94     /// \param graph is the graph, to which we would like to define the 
    95     /// HeapCrossRef.
    96     static HeapCrossRef *createHeapCrossRef(const Graph &graph) {
    97       return new HeapCrossRef(graph);
    98     }
    99     
   100     /// \brief The heap type used by MaxCardinalitySearch algorithm.
   101     ///
   102     /// The heap type used by MaxCardinalitySearch algorithm. It should
   103     /// maximalize the priorities. The default heap type is
   104     /// the \ref BinHeap, but it is specialized when the
   105     /// CapacityMap is ConstMap<Graph::Node, Const<int, 1> >
   106     /// to BucketHeap.
   107     ///
   108     /// \sa MaxCardinalitySearch
   109     typedef typename _min_cut_bits
   110     ::HeapSelector<CapacityMap>
   111     ::template Selector<Value, HeapCrossRef>
   112     ::Heap Heap;
   113 
   114     /// \brief Instantiates a Heap.
   115     ///
   116     /// This function instantiates a \ref Heap. 
   117     /// \param crossref The cross reference of the heap.
   118     static Heap *createHeap(HeapCrossRef& crossref) {
   119       return new Heap(crossref);
   120     }
   121 
   122     /// \brief The type of the map that stores whether a nodes is processed.
   123     ///
   124     /// The type of the map that stores whether a nodes is processed.
   125     /// It must meet the \ref concepts::WriteMap "WriteMap" concept.
   126     /// By default it is a NullMap.
   127     typedef NullMap<typename Graph::Node, bool> ProcessedMap;
   128 
   129     /// \brief Instantiates a ProcessedMap.
   130     ///
   131     /// This function instantiates a \ref ProcessedMap. 
   132     /// \param graph is the graph, to which
   133     /// we would like to define the \ref ProcessedMap
   134 #ifdef DOXYGEN
   135     static ProcessedMap *createProcessedMap(const Graph &graph)
   136 #else
   137     static ProcessedMap *createProcessedMap(const Graph &)
   138 #endif
   139     {
   140       return new ProcessedMap();
   141     }
   142 
   143     /// \brief The type of the map that stores the cardinalties of the nodes.
   144     /// 
   145     /// The type of the map that stores the cardinalities of the nodes.
   146     /// It must meet the \ref concepts::WriteMap "WriteMap" concept.
   147     typedef typename Graph::template NodeMap<Value> CardinalityMap;
   148 
   149     /// \brief Instantiates a CardinalityMap.
   150     ///
   151     /// This function instantiates a \ref CardinalityMap. 
   152     /// \param graph is the graph, to which we would like to define the \ref 
   153     /// CardinalityMap
   154     static CardinalityMap *createCardinalityMap(const Graph &graph) {
   155       return new CardinalityMap(graph);
   156     }
   157 
   158 
   159   };
   160   
   161   /// \ingroup search
   162   ///
   163   /// \brief Maximum Cardinality Search algorithm class.
   164   ///
   165   /// This class provides an efficient implementation of Maximum Cardinality 
   166   /// Search algorithm. The maximum cardinality search chooses first time any 
   167   /// node of the graph. Then every time it chooses the node which is connected
   168   /// to the processed nodes at most in the sum of capacities on the out 
   169   /// edges. If there is a cut in the graph the algorithm should choose
   170   /// again any unprocessed node of the graph. Each node cardinality is
   171   /// the sum of capacities on the out edges to the nodes which are processed
   172   /// before the given node.
   173   ///
   174   /// The edge capacities are passed to the algorithm using a
   175   /// \ref concepts::ReadMap "ReadMap", so it is easy to change it to any 
   176   /// kind of capacity.
   177   ///
   178   /// The type of the capacity is determined by the \ref 
   179   /// concepts::ReadMap::Value "Value" of the capacity map.
   180   ///
   181   /// It is also possible to change the underlying priority heap.
   182   ///
   183   ///
   184   /// \param _Graph The graph type the algorithm runs on. The default value
   185   /// is \ref ListGraph. The value of Graph is not used directly by
   186   /// the search algorithm, it is only passed to 
   187   /// \ref MaxCardinalitySearchDefaultTraits.
   188   /// \param _CapacityMap This read-only EdgeMap determines the capacities of 
   189   /// the edges. It is read once for each edge, so the map may involve in
   190   /// relatively time consuming process to compute the edge capacity if
   191   /// it is necessary. The default map type is \ref
   192   /// concepts::Graph::EdgeMap "Graph::EdgeMap<int>".  The value
   193   /// of CapacityMap is not used directly by search algorithm, it is only 
   194   /// passed to \ref MaxCardinalitySearchDefaultTraits.  
   195   /// \param _Traits Traits class to set various data types used by the 
   196   /// algorithm.  The default traits class is 
   197   /// \ref MaxCardinalitySearchDefaultTraits 
   198   /// "MaxCardinalitySearchDefaultTraits<_Graph, _CapacityMap>".  
   199   /// See \ref MaxCardinalitySearchDefaultTraits 
   200   /// for the documentation of a MaxCardinalitySearch traits class.
   201   ///
   202   /// \author Balazs Dezso
   203 
   204 #ifdef DOXYGEN
   205   template <typename _Graph, typename _CapacityMap, typename _Traits>
   206 #else
   207   template <typename _Graph = ListUGraph,
   208 	    typename _CapacityMap = typename _Graph::template EdgeMap<int>,
   209 	    typename _Traits = 
   210             MaxCardinalitySearchDefaultTraits<_Graph, _CapacityMap> >
   211 #endif
   212   class MaxCardinalitySearch {
   213   public:
   214     /// \brief \ref Exception for uninitialized parameters.
   215     ///
   216     /// This error represents problems in the initialization
   217     /// of the parameters of the algorithms.
   218     class UninitializedParameter : public lemon::UninitializedParameter {
   219     public:
   220       virtual const char* what() const throw() {
   221 	return "lemon::MaxCardinalitySearch::UninitializedParameter";
   222       }
   223     };
   224 
   225     typedef _Traits Traits;
   226     ///The type of the underlying graph.
   227     typedef typename Traits::Graph Graph;
   228     
   229     ///The type of the capacity of the edges.
   230     typedef typename Traits::CapacityMap::Value Value;
   231     ///The type of the map that stores the edge capacities.
   232     typedef typename Traits::CapacityMap CapacityMap;
   233     ///The type of the map indicating if a node is processed.
   234     typedef typename Traits::ProcessedMap ProcessedMap;
   235     ///The type of the map that stores the cardinalities of the nodes.
   236     typedef typename Traits::CardinalityMap CardinalityMap;
   237     ///The cross reference type used for the current heap.
   238     typedef typename Traits::HeapCrossRef HeapCrossRef;
   239     ///The heap type used by the algorithm. It maximize the priorities.
   240     typedef typename Traits::Heap Heap;
   241   private:
   242     /// Pointer to the underlying graph.
   243     const Graph *_graph;
   244     /// Pointer to the capacity map
   245     const CapacityMap *_capacity;
   246     ///Pointer to the map of cardinality.
   247     CardinalityMap *_cardinality;
   248     ///Indicates if \ref _cardinality is locally allocated (\c true) or not.
   249     bool local_cardinality;
   250     ///Pointer to the map of processed status of the nodes.
   251     ProcessedMap *_processed;
   252     ///Indicates if \ref _processed is locally allocated (\c true) or not.
   253     bool local_processed;
   254     ///Pointer to the heap cross references.
   255     HeapCrossRef *_heap_cross_ref;
   256     ///Indicates if \ref _heap_cross_ref is locally allocated (\c true) or not.
   257     bool local_heap_cross_ref;
   258     ///Pointer to the heap.
   259     Heap *_heap;
   260     ///Indicates if \ref _heap is locally allocated (\c true) or not.
   261     bool local_heap;
   262 
   263   public :
   264 
   265     typedef MaxCardinalitySearch Create;
   266  
   267     ///\name Named template parameters
   268 
   269     ///@{
   270 
   271     template <class T>
   272     struct DefCardinalityMapTraits : public Traits {
   273       typedef T CardinalityMap;
   274       static CardinalityMap *createCardinalityMap(const Graph &) 
   275       {
   276 	throw UninitializedParameter();
   277       }
   278     };
   279     /// \brief \ref named-templ-param "Named parameter" for setting 
   280     /// CardinalityMap type
   281     ///
   282     /// \ref named-templ-param "Named parameter" for setting CardinalityMap 
   283     /// type
   284     template <class T>
   285     struct DefCardinalityMap 
   286       : public MaxCardinalitySearch<Graph, CapacityMap, 
   287                                     DefCardinalityMapTraits<T> > { 
   288       typedef MaxCardinalitySearch<Graph, CapacityMap, 
   289                                    DefCardinalityMapTraits<T> > Create;
   290     };
   291     
   292     template <class T>
   293     struct DefProcessedMapTraits : public Traits {
   294       typedef T ProcessedMap;
   295       static ProcessedMap *createProcessedMap(const Graph &) {
   296 	throw UninitializedParameter();
   297       }
   298     };
   299     /// \brief \ref named-templ-param "Named parameter" for setting 
   300     /// ProcessedMap type
   301     ///
   302     /// \ref named-templ-param "Named parameter" for setting ProcessedMap type
   303     ///
   304     template <class T>
   305     struct DefProcessedMap 
   306       : public MaxCardinalitySearch<Graph, CapacityMap, 
   307                                     DefProcessedMapTraits<T> > { 
   308       typedef MaxCardinalitySearch<Graph, CapacityMap, 
   309                                    DefProcessedMapTraits<T> > Create;
   310     };
   311     
   312     template <class H, class CR>
   313     struct DefHeapTraits : public Traits {
   314       typedef CR HeapCrossRef;
   315       typedef H Heap;
   316       static HeapCrossRef *createHeapCrossRef(const Graph &) {
   317 	throw UninitializedParameter();
   318       }
   319       static Heap *createHeap(HeapCrossRef &) {
   320 	throw UninitializedParameter();
   321       }
   322     };
   323     /// \brief \ref named-templ-param "Named parameter" for setting heap 
   324     /// and cross reference type
   325     ///
   326     /// \ref named-templ-param "Named parameter" for setting heap and cross 
   327     /// reference type
   328     template <class H, class CR = typename Graph::template NodeMap<int> >
   329     struct DefHeap
   330       : public MaxCardinalitySearch<Graph, CapacityMap, 
   331                                     DefHeapTraits<H, CR> > { 
   332       typedef MaxCardinalitySearch< Graph, CapacityMap, 
   333                                     DefHeapTraits<H, CR> > Create;
   334     };
   335 
   336     template <class H, class CR>
   337     struct DefStandardHeapTraits : public Traits {
   338       typedef CR HeapCrossRef;
   339       typedef H Heap;
   340       static HeapCrossRef *createHeapCrossRef(const Graph &graph) {
   341 	return new HeapCrossRef(graph);
   342       }
   343       static Heap *createHeap(HeapCrossRef &crossref) {
   344 	return new Heap(crossref);
   345       }
   346     };
   347 
   348     /// \brief \ref named-templ-param "Named parameter" for setting heap and 
   349     /// cross reference type with automatic allocation
   350     ///
   351     /// \ref named-templ-param "Named parameter" for setting heap and cross 
   352     /// reference type. It can allocate the heap and the cross reference 
   353     /// object if the cross reference's constructor waits for the graph as 
   354     /// parameter and the heap's constructor waits for the cross reference.
   355     template <class H, class CR = typename Graph::template NodeMap<int> >
   356     struct DefStandardHeap
   357       : public MaxCardinalitySearch<Graph, CapacityMap, 
   358                                     DefStandardHeapTraits<H, CR> > { 
   359       typedef MaxCardinalitySearch<Graph, CapacityMap, 
   360                                    DefStandardHeapTraits<H, CR> > 
   361       Create;
   362     };
   363     
   364     ///@}
   365 
   366 
   367   protected:
   368 
   369     MaxCardinalitySearch() {}
   370 
   371   public:      
   372     
   373     /// \brief Constructor.
   374     ///
   375     ///\param graph the graph the algorithm will run on.
   376     ///\param capacity the capacity map used by the algorithm.
   377     MaxCardinalitySearch(const Graph& graph, const CapacityMap& capacity) :
   378       _graph(&graph), _capacity(&capacity),
   379       _cardinality(0), local_cardinality(false),
   380       _processed(0), local_processed(false),
   381       _heap_cross_ref(0), local_heap_cross_ref(false),
   382       _heap(0), local_heap(false)
   383     { }
   384     
   385     /// \brief Destructor.
   386     ~MaxCardinalitySearch() {
   387       if(local_cardinality) delete _cardinality;
   388       if(local_processed) delete _processed;
   389       if(local_heap_cross_ref) delete _heap_cross_ref;
   390       if(local_heap) delete _heap;
   391     }
   392 
   393     /// \brief Sets the capacity map.
   394     ///
   395     /// Sets the capacity map.
   396     /// \return <tt> (*this) </tt>
   397     MaxCardinalitySearch &capacityMap(const CapacityMap &m) {
   398       _capacity = &m;
   399       return *this;
   400     }
   401 
   402     /// \brief Sets the map storing the cardinalities calculated by the 
   403     /// algorithm.
   404     ///
   405     /// Sets the map storing the cardinalities calculated by the algorithm.
   406     /// If you don't use this function before calling \ref run(),
   407     /// it will allocate one. The destuctor deallocates this
   408     /// automatically allocated map, of course.
   409     /// \return <tt> (*this) </tt>
   410     MaxCardinalitySearch &cardinalityMap(CardinalityMap &m) {
   411       if(local_cardinality) {
   412 	delete _cardinality;
   413 	local_cardinality=false;
   414       }
   415       _cardinality = &m;
   416       return *this;
   417     }
   418 
   419     /// \brief Sets the map storing the processed nodes.
   420     ///
   421     /// Sets the map storing the processed nodes.
   422     /// If you don't use this function before calling \ref run(),
   423     /// it will allocate one. The destuctor deallocates this
   424     /// automatically allocated map, of course.
   425     /// \return <tt> (*this) </tt>
   426     MaxCardinalitySearch &processedMap(ProcessedMap &m) 
   427     {
   428       if(local_processed) {
   429 	delete _processed;
   430 	local_processed=false;
   431       }
   432       _processed = &m;
   433       return *this;
   434     }
   435 
   436     /// \brief Sets the heap and the cross reference used by algorithm.
   437     ///
   438     /// Sets the heap and the cross reference used by algorithm.
   439     /// If you don't use this function before calling \ref run(),
   440     /// it will allocate one. The destuctor deallocates this
   441     /// automatically allocated map, of course.
   442     /// \return <tt> (*this) </tt>
   443     MaxCardinalitySearch &heap(Heap& hp, HeapCrossRef &cr) {
   444       if(local_heap_cross_ref) {
   445 	delete _heap_cross_ref;
   446 	local_heap_cross_ref = false;
   447       }
   448       _heap_cross_ref = &cr;
   449       if(local_heap) {
   450 	delete _heap;
   451 	local_heap = false;
   452       }
   453       _heap = &hp;
   454       return *this;
   455     }
   456 
   457   private:
   458 
   459     typedef typename Graph::Node Node;
   460     typedef typename Graph::NodeIt NodeIt;
   461     typedef typename Graph::Edge Edge;
   462     typedef typename Graph::InEdgeIt InEdgeIt;
   463 
   464     void create_maps() {
   465       if(!_cardinality) {
   466 	local_cardinality = true;
   467 	_cardinality = Traits::createCardinalityMap(*_graph);
   468       }
   469       if(!_processed) {
   470 	local_processed = true;
   471 	_processed = Traits::createProcessedMap(*_graph);
   472       }
   473       if (!_heap_cross_ref) {
   474 	local_heap_cross_ref = true;
   475 	_heap_cross_ref = Traits::createHeapCrossRef(*_graph);
   476       }
   477       if (!_heap) {
   478 	local_heap = true;
   479 	_heap = Traits::createHeap(*_heap_cross_ref);
   480       }
   481     }
   482     
   483     void finalizeNodeData(Node node, Value capacity) {
   484       _processed->set(node, true);
   485       _cardinality->set(node, capacity);
   486     }
   487 
   488   public:
   489     /// \name Execution control
   490     /// The simplest way to execute the algorithm is to use
   491     /// one of the member functions called \c run(...).
   492     /// \n
   493     /// If you need more control on the execution,
   494     /// first you must call \ref init(), then you can add several source nodes
   495     /// with \ref addSource().
   496     /// Finally \ref start() will perform the actual path
   497     /// computation.
   498 
   499     ///@{
   500 
   501     /// \brief Initializes the internal data structures.
   502     ///
   503     /// Initializes the internal data structures.
   504     void init() {
   505       create_maps();
   506       _heap->clear();
   507       for (NodeIt it(*_graph) ; it != INVALID ; ++it) {
   508 	_processed->set(it, false);
   509 	_heap_cross_ref->set(it, Heap::PRE_HEAP);
   510       }
   511     }
   512     
   513     /// \brief Adds a new source node.
   514     /// 
   515     /// Adds a new source node to the priority heap.
   516     ///
   517     /// It checks if the node has not yet been added to the heap.
   518     void addSource(Node source, Value capacity = 0) {
   519       if(_heap->state(source) == Heap::PRE_HEAP) {
   520 	_heap->push(source, capacity);
   521       } 
   522     }
   523     
   524     /// \brief Processes the next node in the priority heap
   525     ///
   526     /// Processes the next node in the priority heap.
   527     ///
   528     /// \return The processed node.
   529     ///
   530     /// \warning The priority heap must not be empty!
   531     Node processNextNode() {
   532       Node node = _heap->top(); 
   533       finalizeNodeData(node, _heap->prio());
   534       _heap->pop();
   535       
   536       for (InEdgeIt it(*_graph, node); it != INVALID; ++it) {
   537 	Node source = _graph->source(it); 
   538 	switch (_heap->state(source)) {
   539 	case Heap::PRE_HEAP:
   540 	  _heap->push(source, (*_capacity)[it]); 
   541 	  break;
   542 	case Heap::IN_HEAP:
   543 	  _heap->decrease(source, (*_heap)[source] + (*_capacity)[it]); 
   544 	  break;
   545 	case Heap::POST_HEAP:
   546 	  break;
   547 	}
   548       }
   549       return node;
   550     }
   551 
   552     /// \brief Next node to be processed.
   553     ///
   554     /// Next node to be processed.
   555     ///
   556     /// \return The next node to be processed or INVALID if the 
   557     /// priority heap is empty.
   558     Node nextNode() { 
   559       return _heap->empty() ? _heap->top() : INVALID;
   560     }
   561  
   562     /// \brief Returns \c false if there are nodes
   563     /// to be processed in the priority heap
   564     ///
   565     /// Returns \c false if there are nodes
   566     /// to be processed in the priority heap
   567     bool emptyQueue() { return _heap->empty(); }
   568     /// \brief Returns the number of the nodes to be processed 
   569     /// in the priority heap
   570     ///
   571     /// Returns the number of the nodes to be processed in the priority heap
   572     int queueSize() { return _heap->size(); }
   573     
   574     /// \brief Executes the algorithm.
   575     ///
   576     /// Executes the algorithm.
   577     ///
   578     ///\pre init() must be called and at least one node should be added
   579     /// with addSource() before using this function.
   580     ///
   581     /// This method runs the Maximum Cardinality Search algorithm from the 
   582     /// source node(s).
   583     void start() {
   584       while ( !_heap->empty() ) processNextNode();
   585     }
   586     
   587     /// \brief Executes the algorithm until \c dest is reached.
   588     ///
   589     /// Executes the algorithm until \c dest is reached.
   590     ///
   591     /// \pre init() must be called and at least one node should be added
   592     /// with addSource() before using this function.
   593     ///
   594     /// This method runs the %MaxCardinalitySearch algorithm from the source 
   595     /// nodes.
   596     void start(Node dest) {
   597       while ( !_heap->empty() && _heap->top()!=dest ) processNextNode();
   598       if ( !_heap->empty() ) finalizeNodeData(_heap->top(), _heap->prio());
   599     }
   600     
   601     /// \brief Executes the algorithm until a condition is met.
   602     ///
   603     /// Executes the algorithm until a condition is met.
   604     ///
   605     /// \pre init() must be called and at least one node should be added
   606     /// with addSource() before using this function.
   607     ///
   608     /// \param nm must be a bool (or convertible) node map. The algorithm
   609     /// will stop when it reaches a node \c v with <tt>nm[v]==true</tt>.
   610     template <typename NodeBoolMap>
   611     void start(const NodeBoolMap &nm) {
   612       while ( !_heap->empty() && !nm[_heap->top()] ) processNextNode();
   613       if ( !_heap->empty() ) finalizeNodeData(_heap->top(),_heap->prio());
   614     }
   615     
   616     /// \brief Runs the maximal cardinality search algorithm from node \c s.
   617     ///
   618     /// This method runs the %MaxCardinalitySearch algorithm from a root 
   619     /// node \c s.
   620     ///
   621     ///\note d.run(s) is just a shortcut of the following code.
   622     ///\code
   623     ///  d.init();
   624     ///  d.addSource(s);
   625     ///  d.start();
   626     ///\endcode
   627     void run(Node s) {
   628       init();
   629       addSource(s);
   630       start();
   631     }
   632 
   633     /// \brief Runs the maximal cardinality search algorithm for the 
   634     /// whole graph.
   635     ///
   636     /// This method runs the %MaxCardinalitySearch algorithm from all 
   637     /// unprocessed node of the graph.
   638     ///
   639     ///\note d.run(s) is just a shortcut of the following code.
   640     ///\code
   641     ///  d.init();
   642     ///  for (NodeIt it(graph); it != INVALID; ++it) {
   643     ///    if (!d.reached(it)) {
   644     ///      d.addSource(s);
   645     ///      d.start();
   646     ///    }
   647     ///  }
   648     ///\endcode
   649     void run() {
   650       init();
   651       for (NodeIt it(*_graph); it != INVALID; ++it) {
   652         if (!reached(it)) {
   653           addSource(it);
   654           start();
   655         }
   656       }
   657     }
   658     
   659     ///@}
   660 
   661     /// \name Query Functions
   662     /// The result of the maximum cardinality search algorithm can be 
   663     /// obtained using these functions.
   664     /// \n
   665     /// Before the use of these functions, either run() or start() must be 
   666     /// called.
   667     
   668     ///@{
   669 
   670     /// \brief The cardinality of a node.
   671     ///
   672     /// Returns the cardinality of a node.
   673     /// \pre \ref run() must be called before using this function.
   674     /// \warning If node \c v in unreachable from the root the return value
   675     /// of this funcion is undefined.
   676     Value cardinality(Node node) const { return (*_cardinality)[node]; }
   677 
   678     /// \brief The current cardinality of a node.
   679     ///
   680     /// Returns the current cardinality of a node.
   681     /// \pre the given node should be reached but not processed
   682     Value currentCardinality(Node node) const { return (*_heap)[node]; }
   683 
   684     /// \brief Returns a reference to the NodeMap of cardinalities.
   685     ///
   686     /// Returns a reference to the NodeMap of cardinalities. \pre \ref run() 
   687     /// must be called before using this function.
   688     const CardinalityMap &cardinalityMap() const { return *_cardinality;}
   689  
   690     /// \brief Checks if a node is reachable from the root.
   691     ///
   692     /// Returns \c true if \c v is reachable from the root.
   693     /// \warning The source nodes are inditated as unreached.
   694     /// \pre \ref run() must be called before using this function.
   695     bool reached(Node v) { return (*_heap_cross_ref)[v] != Heap::PRE_HEAP; }
   696 
   697     /// \brief Checks if a node is processed.
   698     ///
   699     /// Returns \c true if \c v is processed, i.e. the shortest
   700     /// path to \c v has already found.
   701     /// \pre \ref run() must be called before using this function.
   702     bool processed(Node v) { return (*_heap_cross_ref)[v] == Heap::POST_HEAP; }
   703     
   704     ///@}
   705   };
   706 
   707   /// \brief Default traits class of NagamochiIbaraki class.
   708   ///
   709   /// Default traits class of NagamochiIbaraki class.
   710   /// \param Graph Graph type.
   711   /// \param CapacityMap Type of length map.
   712   template <typename _Graph, typename _CapacityMap>
   713   struct NagamochiIbarakiDefaultTraits {
   714     /// \brief The type of the capacity of the edges.
   715     typedef typename _CapacityMap::Value Value;
   716 
   717     /// The graph type the algorithm runs on. 
   718     typedef _Graph Graph;
   719 
   720     /// The AuxGraph type which is an Graph
   721     typedef ListUGraph AuxGraph;
   722 
   723     /// \brief Instantiates a AuxGraph.
   724     ///
   725     /// This function instantiates a \ref AuxGraph. 
   726     static AuxGraph *createAuxGraph() {
   727       return new AuxGraph();
   728     }
   729 
   730     /// \brief The type of the map that stores the edge capacities.
   731     ///
   732     /// The type of the map that stores the edge capacities.
   733     /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
   734     typedef _CapacityMap CapacityMap;
   735 
   736     /// \brief Instantiates a CapacityMap.
   737     ///
   738     /// This function instantiates a \ref CapacityMap.
   739 #ifdef DOXYGEN
   740     static CapacityMap *createCapacityMap(const Graph& graph) 
   741 #else
   742     static CapacityMap *createCapacityMap(const Graph&)
   743 #endif
   744     {
   745       throw UninitializedParameter();
   746     }
   747 
   748     /// \brief The CutValueMap type
   749     ///
   750     /// The type of the map that stores the cut value of a node.
   751     typedef AuxGraph::NodeMap<Value> AuxCutValueMap;
   752 
   753     /// \brief Instantiates a AuxCutValueMap.
   754     ///
   755     /// This function instantiates a \ref AuxCutValueMap. 
   756     static AuxCutValueMap *createAuxCutValueMap(const AuxGraph& graph) {
   757       return new AuxCutValueMap(graph);
   758     }
   759 
   760     /// \brief The AuxCapacityMap type
   761     ///
   762     /// The type of the map that stores the auxiliary edge capacities.
   763     typedef AuxGraph::UEdgeMap<Value> AuxCapacityMap;
   764 
   765     /// \brief Instantiates a AuxCapacityMap.
   766     ///
   767     /// This function instantiates a \ref AuxCapacityMap. 
   768     static AuxCapacityMap *createAuxCapacityMap(const AuxGraph& graph) {
   769       return new AuxCapacityMap(graph);
   770     }
   771 
   772     /// \brief The cross reference type used by heap.
   773     ///
   774     /// The cross reference type used by heap.
   775     /// Usually it is \c Graph::NodeMap<int>.
   776     typedef AuxGraph::NodeMap<int> HeapCrossRef;
   777 
   778     /// \brief Instantiates a HeapCrossRef.
   779     ///
   780     /// This function instantiates a \ref HeapCrossRef. 
   781     /// \param graph is the graph, to which we would like to define the 
   782     /// HeapCrossRef.
   783     static HeapCrossRef *createHeapCrossRef(const AuxGraph &graph) {
   784       return new HeapCrossRef(graph);
   785     }
   786     
   787     /// \brief The heap type used by NagamochiIbaraki algorithm.
   788     ///
   789     /// The heap type used by NagamochiIbaraki algorithm. It should
   790     /// maximalize the priorities and the heap's key type is
   791     /// the aux graph's node.
   792     ///
   793     /// \sa BinHeap
   794     /// \sa NagamochiIbaraki
   795     typedef typename _min_cut_bits
   796     ::HeapSelector<CapacityMap>
   797     ::template Selector<Value, HeapCrossRef>
   798     ::Heap Heap;
   799     
   800     /// \brief Instantiates a Heap.
   801     ///
   802     /// This function instantiates a \ref Heap. 
   803     /// \param crossref The cross reference of the heap.
   804     static Heap *createHeap(HeapCrossRef& crossref) {
   805       return new Heap(crossref);
   806     }
   807 
   808     /// \brief Map from the AuxGraph's node type to the Graph's node type.
   809     ///
   810     /// Map from the AuxGraph's node type to the Graph's node type.
   811     typedef typename AuxGraph
   812     ::template NodeMap<typename Graph::Node> NodeRefMap;
   813 
   814     /// \brief Instantiates a NodeRefMap.
   815     ///
   816     /// This function instantiates a \ref NodeRefMap. 
   817     static NodeRefMap *createNodeRefMap(const AuxGraph& graph) {
   818       return new NodeRefMap(graph);
   819     }
   820 
   821     /// \brief Map from the Graph's node type to the Graph's node type.
   822     ///
   823     /// Map from the Graph's node type to the Graph's node type.
   824     typedef typename Graph
   825     ::template NodeMap<typename Graph::Node> ListRefMap;
   826 
   827     /// \brief Instantiates a ListRefMap.
   828     ///
   829     /// This function instantiates a \ref ListRefMap. 
   830     static ListRefMap *createListRefMap(const Graph& graph) {
   831       return new ListRefMap(graph);
   832     }
   833     
   834 
   835   };
   836 
   837   /// \ingroup min_cut
   838   ///
   839   /// \brief Calculates the minimum cut in an undirected graph.
   840   ///
   841   /// Calculates the minimum cut in an undirected graph with the
   842   /// Nagamochi-Ibaraki algorithm. The algorithm separates the graph's
   843   /// nodes into two partitions with the minimum sum of edge capacities
   844   /// between the two partitions. The algorithm can be used to test
   845   /// the network reliability specifically to test how many links have
   846   /// to be destroyed in the network to split it at least two
   847   /// distinict subnetwork.
   848   ///
   849   /// The complexity of the algorithm is \f$ O(ne\log(n)) \f$ but with
   850   /// Fibonacci heap it can be decreased to \f$ O(ne+n^2\log(n)) \f$. 
   851   /// When capacity map is neutral then it uses BucketHeap which
   852   /// results \f$ O(ne) \f$ time complexity.
   853   ///
   854   /// \warning The value type of the capacity map should be able to hold
   855   /// any cut value of the graph, otherwise the result can overflow.
   856 #ifdef DOXYGEN
   857   template <typename _Graph, typename _CapacityMap, typename _Traits>
   858 #else
   859   template <typename _Graph = ListUGraph, 
   860 	    typename _CapacityMap = typename _Graph::template UEdgeMap<int>, 
   861 	    typename _Traits 
   862             = NagamochiIbarakiDefaultTraits<_Graph, _CapacityMap> >
   863 #endif
   864   class NagamochiIbaraki {
   865   public:
   866     /// \brief \ref Exception for uninitialized parameters.
   867     ///
   868     /// This error represents problems in the initialization
   869     /// of the parameters of the algorithms.
   870     class UninitializedParameter : public lemon::UninitializedParameter {
   871     public:
   872       virtual const char* what() const throw() {
   873 	return "lemon::NagamochiIbaraki::UninitializedParameter";
   874       }
   875     };
   876 
   877 
   878   private:
   879 
   880     typedef _Traits Traits;
   881     /// The type of the underlying graph.
   882     typedef typename Traits::Graph Graph;
   883     
   884     /// The type of the capacity of the edges.
   885     typedef typename Traits::CapacityMap::Value Value;
   886     /// The type of the map that stores the edge capacities.
   887     typedef typename Traits::CapacityMap CapacityMap;
   888     /// The type of the aux graph
   889     typedef typename Traits::AuxGraph AuxGraph;
   890     /// The type of the aux capacity map
   891     typedef typename Traits::AuxCapacityMap AuxCapacityMap;
   892     /// The type of the aux cut value map
   893     typedef typename Traits::AuxCutValueMap AuxCutValueMap;
   894     /// The cross reference type used for the current heap.
   895     typedef typename Traits::HeapCrossRef HeapCrossRef;
   896     /// The heap type used by the max cardinality algorithm.
   897     typedef typename Traits::Heap Heap;
   898     /// The node refrefernces between the original and aux graph type.
   899     typedef typename Traits::NodeRefMap NodeRefMap;
   900     /// The list node refrefernces in the original graph type.
   901     typedef typename Traits::ListRefMap ListRefMap;
   902 
   903   public:
   904 
   905     ///\name Named template parameters
   906 
   907     ///@{
   908 
   909     struct DefNeutralCapacityTraits : public Traits {
   910       typedef ConstMap<typename Graph::UEdge, Const<int, 1> > CapacityMap;
   911       static CapacityMap *createCapacityMap(const Graph&) {
   912 	return new CapacityMap();
   913       }
   914     };
   915     /// \brief \ref named-templ-param "Named parameter" for setting  
   916     /// the capacity type to constMap<UEdge, int, 1>()
   917     ///
   918     /// \ref named-templ-param "Named parameter" for setting 
   919     /// the capacity type to constMap<UEdge, int, 1>()
   920     struct DefNeutralCapacity
   921       : public NagamochiIbaraki<Graph, CapacityMap, 
   922                                 DefNeutralCapacityTraits> { 
   923       typedef NagamochiIbaraki<Graph, CapacityMap, 
   924                                DefNeutralCapacityTraits> Create;
   925     };
   926 
   927 
   928     template <class H, class CR>
   929     struct DefHeapTraits : public Traits {
   930       typedef CR HeapCrossRef;
   931       typedef H Heap;
   932       static HeapCrossRef *createHeapCrossRef(const AuxGraph &) {
   933 	throw UninitializedParameter();
   934       }
   935       static Heap *createHeap(HeapCrossRef &) {
   936 	throw UninitializedParameter();
   937       }
   938     };
   939     /// \brief \ref named-templ-param "Named parameter" for setting heap 
   940     /// and cross reference type
   941     ///
   942     /// \ref named-templ-param "Named parameter" for setting heap and cross 
   943     /// reference type
   944     template <class H, class CR = typename Graph::template NodeMap<int> >
   945     struct DefHeap
   946       : public NagamochiIbaraki<Graph, CapacityMap, 
   947                                 DefHeapTraits<H, CR> > { 
   948       typedef NagamochiIbaraki< Graph, CapacityMap, 
   949                                 DefHeapTraits<H, CR> > Create;
   950     };
   951 
   952     template <class H, class CR>
   953     struct DefStandardHeapTraits : public Traits {
   954       typedef CR HeapCrossRef;
   955       typedef H Heap;
   956       static HeapCrossRef *createHeapCrossRef(const AuxGraph &graph) {
   957 	return new HeapCrossRef(graph);
   958       }
   959       static Heap *createHeap(HeapCrossRef &crossref) {
   960 	return new Heap(crossref);
   961       }
   962     };
   963 
   964     /// \brief \ref named-templ-param "Named parameter" for setting heap and 
   965     /// cross reference type with automatic allocation
   966     ///
   967     /// \ref named-templ-param "Named parameter" for setting heap and cross 
   968     /// reference type. It can allocate the heap and the cross reference 
   969     /// object if the cross reference's constructor waits for the graph as 
   970     /// parameter and the heap's constructor waits for the cross reference.
   971     template <class H, class CR = typename Graph::template NodeMap<int> >
   972     struct DefStandardHeap
   973       : public NagamochiIbaraki<Graph, CapacityMap, 
   974                                 DefStandardHeapTraits<H, CR> > { 
   975       typedef NagamochiIbaraki<Graph, CapacityMap, 
   976                                DefStandardHeapTraits<H, CR> > 
   977       Create;
   978     };
   979 
   980     ///@}
   981 
   982 
   983   private:
   984     /// Pointer to the underlying graph.
   985     const Graph *_graph;
   986     /// Pointer to the capacity map
   987     const CapacityMap *_capacity;
   988     /// \brief Indicates if \ref _capacity is locally allocated 
   989     /// (\c true) or not.
   990     bool local_capacity;
   991 
   992     /// Pointer to the aux graph.
   993     AuxGraph *_aux_graph;
   994     /// \brief Indicates if \ref _aux_graph is locally allocated 
   995     /// (\c true) or not.
   996     bool local_aux_graph;
   997     /// Pointer to the aux capacity map
   998     AuxCapacityMap *_aux_capacity;
   999     /// \brief Indicates if \ref _aux_capacity is locally allocated 
  1000     /// (\c true) or not.
  1001     bool local_aux_capacity;
  1002     /// Pointer to the aux cut value map
  1003     AuxCutValueMap *_aux_cut_value;
  1004     /// \brief Indicates if \ref _aux_cut_value is locally allocated 
  1005     /// (\c true) or not.
  1006     bool local_aux_cut_value;
  1007     /// Pointer to the heap cross references.
  1008     HeapCrossRef *_heap_cross_ref;
  1009     /// \brief Indicates if \ref _heap_cross_ref is locally allocated 
  1010     /// (\c true) or not.
  1011     bool local_heap_cross_ref;
  1012     /// Pointer to the heap.
  1013     Heap *_heap;
  1014     /// Indicates if \ref _heap is locally allocated (\c true) or not.
  1015     bool local_heap;
  1016 
  1017     /// The min cut value.
  1018     Value _min_cut;
  1019     /// The number of the nodes of the aux graph.
  1020     int _node_num;
  1021     /// The first and last node of the min cut in the next list.
  1022     std::vector<typename Graph::Node> _cut;
  1023 
  1024     /// \brief The first and last element in the list associated
  1025     /// to the aux graph node.
  1026     NodeRefMap *_first, *_last;
  1027     /// \brief The next node in the node lists.
  1028     ListRefMap *_next;
  1029     
  1030     void createStructures() {
  1031       if (!_capacity) {
  1032         local_capacity = true;
  1033         _capacity = Traits::createCapacityMap(*_graph);
  1034       }
  1035       if(!_aux_graph) {
  1036 	local_aux_graph = true;
  1037 	_aux_graph = Traits::createAuxGraph();
  1038       }
  1039       if(!_aux_capacity) {
  1040 	local_aux_capacity = true;
  1041 	_aux_capacity = Traits::createAuxCapacityMap(*_aux_graph);
  1042       }
  1043       if(!_aux_cut_value) {
  1044 	local_aux_cut_value = true;
  1045 	_aux_cut_value = Traits::createAuxCutValueMap(*_aux_graph);
  1046       }
  1047 
  1048       _first = Traits::createNodeRefMap(*_aux_graph);
  1049       _last = Traits::createNodeRefMap(*_aux_graph);
  1050 
  1051       _next = Traits::createListRefMap(*_graph);
  1052 
  1053       if (!_heap_cross_ref) {
  1054 	local_heap_cross_ref = true;
  1055 	_heap_cross_ref = Traits::createHeapCrossRef(*_aux_graph);
  1056       }
  1057       if (!_heap) {
  1058 	local_heap = true;
  1059 	_heap = Traits::createHeap(*_heap_cross_ref);
  1060       }
  1061     }
  1062 
  1063     void createAuxGraph() {
  1064       typename Graph::template NodeMap<typename AuxGraph::Node> ref(*_graph);
  1065 
  1066       for (typename Graph::NodeIt n(*_graph); n != INVALID; ++n) {
  1067         _next->set(n, INVALID);
  1068         typename AuxGraph::Node node = _aux_graph->addNode();
  1069         _first->set(node, n);
  1070         _last->set(node, n);
  1071         ref.set(n, node);
  1072       }
  1073 
  1074       typename AuxGraph::template NodeMap<typename AuxGraph::UEdge> 
  1075       edges(*_aux_graph, INVALID);
  1076 
  1077       for (typename Graph::NodeIt n(*_graph); n != INVALID; ++n) {
  1078         for (typename Graph::IncEdgeIt e(*_graph, n); e != INVALID; ++e) {
  1079           typename Graph::Node tn = _graph->runningNode(e);
  1080           if (n < tn || n == tn) continue;
  1081           if (edges[ref[tn]] != INVALID) {
  1082             Value value = 
  1083               (*_aux_capacity)[edges[ref[tn]]] + (*_capacity)[e];
  1084             _aux_capacity->set(edges[ref[tn]], value);
  1085           } else {
  1086             edges.set(ref[tn], _aux_graph->addEdge(ref[n], ref[tn]));
  1087             Value value = (*_capacity)[e];
  1088             _aux_capacity->set(edges[ref[tn]], value);            
  1089           }
  1090         }
  1091         for (typename Graph::IncEdgeIt e(*_graph, n); e != INVALID; ++e) {
  1092           typename Graph::Node tn = _graph->runningNode(e);
  1093           edges.set(ref[tn], INVALID);
  1094         }
  1095       }
  1096 
  1097       _cut.resize(1, INVALID);
  1098       _min_cut = std::numeric_limits<Value>::max();
  1099       for (typename AuxGraph::NodeIt n(*_aux_graph); n != INVALID; ++n) {
  1100         Value value = 0;
  1101         for (typename AuxGraph::IncEdgeIt e(*_aux_graph, n); 
  1102              e != INVALID; ++e) {
  1103           value += (*_aux_capacity)[e];
  1104         }
  1105         if (_min_cut > value) {
  1106           _min_cut = value;
  1107           _cut[0] = (*_first)[n];
  1108         } 
  1109         (*_aux_cut_value)[n] = value;
  1110       }
  1111     }
  1112     
  1113 
  1114   public :
  1115 
  1116     typedef NagamochiIbaraki Create;
  1117 
  1118 
  1119     /// \brief Constructor.
  1120     ///
  1121     ///\param graph the graph the algorithm will run on.
  1122     ///\param capacity the capacity map used by the algorithm.
  1123     NagamochiIbaraki(const Graph& graph, const CapacityMap& capacity) 
  1124       : _graph(&graph), 
  1125         _capacity(&capacity), local_capacity(false),
  1126         _aux_graph(0), local_aux_graph(false),
  1127         _aux_capacity(0), local_aux_capacity(false),
  1128         _aux_cut_value(0), local_aux_cut_value(false),
  1129         _heap_cross_ref(0), local_heap_cross_ref(false),
  1130         _heap(0), local_heap(false),
  1131         _first(0), _last(0), _next(0) {}
  1132 
  1133     /// \brief Constructor.
  1134     ///
  1135     /// This constructor can be used only when the Traits class
  1136     /// defines how can we instantiate a local capacity map.
  1137     /// If the DefNeutralCapacity used the algorithm automatically
  1138     /// construct the capacity map.
  1139     ///
  1140     ///\param graph the graph the algorithm will run on.
  1141     NagamochiIbaraki(const Graph& graph) 
  1142       : _graph(&graph), 
  1143         _capacity(0), local_capacity(false),
  1144         _aux_graph(0), local_aux_graph(false),
  1145         _aux_capacity(0), local_aux_capacity(false),
  1146         _aux_cut_value(0), local_aux_cut_value(false),
  1147         _heap_cross_ref(0), local_heap_cross_ref(false),
  1148         _heap(0), local_heap(false),
  1149         _first(0), _last(0), _next(0) {}
  1150 
  1151     /// \brief Destructor.
  1152     ///
  1153     /// Destructor.
  1154     ~NagamochiIbaraki() {
  1155       if (local_heap) delete _heap;
  1156       if (local_heap_cross_ref) delete _heap_cross_ref;
  1157       if (_first) delete _first;
  1158       if (_last) delete _last;
  1159       if (_next) delete _next;
  1160       if (local_aux_capacity) delete _aux_capacity;
  1161       if (local_aux_cut_value) delete _aux_cut_value;
  1162       if (local_aux_graph) delete _aux_graph;
  1163       if (local_capacity) delete _capacity;
  1164     }
  1165 
  1166     /// \brief Sets the heap and the cross reference used by algorithm.
  1167     ///
  1168     /// Sets the heap and the cross reference used by algorithm.
  1169     /// If you don't use this function before calling \ref run(),
  1170     /// it will allocate one. The destuctor deallocates this
  1171     /// automatically allocated heap and cross reference, of course.
  1172     /// \return <tt> (*this) </tt>
  1173     NagamochiIbaraki &heap(Heap& hp, HeapCrossRef &cr)
  1174     {
  1175       if (local_heap_cross_ref) {
  1176 	delete _heap_cross_ref;
  1177 	local_heap_cross_ref=false;
  1178       }
  1179       _heap_cross_ref = &cr;
  1180       if (local_heap) {
  1181 	delete _heap;
  1182 	local_heap=false;
  1183       }
  1184       _heap = &hp;
  1185       return *this;
  1186     }
  1187 
  1188     /// \brief Sets the aux graph.
  1189     ///
  1190     /// Sets the aux graph used by algorithm.
  1191     /// If you don't use this function before calling \ref run(),
  1192     /// it will allocate one. The destuctor deallocates this
  1193     /// automatically allocated graph, of course.
  1194     /// \return <tt> (*this) </tt>
  1195     NagamochiIbaraki &auxGraph(AuxGraph& aux_graph)
  1196     {
  1197       if(local_aux_graph) {
  1198 	delete _aux_graph;
  1199 	local_aux_graph=false;
  1200       }
  1201       _aux_graph = &aux_graph;
  1202       return *this;
  1203     }
  1204 
  1205     /// \brief Sets the aux capacity map.
  1206     ///
  1207     /// Sets the aux capacity map used by algorithm.
  1208     /// If you don't use this function before calling \ref run(),
  1209     /// it will allocate one. The destuctor deallocates this
  1210     /// automatically allocated graph, of course.
  1211     /// \return <tt> (*this) </tt>
  1212     NagamochiIbaraki &auxCapacityMap(AuxCapacityMap& aux_capacity_map)
  1213     {
  1214       if(local_aux_capacity) {
  1215 	delete _aux_capacity;
  1216 	local_aux_capacity=false;
  1217       }
  1218       _aux_capacity = &aux_capacity_map;
  1219       return *this;
  1220     }
  1221 
  1222     /// \name Execution control
  1223     /// The simplest way to execute the algorithm is to use
  1224     /// one of the member functions called \c run().
  1225     /// \n
  1226     /// If you need more control on the execution,
  1227     /// first you must call \ref init() and then call the start()
  1228     /// or proper times the processNextPhase() member functions.
  1229 
  1230     ///@{
  1231 
  1232     /// \brief Initializes the internal data structures.
  1233     ///
  1234     /// Initializes the internal data structures.
  1235     void init() {
  1236       _node_num = countNodes(*_graph);
  1237       createStructures();
  1238       createAuxGraph();
  1239     }
  1240 
  1241   private:
  1242 
  1243     struct EdgeInfo {
  1244       typename AuxGraph::Node source, target;
  1245       Value capacity;
  1246     };
  1247     
  1248     struct EdgeInfoLess {
  1249       bool operator()(const EdgeInfo& left, const EdgeInfo& right) const {
  1250         return (left.source < right.source) || 
  1251           (left.source == right.source && left.target < right.target);
  1252       }
  1253     };
  1254 
  1255   public:
  1256 
  1257 
  1258     /// \brief Processes the next phase
  1259     ///
  1260     /// Processes the next phase in the algorithm. The function should
  1261     /// be called at most countNodes(graph) - 1 times to get surely
  1262     /// the min cut in the graph.
  1263     ///
  1264     ///\return %True when the algorithm finished.
  1265     bool processNextPhase() {
  1266       if (_node_num <= 1) return true;
  1267 
  1268       typedef typename AuxGraph::Node Node;
  1269       typedef typename AuxGraph::NodeIt NodeIt;
  1270       typedef typename AuxGraph::UEdge UEdge;
  1271       typedef typename AuxGraph::UEdgeIt UEdgeIt;
  1272       typedef typename AuxGraph::IncEdgeIt IncEdgeIt;
  1273       
  1274       typename AuxGraph::template UEdgeMap<Value> _edge_cut(*_aux_graph);
  1275 
  1276 
  1277       for (NodeIt n(*_aux_graph); n != INVALID; ++n) {
  1278         _heap_cross_ref->set(n, Heap::PRE_HEAP);
  1279       }
  1280 
  1281       std::vector<Node> nodes;
  1282       nodes.reserve(_node_num);
  1283       int sep = 0;
  1284 
  1285       Value alpha = 0;
  1286       Value emc = std::numeric_limits<Value>::max();
  1287 
  1288       _heap->push(typename AuxGraph::NodeIt(*_aux_graph), 0);
  1289       while (!_heap->empty()) {
  1290         Node node = _heap->top();
  1291         Value value = _heap->prio();
  1292         
  1293         _heap->pop();
  1294         for (typename AuxGraph::IncEdgeIt e(*_aux_graph, node); 
  1295              e != INVALID; ++e) {
  1296           Node tn = _aux_graph->runningNode(e);
  1297           switch (_heap->state(tn)) {
  1298           case Heap::PRE_HEAP:
  1299             _heap->push(tn, (*_aux_capacity)[e]);
  1300             _edge_cut[e] = (*_heap)[tn];
  1301             break;
  1302           case Heap::IN_HEAP:
  1303             _heap->decrease(tn, (*_aux_capacity)[e] + (*_heap)[tn]);
  1304             _edge_cut[e] = (*_heap)[tn];
  1305             break;
  1306           case Heap::POST_HEAP:
  1307             break;
  1308           }
  1309         }
  1310 
  1311         alpha += (*_aux_cut_value)[node];
  1312         alpha -= 2 * value;
  1313 
  1314         nodes.push_back(node);
  1315         if (!_heap->empty()) {
  1316           if (alpha < emc) {
  1317             emc = alpha;
  1318             sep = nodes.size();
  1319           }
  1320         }
  1321       }
  1322 
  1323       if (int(nodes.size()) < _node_num) {
  1324         _aux_graph->clear();
  1325         _node_num = 1;
  1326         _cut.clear();
  1327         for (int i = 0; i < int(nodes.size()); ++i) {
  1328           typename Graph::Node n = (*_first)[nodes[i]];
  1329           while (n != INVALID) {
  1330             _cut.push_back(n);
  1331             n = (*_next)[n];
  1332           }
  1333         }
  1334         _min_cut = 0;
  1335         return true;
  1336       }
  1337 
  1338       if (emc < _min_cut) {
  1339         _cut.clear();
  1340         for (int i = 0; i < sep; ++i) {
  1341           typename Graph::Node n = (*_first)[nodes[i]];
  1342           while (n != INVALID) {
  1343             _cut.push_back(n);
  1344             n = (*_next)[n];
  1345           }
  1346         }
  1347         _min_cut = emc;
  1348       }
  1349 
  1350       typedef typename AuxGraph::template NodeMap<int> UfeCr;
  1351       UfeCr ufecr(*_aux_graph);
  1352       typedef UnionFindEnum<UfeCr> Ufe; 
  1353       Ufe ufe(ufecr);
  1354 
  1355       for (typename AuxGraph::NodeIt n(*_aux_graph); n != INVALID; ++n) {
  1356         ufe.insert(n);
  1357       }
  1358 
  1359       for (typename AuxGraph::UEdgeIt e(*_aux_graph); e != INVALID; ++e) {
  1360         if (_edge_cut[e] >= emc) {
  1361           ufe.join(_aux_graph->source(e), _aux_graph->target(e));
  1362         }
  1363       }
  1364 
  1365       typedef typename Ufe::ClassIt UfeCIt;
  1366       if (ufe.size(UfeCIt(ufe)) == _node_num) {
  1367         _aux_graph->clear();
  1368         _node_num = 1;
  1369         return true;
  1370       }
  1371       
  1372       std::vector<typename AuxGraph::Node> remnodes;
  1373 
  1374       typename AuxGraph::template NodeMap<UEdge> edges(*_aux_graph, INVALID);
  1375       for (typename Ufe::ClassIt c(ufe); c != INVALID; ++c) {
  1376         if (ufe.size(c) == 1) continue;
  1377         for (typename Ufe::ItemIt r(ufe, c); r != INVALID; ++r) {
  1378           if (static_cast<Node>(r) == static_cast<Node>(c)) continue;
  1379           _next->set((*_last)[c], (*_first)[r]);
  1380           _last->set(c, (*_last)[r]);
  1381           remnodes.push_back(r);
  1382           --_node_num;
  1383         }
  1384       }
  1385 
  1386       std::vector<EdgeInfo> addedges;
  1387       std::vector<UEdge> remedges;
  1388 
  1389       for (typename AuxGraph::UEdgeIt e(*_aux_graph);
  1390            e != INVALID; ++e) {
  1391         Node sn = ufe.find(_aux_graph->source(e));
  1392         Node tn = ufe.find(_aux_graph->target(e));
  1393         if ((ufe.size(sn) == 1 && ufe.size(tn) == 1)) {
  1394           continue;
  1395         }
  1396         if (sn == tn) {
  1397           remedges.push_back(e);
  1398           continue;
  1399         }
  1400         EdgeInfo info;
  1401         if (sn < tn) {
  1402           info.source = sn;
  1403           info.target = tn;
  1404         } else {
  1405           info.source = tn;
  1406           info.target = sn;
  1407         }
  1408         info.capacity = (*_aux_capacity)[e];
  1409         addedges.push_back(info);
  1410         remedges.push_back(e);
  1411       }
  1412 
  1413       for (int i = 0; i < int(remedges.size()); ++i) {
  1414         _aux_graph->erase(remedges[i]);
  1415       }
  1416 
  1417       sort(addedges.begin(), addedges.end(), EdgeInfoLess());
  1418 
  1419       {
  1420         int i = 0;
  1421         while (i < int(addedges.size())) {
  1422           Node sn = addedges[i].source;
  1423           Node tn = addedges[i].target;
  1424           Value ec = addedges[i].capacity;
  1425           ++i;
  1426           while (i < int(addedges.size()) && 
  1427                  sn == addedges[i].source && tn == addedges[i].target) {
  1428             ec += addedges[i].capacity;
  1429             ++i;
  1430           }
  1431           typename AuxGraph::UEdge ne = _aux_graph->addEdge(sn, tn);
  1432           (*_aux_capacity)[ne] = ec;
  1433         }
  1434       }
  1435 
  1436       for (typename Ufe::ClassIt c(ufe); c != INVALID; ++c) {
  1437         if (ufe.size(c) == 1) continue;
  1438         Value cutvalue = 0;
  1439         for (typename AuxGraph::IncEdgeIt e(*_aux_graph, c);
  1440              e != INVALID; ++e) {
  1441           cutvalue += (*_aux_capacity)[e];
  1442         }
  1443         
  1444         (*_aux_cut_value)[c] = cutvalue;
  1445         
  1446       }
  1447 
  1448       for (int i = 0; i < int(remnodes.size()); ++i) {
  1449         _aux_graph->erase(remnodes[i]);
  1450       }
  1451 
  1452       return _node_num == 1;
  1453     }
  1454 
  1455     /// \brief Executes the algorithm.
  1456     ///
  1457     /// Executes the algorithm.
  1458     ///
  1459     /// \pre init() must be called
  1460     void start() {
  1461       while (!processNextPhase());
  1462     }
  1463 
  1464 
  1465     /// \brief Runs %NagamochiIbaraki algorithm.
  1466     ///
  1467     /// This method runs the %Min cut algorithm
  1468     ///
  1469     /// \note mc.run(s) is just a shortcut of the following code.
  1470     ///\code
  1471     ///  mc.init();
  1472     ///  mc.start();
  1473     ///\endcode
  1474     void run() {
  1475       init();
  1476       start();
  1477     }
  1478 
  1479     ///@}
  1480 
  1481     /// \name Query Functions 
  1482     ///
  1483     /// The result of the %NagamochiIbaraki
  1484     /// algorithm can be obtained using these functions.\n 
  1485     /// Before the use of these functions, either run() or start()
  1486     /// must be called.
  1487     
  1488     ///@{
  1489 
  1490     /// \brief Returns the min cut value.
  1491     ///
  1492     /// Returns the min cut value if the algorithm finished.
  1493     /// After the first processNextPhase() it is a value of a
  1494     /// valid cut in the graph.
  1495     Value minCut() const {
  1496       return _min_cut;
  1497     }
  1498 
  1499     /// \brief Returns a min cut in a NodeMap.
  1500     ///
  1501     /// It sets the nodes of one of the two partitions to true in
  1502     /// the given BoolNodeMap. The map contains a valid cut if the
  1503     /// map have been set false previously. 
  1504     template <typename NodeMap>
  1505     Value quickMinCut(NodeMap& nodeMap) const { 
  1506       for (int i = 0; i < int(_cut.size()); ++i) {
  1507         nodeMap.set(_cut[i], true);
  1508       }
  1509       return minCut();
  1510     }
  1511 
  1512     /// \brief Returns a min cut in a NodeMap.
  1513     ///
  1514     /// It sets the nodes of one of the two partitions to true and
  1515     /// the other partition to false. The function first set all of the
  1516     /// nodes to false and after it call the quickMinCut() member.
  1517     template <typename NodeMap>
  1518     Value minCut(NodeMap& nodeMap) const { 
  1519       for (typename Graph::NodeIt it(*_graph); it != INVALID; ++it) {
  1520         nodeMap.set(it, false);      
  1521       }
  1522       quickMinCut(nodeMap);
  1523       return minCut();
  1524     }
  1525 
  1526     /// \brief Returns a min cut in an EdgeMap.
  1527     ///
  1528     /// If an undirected edge is in a min cut then it will be
  1529     /// set to true and the others will be set to false in the given map.
  1530     template <typename EdgeMap>
  1531     Value cutEdges(EdgeMap& edgeMap) const {
  1532       typename Graph::template NodeMap<bool> cut(*_graph, false);
  1533       quickMinCut(cut);
  1534       for (typename Graph::EdgeIt it(*_graph); it != INVALID; ++it) {
  1535         edgeMap.set(it, cut[_graph->source(it)] ^ cut[_graph->target(it)]);
  1536       }
  1537       return minCut();
  1538     }
  1539 
  1540     ///@}
  1541 
  1542   private:
  1543 
  1544 
  1545   };
  1546 }
  1547 
  1548 #endif