COIN-OR::LEMON - Graph Library

source: lemon-0.x/lemon/nagamochi_ibaraki.h @ 2579:691ce54544c5

Last change on this file since 2579:691ce54544c5 was 2553:bfced05fa852, checked in by Alpar Juttner, 16 years ago

Happy New Year to LEMON (+ better update-copyright-header script)

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