deba@1017: /* -*- mode: C++; indent-tabs-mode: nil; -*-
deba@1017:  *
deba@1017:  * This file is a part of LEMON, a generic C++ optimization library.
deba@1017:  *
alpar@1270:  * Copyright (C) 2003-2013
deba@1017:  * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
deba@1017:  * (Egervary Research Group on Combinatorial Optimization, EGRES).
deba@1017:  *
deba@1017:  * Permission to use, modify and distribute this software is granted
deba@1017:  * provided that this copyright notice appears in all copies. For
deba@1017:  * precise terms see the accompanying LICENSE file.
deba@1017:  *
deba@1017:  * This software is provided "AS IS" with no warranty of any kind,
deba@1017:  * express or implied, and with no claim as to its suitability for any
deba@1017:  * purpose.
deba@1017:  *
deba@1017:  */
deba@1017: 
deba@1017: #ifndef LEMON_NAGAMOCHI_IBARAKI_H
deba@1017: #define LEMON_NAGAMOCHI_IBARAKI_H
deba@1017: 
deba@1017: 
deba@1017: /// \ingroup min_cut
deba@1017: /// \file
deba@1017: /// \brief Implementation of the Nagamochi-Ibaraki algorithm.
deba@1017: 
deba@1017: #include <lemon/core.h>
deba@1017: #include <lemon/bin_heap.h>
deba@1017: #include <lemon/bucket_heap.h>
deba@1017: #include <lemon/maps.h>
deba@1017: #include <lemon/radix_sort.h>
deba@1017: #include <lemon/unionfind.h>
deba@1017: 
deba@1017: #include <cassert>
deba@1017: 
deba@1017: namespace lemon {
deba@1017: 
deba@1017:   /// \brief Default traits class for NagamochiIbaraki class.
deba@1017:   ///
deba@1017:   /// Default traits class for NagamochiIbaraki class.
deba@1017:   /// \param GR The undirected graph type.
deba@1017:   /// \param CM Type of capacity map.
deba@1017:   template <typename GR, typename CM>
deba@1017:   struct NagamochiIbarakiDefaultTraits {
deba@1017:     /// The type of the capacity map.
deba@1017:     typedef typename CM::Value Value;
deba@1017: 
deba@1017:     /// The undirected graph type the algorithm runs on.
deba@1017:     typedef GR Graph;
deba@1017: 
deba@1017:     /// \brief The type of the map that stores the edge capacities.
deba@1017:     ///
deba@1017:     /// The type of the map that stores the edge capacities.
deba@1017:     /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
deba@1017:     typedef CM CapacityMap;
deba@1017: 
deba@1017:     /// \brief Instantiates a CapacityMap.
deba@1017:     ///
deba@1017:     /// This function instantiates a \ref CapacityMap.
deba@1017: #ifdef DOXYGEN
deba@1017:     static CapacityMap *createCapacityMap(const Graph& graph)
deba@1017: #else
deba@1017:     static CapacityMap *createCapacityMap(const Graph&)
deba@1017: #endif
deba@1017:     {
deba@1017:         LEMON_ASSERT(false, "CapacityMap is not initialized");
deba@1017:         return 0; // ignore warnings
deba@1017:     }
deba@1017: 
deba@1017:     /// \brief The cross reference type used by heap.
deba@1017:     ///
deba@1017:     /// The cross reference type used by heap.
deba@1017:     /// Usually \c Graph::NodeMap<int>.
deba@1017:     typedef typename Graph::template NodeMap<int> HeapCrossRef;
deba@1017: 
deba@1017:     /// \brief Instantiates a HeapCrossRef.
deba@1017:     ///
deba@1017:     /// This function instantiates a \ref HeapCrossRef.
deba@1017:     /// \param g is the graph, to which we would like to define the
deba@1017:     /// \ref HeapCrossRef.
deba@1017:     static HeapCrossRef *createHeapCrossRef(const Graph& g) {
deba@1017:       return new HeapCrossRef(g);
deba@1017:     }
deba@1017: 
deba@1017:     /// \brief The heap type used by NagamochiIbaraki algorithm.
deba@1017:     ///
deba@1017:     /// The heap type used by NagamochiIbaraki algorithm. It has to
deba@1017:     /// maximize the priorities.
deba@1017:     ///
deba@1017:     /// \sa BinHeap
deba@1017:     /// \sa NagamochiIbaraki
deba@1017:     typedef BinHeap<Value, HeapCrossRef, std::greater<Value> > Heap;
deba@1017: 
deba@1017:     /// \brief Instantiates a Heap.
deba@1017:     ///
deba@1017:     /// This function instantiates a \ref Heap.
deba@1017:     /// \param r is the cross reference of the heap.
deba@1017:     static Heap *createHeap(HeapCrossRef& r) {
deba@1017:       return new Heap(r);
deba@1017:     }
deba@1017:   };
deba@1017: 
deba@1017:   /// \ingroup min_cut
deba@1017:   ///
deba@1017:   /// \brief Calculates the minimum cut in an undirected graph.
deba@1017:   ///
deba@1017:   /// Calculates the minimum cut in an undirected graph with the
deba@1017:   /// Nagamochi-Ibaraki algorithm. The algorithm separates the graph's
deba@1017:   /// nodes into two partitions with the minimum sum of edge capacities
deba@1017:   /// between the two partitions. The algorithm can be used to test
deba@1017:   /// the network reliability, especially to test how many links have
deba@1017:   /// to be destroyed in the network to split it to at least two
deba@1017:   /// distinict subnetworks.
deba@1017:   ///
deba@1017:   /// The complexity of the algorithm is \f$ O(nm\log(n)) \f$ but with
deba@1017:   /// \ref FibHeap "Fibonacci heap" it can be decreased to
deba@1017:   /// \f$ O(nm+n^2\log(n)) \f$.  When the edges have unit capacities,
deba@1017:   /// \c BucketHeap can be used which yields \f$ O(nm) \f$ time
deba@1017:   /// complexity.
deba@1017:   ///
deba@1017:   /// \warning The value type of the capacity map should be able to
deba@1017:   /// hold any cut value of the graph, otherwise the result can
deba@1017:   /// overflow.
deba@1017:   /// \note This capacity is supposed to be integer type.
deba@1017: #ifdef DOXYGEN
deba@1017:   template <typename GR, typename CM, typename TR>
deba@1017: #else
deba@1017:   template <typename GR,
deba@1017:             typename CM = typename GR::template EdgeMap<int>,
deba@1017:             typename TR = NagamochiIbarakiDefaultTraits<GR, CM> >
deba@1017: #endif
deba@1017:   class NagamochiIbaraki {
deba@1017:   public:
deba@1017: 
deba@1017:     typedef TR Traits;
deba@1017:     /// The type of the underlying graph.
deba@1017:     typedef typename Traits::Graph Graph;
deba@1017: 
deba@1017:     /// The type of the capacity map.
deba@1017:     typedef typename Traits::CapacityMap CapacityMap;
deba@1017:     /// The value type of the capacity map.
deba@1017:     typedef typename Traits::CapacityMap::Value Value;
deba@1017: 
deba@1017:     /// The heap type used by the algorithm.
deba@1017:     typedef typename Traits::Heap Heap;
deba@1017:     /// The cross reference type used for the heap.
deba@1017:     typedef typename Traits::HeapCrossRef HeapCrossRef;
deba@1017: 
deba@1017:     ///\name Named template parameters
deba@1017: 
deba@1017:     ///@{
deba@1017: 
deba@1017:     struct SetUnitCapacityTraits : public Traits {
deba@1017:       typedef ConstMap<typename Graph::Edge, Const<int, 1> > CapacityMap;
deba@1017:       static CapacityMap *createCapacityMap(const Graph&) {
deba@1017:         return new CapacityMap();
deba@1017:       }
deba@1017:     };
deba@1017: 
deba@1017:     /// \brief \ref named-templ-param "Named parameter" for setting
deba@1017:     /// the capacity map to a constMap<Edge, int, 1>() instance
deba@1017:     ///
deba@1017:     /// \ref named-templ-param "Named parameter" for setting
deba@1017:     /// the capacity map to a constMap<Edge, int, 1>() instance
deba@1017:     struct SetUnitCapacity
deba@1017:       : public NagamochiIbaraki<Graph, CapacityMap,
deba@1017:                                 SetUnitCapacityTraits> {
deba@1017:       typedef NagamochiIbaraki<Graph, CapacityMap,
deba@1017:                                SetUnitCapacityTraits> Create;
deba@1017:     };
deba@1017: 
deba@1017: 
deba@1017:     template <class H, class CR>
deba@1017:     struct SetHeapTraits : public Traits {
deba@1017:       typedef CR HeapCrossRef;
deba@1017:       typedef H Heap;
kpeter@1416:       static HeapCrossRef *createHeapCrossRef(int) {
deba@1017:         LEMON_ASSERT(false, "HeapCrossRef is not initialized");
deba@1017:         return 0; // ignore warnings
deba@1017:       }
deba@1017:       static Heap *createHeap(HeapCrossRef &) {
deba@1017:         LEMON_ASSERT(false, "Heap is not initialized");
deba@1017:         return 0; // ignore warnings
deba@1017:       }
deba@1017:     };
deba@1017: 
deba@1017:     /// \brief \ref named-templ-param "Named parameter" for setting
deba@1017:     /// heap and cross reference type
deba@1017:     ///
deba@1017:     /// \ref named-templ-param "Named parameter" for setting heap and
deba@1017:     /// cross reference type. The heap has to maximize the priorities.
deba@1017:     template <class H, class CR = RangeMap<int> >
deba@1017:     struct SetHeap
deba@1017:       : public NagamochiIbaraki<Graph, CapacityMap, SetHeapTraits<H, CR> > {
deba@1017:       typedef NagamochiIbaraki< Graph, CapacityMap, SetHeapTraits<H, CR> >
deba@1017:       Create;
deba@1017:     };
deba@1017: 
deba@1017:     template <class H, class CR>
deba@1017:     struct SetStandardHeapTraits : public Traits {
deba@1017:       typedef CR HeapCrossRef;
deba@1017:       typedef H Heap;
deba@1017:       static HeapCrossRef *createHeapCrossRef(int size) {
deba@1017:         return new HeapCrossRef(size);
deba@1017:       }
deba@1017:       static Heap *createHeap(HeapCrossRef &crossref) {
deba@1017:         return new Heap(crossref);
deba@1017:       }
deba@1017:     };
deba@1017: 
deba@1017:     /// \brief \ref named-templ-param "Named parameter" for setting
deba@1017:     /// heap and cross reference type with automatic allocation
deba@1017:     ///
deba@1017:     /// \ref named-templ-param "Named parameter" for setting heap and
deba@1017:     /// cross reference type with automatic allocation. They should
deba@1017:     /// have standard constructor interfaces to be able to
deba@1017:     /// automatically created by the algorithm (i.e. the graph should
deba@1017:     /// be passed to the constructor of the cross reference and the
deba@1017:     /// cross reference should be passed to the constructor of the
deba@1017:     /// heap). However, external heap and cross reference objects
deba@1017:     /// could also be passed to the algorithm using the \ref heap()
deba@1017:     /// function before calling \ref run() or \ref init(). The heap
deba@1017:     /// has to maximize the priorities.
deba@1017:     /// \sa SetHeap
deba@1017:     template <class H, class CR = RangeMap<int> >
deba@1017:     struct SetStandardHeap
deba@1017:       : public NagamochiIbaraki<Graph, CapacityMap,
deba@1017:                                 SetStandardHeapTraits<H, CR> > {
deba@1017:       typedef NagamochiIbaraki<Graph, CapacityMap,
deba@1017:                                SetStandardHeapTraits<H, CR> > Create;
deba@1017:     };
deba@1017: 
deba@1017:     ///@}
deba@1017: 
deba@1017: 
deba@1017:   private:
deba@1017: 
deba@1017:     const Graph &_graph;
deba@1017:     const CapacityMap *_capacity;
deba@1017:     bool _local_capacity; // unit capacity
deba@1017: 
deba@1017:     struct ArcData {
deba@1017:       typename Graph::Node target;
deba@1017:       int prev, next;
deba@1017:     };
deba@1017:     struct EdgeData {
deba@1017:       Value capacity;
deba@1017:       Value cut;
deba@1017:     };
deba@1017: 
deba@1017:     struct NodeData {
deba@1017:       int first_arc;
deba@1017:       typename Graph::Node prev, next;
deba@1017:       int curr_arc;
deba@1017:       typename Graph::Node last_rep;
deba@1017:       Value sum;
deba@1017:     };
deba@1017: 
deba@1017:     typename Graph::template NodeMap<NodeData> *_nodes;
deba@1017:     std::vector<ArcData> _arcs;
deba@1017:     std::vector<EdgeData> _edges;
deba@1017: 
deba@1017:     typename Graph::Node _first_node;
deba@1017:     int _node_num;
deba@1017: 
deba@1017:     Value _min_cut;
deba@1017: 
deba@1017:     HeapCrossRef *_heap_cross_ref;
deba@1017:     bool _local_heap_cross_ref;
deba@1017:     Heap *_heap;
deba@1017:     bool _local_heap;
deba@1017: 
deba@1017:     typedef typename Graph::template NodeMap<typename Graph::Node> NodeList;
deba@1017:     NodeList *_next_rep;
deba@1017: 
deba@1017:     typedef typename Graph::template NodeMap<bool> MinCutMap;
deba@1017:     MinCutMap *_cut_map;
deba@1017: 
deba@1017:     void createStructures() {
deba@1017:       if (!_nodes) {
deba@1017:         _nodes = new (typename Graph::template NodeMap<NodeData>)(_graph);
deba@1017:       }
deba@1017:       if (!_capacity) {
deba@1017:         _local_capacity = true;
deba@1017:         _capacity = Traits::createCapacityMap(_graph);
deba@1017:       }
deba@1017:       if (!_heap_cross_ref) {
deba@1017:         _local_heap_cross_ref = true;
deba@1017:         _heap_cross_ref = Traits::createHeapCrossRef(_graph);
deba@1017:       }
deba@1017:       if (!_heap) {
deba@1017:         _local_heap = true;
deba@1017:         _heap = Traits::createHeap(*_heap_cross_ref);
deba@1017:       }
deba@1017:       if (!_next_rep) {
deba@1017:         _next_rep = new NodeList(_graph);
deba@1017:       }
deba@1017:       if (!_cut_map) {
deba@1017:         _cut_map = new MinCutMap(_graph);
deba@1017:       }
deba@1017:     }
deba@1017: 
alpar@1130:   protected:
alpar@1130:     //This is here to avoid a gcc-3.3 compilation error.
alpar@1130:     //It should never be called.
alpar@1270:     NagamochiIbaraki() {}
alpar@1270: 
alpar@1130:   public:
deba@1017: 
deba@1017:     typedef NagamochiIbaraki Create;
deba@1017: 
deba@1017: 
deba@1017:     /// \brief Constructor.
deba@1017:     ///
deba@1017:     /// \param graph The graph the algorithm runs on.
deba@1017:     /// \param capacity The capacity map used by the algorithm.
deba@1017:     NagamochiIbaraki(const Graph& graph, const CapacityMap& capacity)
deba@1017:       : _graph(graph), _capacity(&capacity), _local_capacity(false),
deba@1017:         _nodes(0), _arcs(), _edges(), _min_cut(),
deba@1017:         _heap_cross_ref(0), _local_heap_cross_ref(false),
deba@1017:         _heap(0), _local_heap(false),
deba@1017:         _next_rep(0), _cut_map(0) {}
deba@1017: 
deba@1017:     /// \brief Constructor.
deba@1017:     ///
deba@1017:     /// This constructor can be used only when the Traits class
deba@1017:     /// defines how can the local capacity map be instantiated.
deba@1017:     /// If the SetUnitCapacity used the algorithm automatically
deba@1017:     /// constructs the capacity map.
deba@1017:     ///
deba@1017:     ///\param graph The graph the algorithm runs on.
deba@1017:     NagamochiIbaraki(const Graph& graph)
deba@1017:       : _graph(graph), _capacity(0), _local_capacity(false),
deba@1017:         _nodes(0), _arcs(), _edges(), _min_cut(),
deba@1017:         _heap_cross_ref(0), _local_heap_cross_ref(false),
deba@1017:         _heap(0), _local_heap(false),
deba@1017:         _next_rep(0), _cut_map(0) {}
deba@1017: 
deba@1017:     /// \brief Destructor.
deba@1017:     ///
deba@1017:     /// Destructor.
deba@1017:     ~NagamochiIbaraki() {
deba@1017:       if (_local_capacity) delete _capacity;
deba@1017:       if (_nodes) delete _nodes;
deba@1017:       if (_local_heap) delete _heap;
deba@1017:       if (_local_heap_cross_ref) delete _heap_cross_ref;
deba@1017:       if (_next_rep) delete _next_rep;
deba@1017:       if (_cut_map) delete _cut_map;
deba@1017:     }
deba@1017: 
deba@1017:     /// \brief Sets the heap and the cross reference used by algorithm.
deba@1017:     ///
deba@1017:     /// Sets the heap and the cross reference used by algorithm.
deba@1017:     /// If you don't use this function before calling \ref run(),
deba@1017:     /// it will allocate one. The destuctor deallocates this
deba@1017:     /// automatically allocated heap and cross reference, of course.
deba@1017:     /// \return <tt> (*this) </tt>
deba@1017:     NagamochiIbaraki &heap(Heap& hp, HeapCrossRef &cr)
deba@1017:     {
deba@1017:       if (_local_heap_cross_ref) {
deba@1017:         delete _heap_cross_ref;
deba@1017:         _local_heap_cross_ref = false;
deba@1017:       }
deba@1017:       _heap_cross_ref = &cr;
deba@1017:       if (_local_heap) {
deba@1017:         delete _heap;
deba@1017:         _local_heap = false;
deba@1017:       }
deba@1017:       _heap = &hp;
deba@1017:       return *this;
deba@1017:     }
deba@1017: 
deba@1017:     /// \name Execution control
deba@1017:     /// The simplest way to execute the algorithm is to use
deba@1017:     /// one of the member functions called \c run().
deba@1017:     /// \n
deba@1017:     /// If you need more control on the execution,
deba@1017:     /// first you must call \ref init() and then call the start()
deba@1017:     /// or proper times the processNextPhase() member functions.
deba@1017: 
deba@1017:     ///@{
deba@1017: 
deba@1017:     /// \brief Initializes the internal data structures.
deba@1017:     ///
deba@1017:     /// Initializes the internal data structures.
deba@1017:     void init() {
deba@1017:       createStructures();
deba@1017: 
deba@1017:       int edge_num = countEdges(_graph);
deba@1017:       _edges.resize(edge_num);
deba@1017:       _arcs.resize(2 * edge_num);
deba@1017: 
deba@1017:       typename Graph::Node prev = INVALID;
deba@1017:       _node_num = 0;
deba@1017:       for (typename Graph::NodeIt n(_graph); n != INVALID; ++n) {
deba@1017:         (*_cut_map)[n] = false;
deba@1017:         (*_next_rep)[n] = INVALID;
deba@1017:         (*_nodes)[n].last_rep = n;
deba@1017:         (*_nodes)[n].first_arc = -1;
deba@1017:         (*_nodes)[n].curr_arc = -1;
deba@1017:         (*_nodes)[n].prev = prev;
deba@1017:         if (prev != INVALID) {
deba@1017:           (*_nodes)[prev].next = n;
deba@1017:         }
deba@1017:         (*_nodes)[n].next = INVALID;
deba@1017:         (*_nodes)[n].sum = 0;
deba@1017:         prev = n;
deba@1017:         ++_node_num;
deba@1017:       }
deba@1017: 
deba@1017:       _first_node = typename Graph::NodeIt(_graph);
deba@1017: 
deba@1017:       int index = 0;
deba@1017:       for (typename Graph::NodeIt n(_graph); n != INVALID; ++n) {
deba@1017:         for (typename Graph::OutArcIt a(_graph, n); a != INVALID; ++a) {
deba@1017:           typename Graph::Node m = _graph.target(a);
alpar@1270: 
deba@1017:           if (!(n < m)) continue;
deba@1017: 
deba@1017:           (*_nodes)[n].sum += (*_capacity)[a];
deba@1017:           (*_nodes)[m].sum += (*_capacity)[a];
alpar@1270: 
deba@1017:           int c = (*_nodes)[m].curr_arc;
deba@1017:           if (c != -1 && _arcs[c ^ 1].target == n) {
deba@1017:             _edges[c >> 1].capacity += (*_capacity)[a];
deba@1017:           } else {
deba@1017:             _edges[index].capacity = (*_capacity)[a];
alpar@1270: 
deba@1017:             _arcs[index << 1].prev = -1;
deba@1017:             if ((*_nodes)[n].first_arc != -1) {
deba@1017:               _arcs[(*_nodes)[n].first_arc].prev = (index << 1);
deba@1017:             }
deba@1017:             _arcs[index << 1].next = (*_nodes)[n].first_arc;
deba@1017:             (*_nodes)[n].first_arc = (index << 1);
deba@1017:             _arcs[index << 1].target = m;
deba@1017: 
deba@1017:             (*_nodes)[m].curr_arc = (index << 1);
alpar@1270: 
deba@1017:             _arcs[(index << 1) | 1].prev = -1;
deba@1017:             if ((*_nodes)[m].first_arc != -1) {
deba@1017:               _arcs[(*_nodes)[m].first_arc].prev = ((index << 1) | 1);
deba@1017:             }
deba@1017:             _arcs[(index << 1) | 1].next = (*_nodes)[m].first_arc;
deba@1017:             (*_nodes)[m].first_arc = ((index << 1) | 1);
deba@1017:             _arcs[(index << 1) | 1].target = n;
alpar@1270: 
deba@1017:             ++index;
deba@1017:           }
deba@1017:         }
deba@1017:       }
deba@1017: 
deba@1017:       typename Graph::Node cut_node = INVALID;
deba@1017:       _min_cut = std::numeric_limits<Value>::max();
deba@1017: 
alpar@1270:       for (typename Graph::Node n = _first_node;
deba@1017:            n != INVALID; n = (*_nodes)[n].next) {
deba@1017:         if ((*_nodes)[n].sum < _min_cut) {
deba@1017:           cut_node = n;
deba@1017:           _min_cut = (*_nodes)[n].sum;
deba@1017:         }
deba@1017:       }
deba@1017:       (*_cut_map)[cut_node] = true;
deba@1017:       if (_min_cut == 0) {
deba@1017:         _first_node = INVALID;
deba@1017:       }
deba@1017:     }
deba@1017: 
deba@1017:   public:
deba@1017: 
deba@1017:     /// \brief Processes the next phase
deba@1017:     ///
deba@1017:     /// Processes the next phase in the algorithm. It must be called
deba@1017:     /// at most one less the number of the nodes in the graph.
deba@1017:     ///
deba@1017:     ///\return %True when the algorithm finished.
deba@1017:     bool processNextPhase() {
deba@1017:       if (_first_node == INVALID) return true;
deba@1017: 
deba@1017:       _heap->clear();
alpar@1270:       for (typename Graph::Node n = _first_node;
deba@1017:            n != INVALID; n = (*_nodes)[n].next) {
deba@1017:         (*_heap_cross_ref)[n] = Heap::PRE_HEAP;
deba@1017:       }
deba@1017: 
deba@1017:       std::vector<typename Graph::Node> order;
deba@1017:       order.reserve(_node_num);
deba@1017:       int sep = 0;
deba@1017: 
deba@1017:       Value alpha = 0;
deba@1017:       Value pmc = std::numeric_limits<Value>::max();
deba@1017: 
deba@1017:       _heap->push(_first_node, static_cast<Value>(0));
deba@1017:       while (!_heap->empty()) {
deba@1017:         typename Graph::Node n = _heap->top();
deba@1017:         Value v = _heap->prio();
deba@1017: 
deba@1017:         _heap->pop();
deba@1017:         for (int a = (*_nodes)[n].first_arc; a != -1; a = _arcs[a].next) {
deba@1017:           switch (_heap->state(_arcs[a].target)) {
alpar@1270:           case Heap::PRE_HEAP:
deba@1017:             {
deba@1017:               Value nv = _edges[a >> 1].capacity;
deba@1017:               _heap->push(_arcs[a].target, nv);
deba@1017:               _edges[a >> 1].cut = nv;
deba@1017:             } break;
deba@1017:           case Heap::IN_HEAP:
deba@1017:             {
deba@1017:               Value nv = _edges[a >> 1].capacity + (*_heap)[_arcs[a].target];
deba@1017:               _heap->decrease(_arcs[a].target, nv);
deba@1017:               _edges[a >> 1].cut = nv;
deba@1017:             } break;
deba@1017:           case Heap::POST_HEAP:
deba@1017:             break;
deba@1017:           }
deba@1017:         }
deba@1017: 
deba@1017:         alpha += (*_nodes)[n].sum;
deba@1017:         alpha -= 2 * v;
deba@1017: 
deba@1017:         order.push_back(n);
deba@1017:         if (!_heap->empty()) {
deba@1017:           if (alpha < pmc) {
deba@1017:             pmc = alpha;
deba@1017:             sep = order.size();
deba@1017:           }
deba@1017:         }
deba@1017:       }
deba@1017: 
deba@1017:       if (static_cast<int>(order.size()) < _node_num) {
deba@1017:         _first_node = INVALID;
deba@1017:         for (typename Graph::NodeIt n(_graph); n != INVALID; ++n) {
deba@1017:           (*_cut_map)[n] = false;
deba@1017:         }
deba@1017:         for (int i = 0; i < static_cast<int>(order.size()); ++i) {
deba@1017:           typename Graph::Node n = order[i];
deba@1017:           while (n != INVALID) {
deba@1017:             (*_cut_map)[n] = true;
deba@1017:             n = (*_next_rep)[n];
deba@1017:           }
deba@1017:         }
deba@1017:         _min_cut = 0;
deba@1017:         return true;
deba@1017:       }
deba@1017: 
deba@1017:       if (pmc < _min_cut) {
deba@1017:         for (typename Graph::NodeIt n(_graph); n != INVALID; ++n) {
deba@1017:           (*_cut_map)[n] = false;
deba@1017:         }
deba@1017:         for (int i = 0; i < sep; ++i) {
deba@1017:           typename Graph::Node n = order[i];
deba@1017:           while (n != INVALID) {
deba@1017:             (*_cut_map)[n] = true;
deba@1017:             n = (*_next_rep)[n];
deba@1017:           }
deba@1017:         }
deba@1017:         _min_cut = pmc;
deba@1017:       }
deba@1017: 
deba@1017:       for (typename Graph::Node n = _first_node;
deba@1017:            n != INVALID; n = (*_nodes)[n].next) {
deba@1017:         bool merged = false;
deba@1017:         for (int a = (*_nodes)[n].first_arc; a != -1; a = _arcs[a].next) {
deba@1017:           if (!(_edges[a >> 1].cut < pmc)) {
deba@1017:             if (!merged) {
deba@1017:               for (int b = (*_nodes)[n].first_arc; b != -1; b = _arcs[b].next) {
alpar@1270:                 (*_nodes)[_arcs[b].target].curr_arc = b;
deba@1017:               }
deba@1017:               merged = true;
deba@1017:             }
deba@1017:             typename Graph::Node m = _arcs[a].target;
deba@1017:             int nb = 0;
deba@1017:             for (int b = (*_nodes)[m].first_arc; b != -1; b = nb) {
deba@1017:               nb = _arcs[b].next;
deba@1017:               if ((b ^ a) == 1) continue;
deba@1017:               typename Graph::Node o = _arcs[b].target;
alpar@1270:               int c = (*_nodes)[o].curr_arc;
deba@1017:               if (c != -1 && _arcs[c ^ 1].target == n) {
deba@1017:                 _edges[c >> 1].capacity += _edges[b >> 1].capacity;
deba@1017:                 (*_nodes)[n].sum += _edges[b >> 1].capacity;
deba@1017:                 if (_edges[b >> 1].cut < _edges[c >> 1].cut) {
deba@1017:                   _edges[b >> 1].cut = _edges[c >> 1].cut;
deba@1017:                 }
deba@1017:                 if (_arcs[b ^ 1].prev != -1) {
deba@1017:                   _arcs[_arcs[b ^ 1].prev].next = _arcs[b ^ 1].next;
deba@1017:                 } else {
deba@1017:                   (*_nodes)[o].first_arc = _arcs[b ^ 1].next;
deba@1017:                 }
deba@1017:                 if (_arcs[b ^ 1].next != -1) {
deba@1017:                   _arcs[_arcs[b ^ 1].next].prev = _arcs[b ^ 1].prev;
deba@1017:                 }
deba@1017:               } else {
deba@1017:                 if (_arcs[a].next != -1) {
deba@1017:                   _arcs[_arcs[a].next].prev = b;
deba@1017:                 }
deba@1017:                 _arcs[b].next = _arcs[a].next;
deba@1017:                 _arcs[b].prev = a;
deba@1017:                 _arcs[a].next = b;
deba@1017:                 _arcs[b ^ 1].target = n;
deba@1017: 
deba@1017:                 (*_nodes)[n].sum += _edges[b >> 1].capacity;
deba@1017:                 (*_nodes)[o].curr_arc = b;
deba@1017:               }
deba@1017:             }
deba@1017: 
deba@1017:             if (_arcs[a].prev != -1) {
deba@1017:               _arcs[_arcs[a].prev].next = _arcs[a].next;
deba@1017:             } else {
deba@1017:               (*_nodes)[n].first_arc = _arcs[a].next;
alpar@1270:             }
deba@1017:             if (_arcs[a].next != -1) {
deba@1017:               _arcs[_arcs[a].next].prev = _arcs[a].prev;
deba@1017:             }
deba@1017: 
deba@1017:             (*_nodes)[n].sum -= _edges[a >> 1].capacity;
deba@1017:             (*_next_rep)[(*_nodes)[n].last_rep] = m;
deba@1017:             (*_nodes)[n].last_rep = (*_nodes)[m].last_rep;
alpar@1270: 
deba@1017:             if ((*_nodes)[m].prev != INVALID) {
deba@1017:               (*_nodes)[(*_nodes)[m].prev].next = (*_nodes)[m].next;
deba@1017:             } else{
deba@1017:               _first_node = (*_nodes)[m].next;
deba@1017:             }
deba@1017:             if ((*_nodes)[m].next != INVALID) {
deba@1017:               (*_nodes)[(*_nodes)[m].next].prev = (*_nodes)[m].prev;
deba@1017:             }
deba@1017:             --_node_num;
deba@1017:           }
deba@1017:         }
deba@1017:       }
deba@1017: 
deba@1017:       if (_node_num == 1) {
deba@1017:         _first_node = INVALID;
deba@1017:         return true;
deba@1017:       }
deba@1017: 
deba@1017:       return false;
deba@1017:     }
deba@1017: 
deba@1017:     /// \brief Executes the algorithm.
deba@1017:     ///
deba@1017:     /// Executes the algorithm.
deba@1017:     ///
deba@1017:     /// \pre init() must be called
deba@1017:     void start() {
deba@1017:       while (!processNextPhase()) {}
deba@1017:     }
deba@1017: 
deba@1017: 
deba@1017:     /// \brief Runs %NagamochiIbaraki algorithm.
deba@1017:     ///
deba@1017:     /// This method runs the %Min cut algorithm
deba@1017:     ///
deba@1017:     /// \note mc.run(s) is just a shortcut of the following code.
deba@1017:     ///\code
deba@1017:     ///  mc.init();
deba@1017:     ///  mc.start();
deba@1017:     ///\endcode
deba@1017:     void run() {
deba@1017:       init();
deba@1017:       start();
deba@1017:     }
deba@1017: 
deba@1017:     ///@}
deba@1017: 
deba@1017:     /// \name Query Functions
deba@1017:     ///
deba@1017:     /// The result of the %NagamochiIbaraki
deba@1017:     /// algorithm can be obtained using these functions.\n
deba@1017:     /// Before the use of these functions, either run() or start()
deba@1017:     /// must be called.
deba@1017: 
deba@1017:     ///@{
deba@1017: 
deba@1017:     /// \brief Returns the min cut value.
deba@1017:     ///
deba@1017:     /// Returns the min cut value if the algorithm finished.
deba@1017:     /// After the first processNextPhase() it is a value of a
deba@1017:     /// valid cut in the graph.
deba@1017:     Value minCutValue() const {
deba@1017:       return _min_cut;
deba@1017:     }
deba@1017: 
deba@1017:     /// \brief Returns a min cut in a NodeMap.
deba@1017:     ///
deba@1017:     /// It sets the nodes of one of the two partitions to true and
deba@1017:     /// the other partition to false.
deba@1017:     /// \param cutMap A \ref concepts::WriteMap "writable" node map with
deba@1017:     /// \c bool (or convertible) value type.
deba@1017:     template <typename CutMap>
deba@1017:     Value minCutMap(CutMap& cutMap) const {
deba@1017:       for (typename Graph::NodeIt n(_graph); n != INVALID; ++n) {
deba@1017:         cutMap.set(n, (*_cut_map)[n]);
deba@1017:       }
deba@1017:       return minCutValue();
deba@1017:     }
deba@1017: 
deba@1017:     ///@}
deba@1017: 
deba@1017:   };
deba@1017: }
deba@1017: 
deba@1017: #endif