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