lemon/nagamochi_ibaraki.h
changeset 1184 3c00344f49c9
parent 978 eb252f805431
child 1197 f179aa1045a4
     1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/lemon/nagamochi_ibaraki.h	Wed Oct 17 19:14:07 2018 +0200
     1.3 @@ -0,0 +1,702 @@
     1.4 +/* -*- mode: C++; indent-tabs-mode: nil; -*-
     1.5 + *
     1.6 + * This file is a part of LEMON, a generic C++ optimization library.
     1.7 + *
     1.8 + * Copyright (C) 2003-2013
     1.9 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
    1.10 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
    1.11 + *
    1.12 + * Permission to use, modify and distribute this software is granted
    1.13 + * provided that this copyright notice appears in all copies. For
    1.14 + * precise terms see the accompanying LICENSE file.
    1.15 + *
    1.16 + * This software is provided "AS IS" with no warranty of any kind,
    1.17 + * express or implied, and with no claim as to its suitability for any
    1.18 + * purpose.
    1.19 + *
    1.20 + */
    1.21 +
    1.22 +#ifndef LEMON_NAGAMOCHI_IBARAKI_H
    1.23 +#define LEMON_NAGAMOCHI_IBARAKI_H
    1.24 +
    1.25 +
    1.26 +/// \ingroup min_cut
    1.27 +/// \file
    1.28 +/// \brief Implementation of the Nagamochi-Ibaraki algorithm.
    1.29 +
    1.30 +#include <lemon/core.h>
    1.31 +#include <lemon/bin_heap.h>
    1.32 +#include <lemon/bucket_heap.h>
    1.33 +#include <lemon/maps.h>
    1.34 +#include <lemon/radix_sort.h>
    1.35 +#include <lemon/unionfind.h>
    1.36 +
    1.37 +#include <cassert>
    1.38 +
    1.39 +namespace lemon {
    1.40 +
    1.41 +  /// \brief Default traits class for NagamochiIbaraki class.
    1.42 +  ///
    1.43 +  /// Default traits class for NagamochiIbaraki class.
    1.44 +  /// \param GR The undirected graph type.
    1.45 +  /// \param CM Type of capacity map.
    1.46 +  template <typename GR, typename CM>
    1.47 +  struct NagamochiIbarakiDefaultTraits {
    1.48 +    /// The type of the capacity map.
    1.49 +    typedef typename CM::Value Value;
    1.50 +
    1.51 +    /// The undirected graph type the algorithm runs on.
    1.52 +    typedef GR Graph;
    1.53 +
    1.54 +    /// \brief The type of the map that stores the edge capacities.
    1.55 +    ///
    1.56 +    /// The type of the map that stores the edge capacities.
    1.57 +    /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
    1.58 +    typedef CM CapacityMap;
    1.59 +
    1.60 +    /// \brief Instantiates a CapacityMap.
    1.61 +    ///
    1.62 +    /// This function instantiates a \ref CapacityMap.
    1.63 +#ifdef DOXYGEN
    1.64 +    static CapacityMap *createCapacityMap(const Graph& graph)
    1.65 +#else
    1.66 +    static CapacityMap *createCapacityMap(const Graph&)
    1.67 +#endif
    1.68 +    {
    1.69 +        LEMON_ASSERT(false, "CapacityMap is not initialized");
    1.70 +        return 0; // ignore warnings
    1.71 +    }
    1.72 +
    1.73 +    /// \brief The cross reference type used by heap.
    1.74 +    ///
    1.75 +    /// The cross reference type used by heap.
    1.76 +    /// Usually \c Graph::NodeMap<int>.
    1.77 +    typedef typename Graph::template NodeMap<int> HeapCrossRef;
    1.78 +
    1.79 +    /// \brief Instantiates a HeapCrossRef.
    1.80 +    ///
    1.81 +    /// This function instantiates a \ref HeapCrossRef.
    1.82 +    /// \param g is the graph, to which we would like to define the
    1.83 +    /// \ref HeapCrossRef.
    1.84 +    static HeapCrossRef *createHeapCrossRef(const Graph& g) {
    1.85 +      return new HeapCrossRef(g);
    1.86 +    }
    1.87 +
    1.88 +    /// \brief The heap type used by NagamochiIbaraki algorithm.
    1.89 +    ///
    1.90 +    /// The heap type used by NagamochiIbaraki algorithm. It has to
    1.91 +    /// maximize the priorities.
    1.92 +    ///
    1.93 +    /// \sa BinHeap
    1.94 +    /// \sa NagamochiIbaraki
    1.95 +    typedef BinHeap<Value, HeapCrossRef, std::greater<Value> > Heap;
    1.96 +
    1.97 +    /// \brief Instantiates a Heap.
    1.98 +    ///
    1.99 +    /// This function instantiates a \ref Heap.
   1.100 +    /// \param r is the cross reference of the heap.
   1.101 +    static Heap *createHeap(HeapCrossRef& r) {
   1.102 +      return new Heap(r);
   1.103 +    }
   1.104 +  };
   1.105 +
   1.106 +  /// \ingroup min_cut
   1.107 +  ///
   1.108 +  /// \brief Calculates the minimum cut in an undirected graph.
   1.109 +  ///
   1.110 +  /// Calculates the minimum cut in an undirected graph with the
   1.111 +  /// Nagamochi-Ibaraki algorithm. The algorithm separates the graph's
   1.112 +  /// nodes into two partitions with the minimum sum of edge capacities
   1.113 +  /// between the two partitions. The algorithm can be used to test
   1.114 +  /// the network reliability, especially to test how many links have
   1.115 +  /// to be destroyed in the network to split it to at least two
   1.116 +  /// distinict subnetworks.
   1.117 +  ///
   1.118 +  /// The complexity of the algorithm is \f$ O(nm\log(n)) \f$ but with
   1.119 +  /// \ref FibHeap "Fibonacci heap" it can be decreased to
   1.120 +  /// \f$ O(nm+n^2\log(n)) \f$.  When the edges have unit capacities,
   1.121 +  /// \c BucketHeap can be used which yields \f$ O(nm) \f$ time
   1.122 +  /// complexity.
   1.123 +  ///
   1.124 +  /// \warning The value type of the capacity map should be able to
   1.125 +  /// hold any cut value of the graph, otherwise the result can
   1.126 +  /// overflow.
   1.127 +  /// \note This capacity is supposed to be integer type.
   1.128 +#ifdef DOXYGEN
   1.129 +  template <typename GR, typename CM, typename TR>
   1.130 +#else
   1.131 +  template <typename GR,
   1.132 +            typename CM = typename GR::template EdgeMap<int>,
   1.133 +            typename TR = NagamochiIbarakiDefaultTraits<GR, CM> >
   1.134 +#endif
   1.135 +  class NagamochiIbaraki {
   1.136 +  public:
   1.137 +
   1.138 +    typedef TR Traits;
   1.139 +    /// The type of the underlying graph.
   1.140 +    typedef typename Traits::Graph Graph;
   1.141 +
   1.142 +    /// The type of the capacity map.
   1.143 +    typedef typename Traits::CapacityMap CapacityMap;
   1.144 +    /// The value type of the capacity map.
   1.145 +    typedef typename Traits::CapacityMap::Value Value;
   1.146 +
   1.147 +    /// The heap type used by the algorithm.
   1.148 +    typedef typename Traits::Heap Heap;
   1.149 +    /// The cross reference type used for the heap.
   1.150 +    typedef typename Traits::HeapCrossRef HeapCrossRef;
   1.151 +
   1.152 +    ///\name Named template parameters
   1.153 +
   1.154 +    ///@{
   1.155 +
   1.156 +    struct SetUnitCapacityTraits : public Traits {
   1.157 +      typedef ConstMap<typename Graph::Edge, Const<int, 1> > CapacityMap;
   1.158 +      static CapacityMap *createCapacityMap(const Graph&) {
   1.159 +        return new CapacityMap();
   1.160 +      }
   1.161 +    };
   1.162 +
   1.163 +    /// \brief \ref named-templ-param "Named parameter" for setting
   1.164 +    /// the capacity map to a constMap<Edge, int, 1>() instance
   1.165 +    ///
   1.166 +    /// \ref named-templ-param "Named parameter" for setting
   1.167 +    /// the capacity map to a constMap<Edge, int, 1>() instance
   1.168 +    struct SetUnitCapacity
   1.169 +      : public NagamochiIbaraki<Graph, CapacityMap,
   1.170 +                                SetUnitCapacityTraits> {
   1.171 +      typedef NagamochiIbaraki<Graph, CapacityMap,
   1.172 +                               SetUnitCapacityTraits> Create;
   1.173 +    };
   1.174 +
   1.175 +
   1.176 +    template <class H, class CR>
   1.177 +    struct SetHeapTraits : public Traits {
   1.178 +      typedef CR HeapCrossRef;
   1.179 +      typedef H Heap;
   1.180 +      static HeapCrossRef *createHeapCrossRef(int num) {
   1.181 +        LEMON_ASSERT(false, "HeapCrossRef is not initialized");
   1.182 +        return 0; // ignore warnings
   1.183 +      }
   1.184 +      static Heap *createHeap(HeapCrossRef &) {
   1.185 +        LEMON_ASSERT(false, "Heap is not initialized");
   1.186 +        return 0; // ignore warnings
   1.187 +      }
   1.188 +    };
   1.189 +
   1.190 +    /// \brief \ref named-templ-param "Named parameter" for setting
   1.191 +    /// heap and cross reference type
   1.192 +    ///
   1.193 +    /// \ref named-templ-param "Named parameter" for setting heap and
   1.194 +    /// cross reference type. The heap has to maximize the priorities.
   1.195 +    template <class H, class CR = RangeMap<int> >
   1.196 +    struct SetHeap
   1.197 +      : public NagamochiIbaraki<Graph, CapacityMap, SetHeapTraits<H, CR> > {
   1.198 +      typedef NagamochiIbaraki< Graph, CapacityMap, SetHeapTraits<H, CR> >
   1.199 +      Create;
   1.200 +    };
   1.201 +
   1.202 +    template <class H, class CR>
   1.203 +    struct SetStandardHeapTraits : public Traits {
   1.204 +      typedef CR HeapCrossRef;
   1.205 +      typedef H Heap;
   1.206 +      static HeapCrossRef *createHeapCrossRef(int size) {
   1.207 +        return new HeapCrossRef(size);
   1.208 +      }
   1.209 +      static Heap *createHeap(HeapCrossRef &crossref) {
   1.210 +        return new Heap(crossref);
   1.211 +      }
   1.212 +    };
   1.213 +
   1.214 +    /// \brief \ref named-templ-param "Named parameter" for setting
   1.215 +    /// heap and cross reference type with automatic allocation
   1.216 +    ///
   1.217 +    /// \ref named-templ-param "Named parameter" for setting heap and
   1.218 +    /// cross reference type with automatic allocation. They should
   1.219 +    /// have standard constructor interfaces to be able to
   1.220 +    /// automatically created by the algorithm (i.e. the graph should
   1.221 +    /// be passed to the constructor of the cross reference and the
   1.222 +    /// cross reference should be passed to the constructor of the
   1.223 +    /// heap). However, external heap and cross reference objects
   1.224 +    /// could also be passed to the algorithm using the \ref heap()
   1.225 +    /// function before calling \ref run() or \ref init(). The heap
   1.226 +    /// has to maximize the priorities.
   1.227 +    /// \sa SetHeap
   1.228 +    template <class H, class CR = RangeMap<int> >
   1.229 +    struct SetStandardHeap
   1.230 +      : public NagamochiIbaraki<Graph, CapacityMap,
   1.231 +                                SetStandardHeapTraits<H, CR> > {
   1.232 +      typedef NagamochiIbaraki<Graph, CapacityMap,
   1.233 +                               SetStandardHeapTraits<H, CR> > Create;
   1.234 +    };
   1.235 +
   1.236 +    ///@}
   1.237 +
   1.238 +
   1.239 +  private:
   1.240 +
   1.241 +    const Graph &_graph;
   1.242 +    const CapacityMap *_capacity;
   1.243 +    bool _local_capacity; // unit capacity
   1.244 +
   1.245 +    struct ArcData {
   1.246 +      typename Graph::Node target;
   1.247 +      int prev, next;
   1.248 +    };
   1.249 +    struct EdgeData {
   1.250 +      Value capacity;
   1.251 +      Value cut;
   1.252 +    };
   1.253 +
   1.254 +    struct NodeData {
   1.255 +      int first_arc;
   1.256 +      typename Graph::Node prev, next;
   1.257 +      int curr_arc;
   1.258 +      typename Graph::Node last_rep;
   1.259 +      Value sum;
   1.260 +    };
   1.261 +
   1.262 +    typename Graph::template NodeMap<NodeData> *_nodes;
   1.263 +    std::vector<ArcData> _arcs;
   1.264 +    std::vector<EdgeData> _edges;
   1.265 +
   1.266 +    typename Graph::Node _first_node;
   1.267 +    int _node_num;
   1.268 +
   1.269 +    Value _min_cut;
   1.270 +
   1.271 +    HeapCrossRef *_heap_cross_ref;
   1.272 +    bool _local_heap_cross_ref;
   1.273 +    Heap *_heap;
   1.274 +    bool _local_heap;
   1.275 +
   1.276 +    typedef typename Graph::template NodeMap<typename Graph::Node> NodeList;
   1.277 +    NodeList *_next_rep;
   1.278 +
   1.279 +    typedef typename Graph::template NodeMap<bool> MinCutMap;
   1.280 +    MinCutMap *_cut_map;
   1.281 +
   1.282 +    void createStructures() {
   1.283 +      if (!_nodes) {
   1.284 +        _nodes = new (typename Graph::template NodeMap<NodeData>)(_graph);
   1.285 +      }
   1.286 +      if (!_capacity) {
   1.287 +        _local_capacity = true;
   1.288 +        _capacity = Traits::createCapacityMap(_graph);
   1.289 +      }
   1.290 +      if (!_heap_cross_ref) {
   1.291 +        _local_heap_cross_ref = true;
   1.292 +        _heap_cross_ref = Traits::createHeapCrossRef(_graph);
   1.293 +      }
   1.294 +      if (!_heap) {
   1.295 +        _local_heap = true;
   1.296 +        _heap = Traits::createHeap(*_heap_cross_ref);
   1.297 +      }
   1.298 +      if (!_next_rep) {
   1.299 +        _next_rep = new NodeList(_graph);
   1.300 +      }
   1.301 +      if (!_cut_map) {
   1.302 +        _cut_map = new MinCutMap(_graph);
   1.303 +      }
   1.304 +    }
   1.305 +
   1.306 +  protected:
   1.307 +    //This is here to avoid a gcc-3.3 compilation error.
   1.308 +    //It should never be called.
   1.309 +    NagamochiIbaraki() {}
   1.310 +
   1.311 +  public:
   1.312 +
   1.313 +    typedef NagamochiIbaraki Create;
   1.314 +
   1.315 +
   1.316 +    /// \brief Constructor.
   1.317 +    ///
   1.318 +    /// \param graph The graph the algorithm runs on.
   1.319 +    /// \param capacity The capacity map used by the algorithm.
   1.320 +    NagamochiIbaraki(const Graph& graph, const CapacityMap& capacity)
   1.321 +      : _graph(graph), _capacity(&capacity), _local_capacity(false),
   1.322 +        _nodes(0), _arcs(), _edges(), _min_cut(),
   1.323 +        _heap_cross_ref(0), _local_heap_cross_ref(false),
   1.324 +        _heap(0), _local_heap(false),
   1.325 +        _next_rep(0), _cut_map(0) {}
   1.326 +
   1.327 +    /// \brief Constructor.
   1.328 +    ///
   1.329 +    /// This constructor can be used only when the Traits class
   1.330 +    /// defines how can the local capacity map be instantiated.
   1.331 +    /// If the SetUnitCapacity used the algorithm automatically
   1.332 +    /// constructs the capacity map.
   1.333 +    ///
   1.334 +    ///\param graph The graph the algorithm runs on.
   1.335 +    NagamochiIbaraki(const Graph& graph)
   1.336 +      : _graph(graph), _capacity(0), _local_capacity(false),
   1.337 +        _nodes(0), _arcs(), _edges(), _min_cut(),
   1.338 +        _heap_cross_ref(0), _local_heap_cross_ref(false),
   1.339 +        _heap(0), _local_heap(false),
   1.340 +        _next_rep(0), _cut_map(0) {}
   1.341 +
   1.342 +    /// \brief Destructor.
   1.343 +    ///
   1.344 +    /// Destructor.
   1.345 +    ~NagamochiIbaraki() {
   1.346 +      if (_local_capacity) delete _capacity;
   1.347 +      if (_nodes) delete _nodes;
   1.348 +      if (_local_heap) delete _heap;
   1.349 +      if (_local_heap_cross_ref) delete _heap_cross_ref;
   1.350 +      if (_next_rep) delete _next_rep;
   1.351 +      if (_cut_map) delete _cut_map;
   1.352 +    }
   1.353 +
   1.354 +    /// \brief Sets the heap and the cross reference used by algorithm.
   1.355 +    ///
   1.356 +    /// Sets the heap and the cross reference used by algorithm.
   1.357 +    /// If you don't use this function before calling \ref run(),
   1.358 +    /// it will allocate one. The destuctor deallocates this
   1.359 +    /// automatically allocated heap and cross reference, of course.
   1.360 +    /// \return <tt> (*this) </tt>
   1.361 +    NagamochiIbaraki &heap(Heap& hp, HeapCrossRef &cr)
   1.362 +    {
   1.363 +      if (_local_heap_cross_ref) {
   1.364 +        delete _heap_cross_ref;
   1.365 +        _local_heap_cross_ref = false;
   1.366 +      }
   1.367 +      _heap_cross_ref = &cr;
   1.368 +      if (_local_heap) {
   1.369 +        delete _heap;
   1.370 +        _local_heap = false;
   1.371 +      }
   1.372 +      _heap = &hp;
   1.373 +      return *this;
   1.374 +    }
   1.375 +
   1.376 +    /// \name Execution control
   1.377 +    /// The simplest way to execute the algorithm is to use
   1.378 +    /// one of the member functions called \c run().
   1.379 +    /// \n
   1.380 +    /// If you need more control on the execution,
   1.381 +    /// first you must call \ref init() and then call the start()
   1.382 +    /// or proper times the processNextPhase() member functions.
   1.383 +
   1.384 +    ///@{
   1.385 +
   1.386 +    /// \brief Initializes the internal data structures.
   1.387 +    ///
   1.388 +    /// Initializes the internal data structures.
   1.389 +    void init() {
   1.390 +      createStructures();
   1.391 +
   1.392 +      int edge_num = countEdges(_graph);
   1.393 +      _edges.resize(edge_num);
   1.394 +      _arcs.resize(2 * edge_num);
   1.395 +
   1.396 +      typename Graph::Node prev = INVALID;
   1.397 +      _node_num = 0;
   1.398 +      for (typename Graph::NodeIt n(_graph); n != INVALID; ++n) {
   1.399 +        (*_cut_map)[n] = false;
   1.400 +        (*_next_rep)[n] = INVALID;
   1.401 +        (*_nodes)[n].last_rep = n;
   1.402 +        (*_nodes)[n].first_arc = -1;
   1.403 +        (*_nodes)[n].curr_arc = -1;
   1.404 +        (*_nodes)[n].prev = prev;
   1.405 +        if (prev != INVALID) {
   1.406 +          (*_nodes)[prev].next = n;
   1.407 +        }
   1.408 +        (*_nodes)[n].next = INVALID;
   1.409 +        (*_nodes)[n].sum = 0;
   1.410 +        prev = n;
   1.411 +        ++_node_num;
   1.412 +      }
   1.413 +
   1.414 +      _first_node = typename Graph::NodeIt(_graph);
   1.415 +
   1.416 +      int index = 0;
   1.417 +      for (typename Graph::NodeIt n(_graph); n != INVALID; ++n) {
   1.418 +        for (typename Graph::OutArcIt a(_graph, n); a != INVALID; ++a) {
   1.419 +          typename Graph::Node m = _graph.target(a);
   1.420 +
   1.421 +          if (!(n < m)) continue;
   1.422 +
   1.423 +          (*_nodes)[n].sum += (*_capacity)[a];
   1.424 +          (*_nodes)[m].sum += (*_capacity)[a];
   1.425 +
   1.426 +          int c = (*_nodes)[m].curr_arc;
   1.427 +          if (c != -1 && _arcs[c ^ 1].target == n) {
   1.428 +            _edges[c >> 1].capacity += (*_capacity)[a];
   1.429 +          } else {
   1.430 +            _edges[index].capacity = (*_capacity)[a];
   1.431 +
   1.432 +            _arcs[index << 1].prev = -1;
   1.433 +            if ((*_nodes)[n].first_arc != -1) {
   1.434 +              _arcs[(*_nodes)[n].first_arc].prev = (index << 1);
   1.435 +            }
   1.436 +            _arcs[index << 1].next = (*_nodes)[n].first_arc;
   1.437 +            (*_nodes)[n].first_arc = (index << 1);
   1.438 +            _arcs[index << 1].target = m;
   1.439 +
   1.440 +            (*_nodes)[m].curr_arc = (index << 1);
   1.441 +
   1.442 +            _arcs[(index << 1) | 1].prev = -1;
   1.443 +            if ((*_nodes)[m].first_arc != -1) {
   1.444 +              _arcs[(*_nodes)[m].first_arc].prev = ((index << 1) | 1);
   1.445 +            }
   1.446 +            _arcs[(index << 1) | 1].next = (*_nodes)[m].first_arc;
   1.447 +            (*_nodes)[m].first_arc = ((index << 1) | 1);
   1.448 +            _arcs[(index << 1) | 1].target = n;
   1.449 +
   1.450 +            ++index;
   1.451 +          }
   1.452 +        }
   1.453 +      }
   1.454 +
   1.455 +      typename Graph::Node cut_node = INVALID;
   1.456 +      _min_cut = std::numeric_limits<Value>::max();
   1.457 +
   1.458 +      for (typename Graph::Node n = _first_node;
   1.459 +           n != INVALID; n = (*_nodes)[n].next) {
   1.460 +        if ((*_nodes)[n].sum < _min_cut) {
   1.461 +          cut_node = n;
   1.462 +          _min_cut = (*_nodes)[n].sum;
   1.463 +        }
   1.464 +      }
   1.465 +      (*_cut_map)[cut_node] = true;
   1.466 +      if (_min_cut == 0) {
   1.467 +        _first_node = INVALID;
   1.468 +      }
   1.469 +    }
   1.470 +
   1.471 +  public:
   1.472 +
   1.473 +    /// \brief Processes the next phase
   1.474 +    ///
   1.475 +    /// Processes the next phase in the algorithm. It must be called
   1.476 +    /// at most one less the number of the nodes in the graph.
   1.477 +    ///
   1.478 +    ///\return %True when the algorithm finished.
   1.479 +    bool processNextPhase() {
   1.480 +      if (_first_node == INVALID) return true;
   1.481 +
   1.482 +      _heap->clear();
   1.483 +      for (typename Graph::Node n = _first_node;
   1.484 +           n != INVALID; n = (*_nodes)[n].next) {
   1.485 +        (*_heap_cross_ref)[n] = Heap::PRE_HEAP;
   1.486 +      }
   1.487 +
   1.488 +      std::vector<typename Graph::Node> order;
   1.489 +      order.reserve(_node_num);
   1.490 +      int sep = 0;
   1.491 +
   1.492 +      Value alpha = 0;
   1.493 +      Value pmc = std::numeric_limits<Value>::max();
   1.494 +
   1.495 +      _heap->push(_first_node, static_cast<Value>(0));
   1.496 +      while (!_heap->empty()) {
   1.497 +        typename Graph::Node n = _heap->top();
   1.498 +        Value v = _heap->prio();
   1.499 +
   1.500 +        _heap->pop();
   1.501 +        for (int a = (*_nodes)[n].first_arc; a != -1; a = _arcs[a].next) {
   1.502 +          switch (_heap->state(_arcs[a].target)) {
   1.503 +          case Heap::PRE_HEAP:
   1.504 +            {
   1.505 +              Value nv = _edges[a >> 1].capacity;
   1.506 +              _heap->push(_arcs[a].target, nv);
   1.507 +              _edges[a >> 1].cut = nv;
   1.508 +            } break;
   1.509 +          case Heap::IN_HEAP:
   1.510 +            {
   1.511 +              Value nv = _edges[a >> 1].capacity + (*_heap)[_arcs[a].target];
   1.512 +              _heap->decrease(_arcs[a].target, nv);
   1.513 +              _edges[a >> 1].cut = nv;
   1.514 +            } break;
   1.515 +          case Heap::POST_HEAP:
   1.516 +            break;
   1.517 +          }
   1.518 +        }
   1.519 +
   1.520 +        alpha += (*_nodes)[n].sum;
   1.521 +        alpha -= 2 * v;
   1.522 +
   1.523 +        order.push_back(n);
   1.524 +        if (!_heap->empty()) {
   1.525 +          if (alpha < pmc) {
   1.526 +            pmc = alpha;
   1.527 +            sep = order.size();
   1.528 +          }
   1.529 +        }
   1.530 +      }
   1.531 +
   1.532 +      if (static_cast<int>(order.size()) < _node_num) {
   1.533 +        _first_node = INVALID;
   1.534 +        for (typename Graph::NodeIt n(_graph); n != INVALID; ++n) {
   1.535 +          (*_cut_map)[n] = false;
   1.536 +        }
   1.537 +        for (int i = 0; i < static_cast<int>(order.size()); ++i) {
   1.538 +          typename Graph::Node n = order[i];
   1.539 +          while (n != INVALID) {
   1.540 +            (*_cut_map)[n] = true;
   1.541 +            n = (*_next_rep)[n];
   1.542 +          }
   1.543 +        }
   1.544 +        _min_cut = 0;
   1.545 +        return true;
   1.546 +      }
   1.547 +
   1.548 +      if (pmc < _min_cut) {
   1.549 +        for (typename Graph::NodeIt n(_graph); n != INVALID; ++n) {
   1.550 +          (*_cut_map)[n] = false;
   1.551 +        }
   1.552 +        for (int i = 0; i < sep; ++i) {
   1.553 +          typename Graph::Node n = order[i];
   1.554 +          while (n != INVALID) {
   1.555 +            (*_cut_map)[n] = true;
   1.556 +            n = (*_next_rep)[n];
   1.557 +          }
   1.558 +        }
   1.559 +        _min_cut = pmc;
   1.560 +      }
   1.561 +
   1.562 +      for (typename Graph::Node n = _first_node;
   1.563 +           n != INVALID; n = (*_nodes)[n].next) {
   1.564 +        bool merged = false;
   1.565 +        for (int a = (*_nodes)[n].first_arc; a != -1; a = _arcs[a].next) {
   1.566 +          if (!(_edges[a >> 1].cut < pmc)) {
   1.567 +            if (!merged) {
   1.568 +              for (int b = (*_nodes)[n].first_arc; b != -1; b = _arcs[b].next) {
   1.569 +                (*_nodes)[_arcs[b].target].curr_arc = b;
   1.570 +              }
   1.571 +              merged = true;
   1.572 +            }
   1.573 +            typename Graph::Node m = _arcs[a].target;
   1.574 +            int nb = 0;
   1.575 +            for (int b = (*_nodes)[m].first_arc; b != -1; b = nb) {
   1.576 +              nb = _arcs[b].next;
   1.577 +              if ((b ^ a) == 1) continue;
   1.578 +              typename Graph::Node o = _arcs[b].target;
   1.579 +              int c = (*_nodes)[o].curr_arc;
   1.580 +              if (c != -1 && _arcs[c ^ 1].target == n) {
   1.581 +                _edges[c >> 1].capacity += _edges[b >> 1].capacity;
   1.582 +                (*_nodes)[n].sum += _edges[b >> 1].capacity;
   1.583 +                if (_edges[b >> 1].cut < _edges[c >> 1].cut) {
   1.584 +                  _edges[b >> 1].cut = _edges[c >> 1].cut;
   1.585 +                }
   1.586 +                if (_arcs[b ^ 1].prev != -1) {
   1.587 +                  _arcs[_arcs[b ^ 1].prev].next = _arcs[b ^ 1].next;
   1.588 +                } else {
   1.589 +                  (*_nodes)[o].first_arc = _arcs[b ^ 1].next;
   1.590 +                }
   1.591 +                if (_arcs[b ^ 1].next != -1) {
   1.592 +                  _arcs[_arcs[b ^ 1].next].prev = _arcs[b ^ 1].prev;
   1.593 +                }
   1.594 +              } else {
   1.595 +                if (_arcs[a].next != -1) {
   1.596 +                  _arcs[_arcs[a].next].prev = b;
   1.597 +                }
   1.598 +                _arcs[b].next = _arcs[a].next;
   1.599 +                _arcs[b].prev = a;
   1.600 +                _arcs[a].next = b;
   1.601 +                _arcs[b ^ 1].target = n;
   1.602 +
   1.603 +                (*_nodes)[n].sum += _edges[b >> 1].capacity;
   1.604 +                (*_nodes)[o].curr_arc = b;
   1.605 +              }
   1.606 +            }
   1.607 +
   1.608 +            if (_arcs[a].prev != -1) {
   1.609 +              _arcs[_arcs[a].prev].next = _arcs[a].next;
   1.610 +            } else {
   1.611 +              (*_nodes)[n].first_arc = _arcs[a].next;
   1.612 +            }
   1.613 +            if (_arcs[a].next != -1) {
   1.614 +              _arcs[_arcs[a].next].prev = _arcs[a].prev;
   1.615 +            }
   1.616 +
   1.617 +            (*_nodes)[n].sum -= _edges[a >> 1].capacity;
   1.618 +            (*_next_rep)[(*_nodes)[n].last_rep] = m;
   1.619 +            (*_nodes)[n].last_rep = (*_nodes)[m].last_rep;
   1.620 +
   1.621 +            if ((*_nodes)[m].prev != INVALID) {
   1.622 +              (*_nodes)[(*_nodes)[m].prev].next = (*_nodes)[m].next;
   1.623 +            } else{
   1.624 +              _first_node = (*_nodes)[m].next;
   1.625 +            }
   1.626 +            if ((*_nodes)[m].next != INVALID) {
   1.627 +              (*_nodes)[(*_nodes)[m].next].prev = (*_nodes)[m].prev;
   1.628 +            }
   1.629 +            --_node_num;
   1.630 +          }
   1.631 +        }
   1.632 +      }
   1.633 +
   1.634 +      if (_node_num == 1) {
   1.635 +        _first_node = INVALID;
   1.636 +        return true;
   1.637 +      }
   1.638 +
   1.639 +      return false;
   1.640 +    }
   1.641 +
   1.642 +    /// \brief Executes the algorithm.
   1.643 +    ///
   1.644 +    /// Executes the algorithm.
   1.645 +    ///
   1.646 +    /// \pre init() must be called
   1.647 +    void start() {
   1.648 +      while (!processNextPhase()) {}
   1.649 +    }
   1.650 +
   1.651 +
   1.652 +    /// \brief Runs %NagamochiIbaraki algorithm.
   1.653 +    ///
   1.654 +    /// This method runs the %Min cut algorithm
   1.655 +    ///
   1.656 +    /// \note mc.run(s) is just a shortcut of the following code.
   1.657 +    ///\code
   1.658 +    ///  mc.init();
   1.659 +    ///  mc.start();
   1.660 +    ///\endcode
   1.661 +    void run() {
   1.662 +      init();
   1.663 +      start();
   1.664 +    }
   1.665 +
   1.666 +    ///@}
   1.667 +
   1.668 +    /// \name Query Functions
   1.669 +    ///
   1.670 +    /// The result of the %NagamochiIbaraki
   1.671 +    /// algorithm can be obtained using these functions.\n
   1.672 +    /// Before the use of these functions, either run() or start()
   1.673 +    /// must be called.
   1.674 +
   1.675 +    ///@{
   1.676 +
   1.677 +    /// \brief Returns the min cut value.
   1.678 +    ///
   1.679 +    /// Returns the min cut value if the algorithm finished.
   1.680 +    /// After the first processNextPhase() it is a value of a
   1.681 +    /// valid cut in the graph.
   1.682 +    Value minCutValue() const {
   1.683 +      return _min_cut;
   1.684 +    }
   1.685 +
   1.686 +    /// \brief Returns a min cut in a NodeMap.
   1.687 +    ///
   1.688 +    /// It sets the nodes of one of the two partitions to true and
   1.689 +    /// the other partition to false.
   1.690 +    /// \param cutMap A \ref concepts::WriteMap "writable" node map with
   1.691 +    /// \c bool (or convertible) value type.
   1.692 +    template <typename CutMap>
   1.693 +    Value minCutMap(CutMap& cutMap) const {
   1.694 +      for (typename Graph::NodeIt n(_graph); n != INVALID; ++n) {
   1.695 +        cutMap.set(n, (*_cut_map)[n]);
   1.696 +      }
   1.697 +      return minCutValue();
   1.698 +    }
   1.699 +
   1.700 +    ///@}
   1.701 +
   1.702 +  };
   1.703 +}
   1.704 +
   1.705 +#endif