COIN-OR::LEMON - Graph Library

source: lemon-0.x/lemon/nagamochi_ibaraki.h @ 2354:3609c77b77be

Last change on this file since 2354:3609c77b77be was 2337:9c3d44ac39fb, checked in by Balazs Dezso, 17 years ago

Adding two heuristics
Based on:
http://www.avglab.com/andrew/pub/neci-tr-96-132.ps

File size: 49.4 KB
Line 
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
42namespace 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
Note: See TracBrowser for help on using the repository browser.