New implementation for Nagamochi-Ibaraki algorithm
authorBalazs Dezso <deba@inf.elte.hu>
Sun, 14 Nov 2010 09:25:03 +0100
changeset 10175087694945e4
parent 1011 1937b6455b7d
child 1018 35ba7236bd67
New implementation for Nagamochi-Ibaraki algorithm
lemon/Makefile.am
lemon/nagamochi_ibaraki.h
test/CMakeLists.txt
test/Makefile.am
test/nagamochi_ibaraki_test.cc
     1.1 --- a/lemon/Makefile.am	Wed Sep 22 09:38:23 2010 +0200
     1.2 +++ b/lemon/Makefile.am	Sun Nov 14 09:25:03 2010 +0100
     1.3 @@ -107,6 +107,7 @@
     1.4  	lemon/matching.h \
     1.5  	lemon/math.h \
     1.6  	lemon/min_cost_arborescence.h \
     1.7 +	lemon/nagamochi_ibaraki.h \
     1.8  	lemon/nauty_reader.h \
     1.9  	lemon/network_simplex.h \
    1.10  	lemon/pairing_heap.h \
     2.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     2.2 +++ b/lemon/nagamochi_ibaraki.h	Sun Nov 14 09:25:03 2010 +0100
     2.3 @@ -0,0 +1,697 @@
     2.4 +/* -*- mode: C++; indent-tabs-mode: nil; -*-
     2.5 + *
     2.6 + * This file is a part of LEMON, a generic C++ optimization library.
     2.7 + *
     2.8 + * Copyright (C) 2003-2010
     2.9 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
    2.10 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
    2.11 + *
    2.12 + * Permission to use, modify and distribute this software is granted
    2.13 + * provided that this copyright notice appears in all copies. For
    2.14 + * precise terms see the accompanying LICENSE file.
    2.15 + *
    2.16 + * This software is provided "AS IS" with no warranty of any kind,
    2.17 + * express or implied, and with no claim as to its suitability for any
    2.18 + * purpose.
    2.19 + *
    2.20 + */
    2.21 +
    2.22 +#ifndef LEMON_NAGAMOCHI_IBARAKI_H
    2.23 +#define LEMON_NAGAMOCHI_IBARAKI_H
    2.24 +
    2.25 +
    2.26 +/// \ingroup min_cut
    2.27 +/// \file
    2.28 +/// \brief Implementation of the Nagamochi-Ibaraki algorithm.
    2.29 +
    2.30 +#include <lemon/core.h>
    2.31 +#include <lemon/bin_heap.h>
    2.32 +#include <lemon/bucket_heap.h>
    2.33 +#include <lemon/maps.h>
    2.34 +#include <lemon/radix_sort.h>
    2.35 +#include <lemon/unionfind.h>
    2.36 +
    2.37 +#include <cassert>
    2.38 +
    2.39 +namespace lemon {
    2.40 +
    2.41 +  /// \brief Default traits class for NagamochiIbaraki class.
    2.42 +  ///
    2.43 +  /// Default traits class for NagamochiIbaraki class.
    2.44 +  /// \param GR The undirected graph type.
    2.45 +  /// \param CM Type of capacity map.
    2.46 +  template <typename GR, typename CM>
    2.47 +  struct NagamochiIbarakiDefaultTraits {
    2.48 +    /// The type of the capacity map.
    2.49 +    typedef typename CM::Value Value;
    2.50 +
    2.51 +    /// The undirected graph type the algorithm runs on.
    2.52 +    typedef GR Graph;
    2.53 +
    2.54 +    /// \brief The type of the map that stores the edge capacities.
    2.55 +    ///
    2.56 +    /// The type of the map that stores the edge capacities.
    2.57 +    /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
    2.58 +    typedef CM CapacityMap;
    2.59 +
    2.60 +    /// \brief Instantiates a CapacityMap.
    2.61 +    ///
    2.62 +    /// This function instantiates a \ref CapacityMap.
    2.63 +#ifdef DOXYGEN
    2.64 +    static CapacityMap *createCapacityMap(const Graph& graph)
    2.65 +#else
    2.66 +    static CapacityMap *createCapacityMap(const Graph&)
    2.67 +#endif
    2.68 +    {
    2.69 +        LEMON_ASSERT(false, "CapacityMap is not initialized");
    2.70 +        return 0; // ignore warnings
    2.71 +    }
    2.72 +
    2.73 +    /// \brief The cross reference type used by heap.
    2.74 +    ///
    2.75 +    /// The cross reference type used by heap.
    2.76 +    /// Usually \c Graph::NodeMap<int>.
    2.77 +    typedef typename Graph::template NodeMap<int> HeapCrossRef;
    2.78 +
    2.79 +    /// \brief Instantiates a HeapCrossRef.
    2.80 +    ///
    2.81 +    /// This function instantiates a \ref HeapCrossRef.
    2.82 +    /// \param g is the graph, to which we would like to define the
    2.83 +    /// \ref HeapCrossRef.
    2.84 +    static HeapCrossRef *createHeapCrossRef(const Graph& g) {
    2.85 +      return new HeapCrossRef(g);
    2.86 +    }
    2.87 +
    2.88 +    /// \brief The heap type used by NagamochiIbaraki algorithm.
    2.89 +    ///
    2.90 +    /// The heap type used by NagamochiIbaraki algorithm. It has to
    2.91 +    /// maximize the priorities.
    2.92 +    ///
    2.93 +    /// \sa BinHeap
    2.94 +    /// \sa NagamochiIbaraki
    2.95 +    typedef BinHeap<Value, HeapCrossRef, std::greater<Value> > Heap;
    2.96 +
    2.97 +    /// \brief Instantiates a Heap.
    2.98 +    ///
    2.99 +    /// This function instantiates a \ref Heap.
   2.100 +    /// \param r is the cross reference of the heap.
   2.101 +    static Heap *createHeap(HeapCrossRef& r) {
   2.102 +      return new Heap(r);
   2.103 +    }
   2.104 +  };
   2.105 +
   2.106 +  /// \ingroup min_cut
   2.107 +  ///
   2.108 +  /// \brief Calculates the minimum cut in an undirected graph.
   2.109 +  ///
   2.110 +  /// Calculates the minimum cut in an undirected graph with the
   2.111 +  /// Nagamochi-Ibaraki algorithm. The algorithm separates the graph's
   2.112 +  /// nodes into two partitions with the minimum sum of edge capacities
   2.113 +  /// between the two partitions. The algorithm can be used to test
   2.114 +  /// the network reliability, especially to test how many links have
   2.115 +  /// to be destroyed in the network to split it to at least two
   2.116 +  /// distinict subnetworks.
   2.117 +  ///
   2.118 +  /// The complexity of the algorithm is \f$ O(nm\log(n)) \f$ but with
   2.119 +  /// \ref FibHeap "Fibonacci heap" it can be decreased to
   2.120 +  /// \f$ O(nm+n^2\log(n)) \f$.  When the edges have unit capacities,
   2.121 +  /// \c BucketHeap can be used which yields \f$ O(nm) \f$ time
   2.122 +  /// complexity.
   2.123 +  ///
   2.124 +  /// \warning The value type of the capacity map should be able to
   2.125 +  /// hold any cut value of the graph, otherwise the result can
   2.126 +  /// overflow.
   2.127 +  /// \note This capacity is supposed to be integer type.
   2.128 +#ifdef DOXYGEN
   2.129 +  template <typename GR, typename CM, typename TR>
   2.130 +#else
   2.131 +  template <typename GR,
   2.132 +            typename CM = typename GR::template EdgeMap<int>,
   2.133 +            typename TR = NagamochiIbarakiDefaultTraits<GR, CM> >
   2.134 +#endif
   2.135 +  class NagamochiIbaraki {
   2.136 +  public:
   2.137 +
   2.138 +    typedef TR Traits;
   2.139 +    /// The type of the underlying graph.
   2.140 +    typedef typename Traits::Graph Graph;
   2.141 +
   2.142 +    /// The type of the capacity map.
   2.143 +    typedef typename Traits::CapacityMap CapacityMap;
   2.144 +    /// The value type of the capacity map.
   2.145 +    typedef typename Traits::CapacityMap::Value Value;
   2.146 +
   2.147 +    /// The heap type used by the algorithm.
   2.148 +    typedef typename Traits::Heap Heap;
   2.149 +    /// The cross reference type used for the heap.
   2.150 +    typedef typename Traits::HeapCrossRef HeapCrossRef;
   2.151 +
   2.152 +    ///\name Named template parameters
   2.153 +
   2.154 +    ///@{
   2.155 +
   2.156 +    struct SetUnitCapacityTraits : public Traits {
   2.157 +      typedef ConstMap<typename Graph::Edge, Const<int, 1> > CapacityMap;
   2.158 +      static CapacityMap *createCapacityMap(const Graph&) {
   2.159 +        return new CapacityMap();
   2.160 +      }
   2.161 +    };
   2.162 +
   2.163 +    /// \brief \ref named-templ-param "Named parameter" for setting
   2.164 +    /// the capacity map to a constMap<Edge, int, 1>() instance
   2.165 +    ///
   2.166 +    /// \ref named-templ-param "Named parameter" for setting
   2.167 +    /// the capacity map to a constMap<Edge, int, 1>() instance
   2.168 +    struct SetUnitCapacity
   2.169 +      : public NagamochiIbaraki<Graph, CapacityMap,
   2.170 +                                SetUnitCapacityTraits> {
   2.171 +      typedef NagamochiIbaraki<Graph, CapacityMap,
   2.172 +                               SetUnitCapacityTraits> Create;
   2.173 +    };
   2.174 +
   2.175 +
   2.176 +    template <class H, class CR>
   2.177 +    struct SetHeapTraits : public Traits {
   2.178 +      typedef CR HeapCrossRef;
   2.179 +      typedef H Heap;
   2.180 +      static HeapCrossRef *createHeapCrossRef(int num) {
   2.181 +        LEMON_ASSERT(false, "HeapCrossRef is not initialized");
   2.182 +        return 0; // ignore warnings
   2.183 +      }
   2.184 +      static Heap *createHeap(HeapCrossRef &) {
   2.185 +        LEMON_ASSERT(false, "Heap is not initialized");
   2.186 +        return 0; // ignore warnings
   2.187 +      }
   2.188 +    };
   2.189 +
   2.190 +    /// \brief \ref named-templ-param "Named parameter" for setting
   2.191 +    /// heap and cross reference type
   2.192 +    ///
   2.193 +    /// \ref named-templ-param "Named parameter" for setting heap and
   2.194 +    /// cross reference type. The heap has to maximize the priorities.
   2.195 +    template <class H, class CR = RangeMap<int> >
   2.196 +    struct SetHeap
   2.197 +      : public NagamochiIbaraki<Graph, CapacityMap, SetHeapTraits<H, CR> > {
   2.198 +      typedef NagamochiIbaraki< Graph, CapacityMap, SetHeapTraits<H, CR> >
   2.199 +      Create;
   2.200 +    };
   2.201 +
   2.202 +    template <class H, class CR>
   2.203 +    struct SetStandardHeapTraits : public Traits {
   2.204 +      typedef CR HeapCrossRef;
   2.205 +      typedef H Heap;
   2.206 +      static HeapCrossRef *createHeapCrossRef(int size) {
   2.207 +        return new HeapCrossRef(size);
   2.208 +      }
   2.209 +      static Heap *createHeap(HeapCrossRef &crossref) {
   2.210 +        return new Heap(crossref);
   2.211 +      }
   2.212 +    };
   2.213 +
   2.214 +    /// \brief \ref named-templ-param "Named parameter" for setting
   2.215 +    /// heap and cross reference type with automatic allocation
   2.216 +    ///
   2.217 +    /// \ref named-templ-param "Named parameter" for setting heap and
   2.218 +    /// cross reference type with automatic allocation. They should
   2.219 +    /// have standard constructor interfaces to be able to
   2.220 +    /// automatically created by the algorithm (i.e. the graph should
   2.221 +    /// be passed to the constructor of the cross reference and the
   2.222 +    /// cross reference should be passed to the constructor of the
   2.223 +    /// heap). However, external heap and cross reference objects
   2.224 +    /// could also be passed to the algorithm using the \ref heap()
   2.225 +    /// function before calling \ref run() or \ref init(). The heap
   2.226 +    /// has to maximize the priorities.
   2.227 +    /// \sa SetHeap
   2.228 +    template <class H, class CR = RangeMap<int> >
   2.229 +    struct SetStandardHeap
   2.230 +      : public NagamochiIbaraki<Graph, CapacityMap,
   2.231 +                                SetStandardHeapTraits<H, CR> > {
   2.232 +      typedef NagamochiIbaraki<Graph, CapacityMap,
   2.233 +                               SetStandardHeapTraits<H, CR> > Create;
   2.234 +    };
   2.235 +
   2.236 +    ///@}
   2.237 +
   2.238 +
   2.239 +  private:
   2.240 +
   2.241 +    const Graph &_graph;
   2.242 +    const CapacityMap *_capacity;
   2.243 +    bool _local_capacity; // unit capacity
   2.244 +
   2.245 +    struct ArcData {
   2.246 +      typename Graph::Node target;
   2.247 +      int prev, next;
   2.248 +    };
   2.249 +    struct EdgeData {
   2.250 +      Value capacity;
   2.251 +      Value cut;
   2.252 +    };
   2.253 +
   2.254 +    struct NodeData {
   2.255 +      int first_arc;
   2.256 +      typename Graph::Node prev, next;
   2.257 +      int curr_arc;
   2.258 +      typename Graph::Node last_rep;
   2.259 +      Value sum;
   2.260 +    };
   2.261 +
   2.262 +    typename Graph::template NodeMap<NodeData> *_nodes;
   2.263 +    std::vector<ArcData> _arcs;
   2.264 +    std::vector<EdgeData> _edges;
   2.265 +
   2.266 +    typename Graph::Node _first_node;
   2.267 +    int _node_num;
   2.268 +
   2.269 +    Value _min_cut;
   2.270 +
   2.271 +    HeapCrossRef *_heap_cross_ref;
   2.272 +    bool _local_heap_cross_ref;
   2.273 +    Heap *_heap;
   2.274 +    bool _local_heap;
   2.275 +
   2.276 +    typedef typename Graph::template NodeMap<typename Graph::Node> NodeList;
   2.277 +    NodeList *_next_rep;
   2.278 +
   2.279 +    typedef typename Graph::template NodeMap<bool> MinCutMap;
   2.280 +    MinCutMap *_cut_map;
   2.281 +
   2.282 +    void createStructures() {
   2.283 +      if (!_nodes) {
   2.284 +        _nodes = new (typename Graph::template NodeMap<NodeData>)(_graph);
   2.285 +      }
   2.286 +      if (!_capacity) {
   2.287 +        _local_capacity = true;
   2.288 +        _capacity = Traits::createCapacityMap(_graph);
   2.289 +      }
   2.290 +      if (!_heap_cross_ref) {
   2.291 +        _local_heap_cross_ref = true;
   2.292 +        _heap_cross_ref = Traits::createHeapCrossRef(_graph);
   2.293 +      }
   2.294 +      if (!_heap) {
   2.295 +        _local_heap = true;
   2.296 +        _heap = Traits::createHeap(*_heap_cross_ref);
   2.297 +      }
   2.298 +      if (!_next_rep) {
   2.299 +        _next_rep = new NodeList(_graph);
   2.300 +      }
   2.301 +      if (!_cut_map) {
   2.302 +        _cut_map = new MinCutMap(_graph);
   2.303 +      }
   2.304 +    }
   2.305 +
   2.306 +  public :
   2.307 +
   2.308 +    typedef NagamochiIbaraki Create;
   2.309 +
   2.310 +
   2.311 +    /// \brief Constructor.
   2.312 +    ///
   2.313 +    /// \param graph The graph the algorithm runs on.
   2.314 +    /// \param capacity The capacity map used by the algorithm.
   2.315 +    NagamochiIbaraki(const Graph& graph, const CapacityMap& capacity)
   2.316 +      : _graph(graph), _capacity(&capacity), _local_capacity(false),
   2.317 +        _nodes(0), _arcs(), _edges(), _min_cut(),
   2.318 +        _heap_cross_ref(0), _local_heap_cross_ref(false),
   2.319 +        _heap(0), _local_heap(false),
   2.320 +        _next_rep(0), _cut_map(0) {}
   2.321 +
   2.322 +    /// \brief Constructor.
   2.323 +    ///
   2.324 +    /// This constructor can be used only when the Traits class
   2.325 +    /// defines how can the local capacity map be instantiated.
   2.326 +    /// If the SetUnitCapacity used the algorithm automatically
   2.327 +    /// constructs the capacity map.
   2.328 +    ///
   2.329 +    ///\param graph The graph the algorithm runs on.
   2.330 +    NagamochiIbaraki(const Graph& graph)
   2.331 +      : _graph(graph), _capacity(0), _local_capacity(false),
   2.332 +        _nodes(0), _arcs(), _edges(), _min_cut(),
   2.333 +        _heap_cross_ref(0), _local_heap_cross_ref(false),
   2.334 +        _heap(0), _local_heap(false),
   2.335 +        _next_rep(0), _cut_map(0) {}
   2.336 +
   2.337 +    /// \brief Destructor.
   2.338 +    ///
   2.339 +    /// Destructor.
   2.340 +    ~NagamochiIbaraki() {
   2.341 +      if (_local_capacity) delete _capacity;
   2.342 +      if (_nodes) delete _nodes;
   2.343 +      if (_local_heap) delete _heap;
   2.344 +      if (_local_heap_cross_ref) delete _heap_cross_ref;
   2.345 +      if (_next_rep) delete _next_rep;
   2.346 +      if (_cut_map) delete _cut_map;
   2.347 +    }
   2.348 +
   2.349 +    /// \brief Sets the heap and the cross reference used by algorithm.
   2.350 +    ///
   2.351 +    /// Sets the heap and the cross reference used by algorithm.
   2.352 +    /// If you don't use this function before calling \ref run(),
   2.353 +    /// it will allocate one. The destuctor deallocates this
   2.354 +    /// automatically allocated heap and cross reference, of course.
   2.355 +    /// \return <tt> (*this) </tt>
   2.356 +    NagamochiIbaraki &heap(Heap& hp, HeapCrossRef &cr)
   2.357 +    {
   2.358 +      if (_local_heap_cross_ref) {
   2.359 +        delete _heap_cross_ref;
   2.360 +        _local_heap_cross_ref = false;
   2.361 +      }
   2.362 +      _heap_cross_ref = &cr;
   2.363 +      if (_local_heap) {
   2.364 +        delete _heap;
   2.365 +        _local_heap = false;
   2.366 +      }
   2.367 +      _heap = &hp;
   2.368 +      return *this;
   2.369 +    }
   2.370 +
   2.371 +    /// \name Execution control
   2.372 +    /// The simplest way to execute the algorithm is to use
   2.373 +    /// one of the member functions called \c run().
   2.374 +    /// \n
   2.375 +    /// If you need more control on the execution,
   2.376 +    /// first you must call \ref init() and then call the start()
   2.377 +    /// or proper times the processNextPhase() member functions.
   2.378 +
   2.379 +    ///@{
   2.380 +
   2.381 +    /// \brief Initializes the internal data structures.
   2.382 +    ///
   2.383 +    /// Initializes the internal data structures.
   2.384 +    void init() {
   2.385 +      createStructures();
   2.386 +
   2.387 +      int edge_num = countEdges(_graph);
   2.388 +      _edges.resize(edge_num);
   2.389 +      _arcs.resize(2 * edge_num);
   2.390 +
   2.391 +      typename Graph::Node prev = INVALID;
   2.392 +      _node_num = 0;
   2.393 +      for (typename Graph::NodeIt n(_graph); n != INVALID; ++n) {
   2.394 +        (*_cut_map)[n] = false;
   2.395 +        (*_next_rep)[n] = INVALID;
   2.396 +        (*_nodes)[n].last_rep = n;
   2.397 +        (*_nodes)[n].first_arc = -1;
   2.398 +        (*_nodes)[n].curr_arc = -1;
   2.399 +        (*_nodes)[n].prev = prev;
   2.400 +        if (prev != INVALID) {
   2.401 +          (*_nodes)[prev].next = n;
   2.402 +        }
   2.403 +        (*_nodes)[n].next = INVALID;
   2.404 +        (*_nodes)[n].sum = 0;
   2.405 +        prev = n;
   2.406 +        ++_node_num;
   2.407 +      }
   2.408 +
   2.409 +      _first_node = typename Graph::NodeIt(_graph);
   2.410 +
   2.411 +      int index = 0;
   2.412 +      for (typename Graph::NodeIt n(_graph); n != INVALID; ++n) {
   2.413 +        for (typename Graph::OutArcIt a(_graph, n); a != INVALID; ++a) {
   2.414 +          typename Graph::Node m = _graph.target(a);
   2.415 +          
   2.416 +          if (!(n < m)) continue;
   2.417 +
   2.418 +          (*_nodes)[n].sum += (*_capacity)[a];
   2.419 +          (*_nodes)[m].sum += (*_capacity)[a];
   2.420 +          
   2.421 +          int c = (*_nodes)[m].curr_arc;
   2.422 +          if (c != -1 && _arcs[c ^ 1].target == n) {
   2.423 +            _edges[c >> 1].capacity += (*_capacity)[a];
   2.424 +          } else {
   2.425 +            _edges[index].capacity = (*_capacity)[a];
   2.426 +            
   2.427 +            _arcs[index << 1].prev = -1;
   2.428 +            if ((*_nodes)[n].first_arc != -1) {
   2.429 +              _arcs[(*_nodes)[n].first_arc].prev = (index << 1);
   2.430 +            }
   2.431 +            _arcs[index << 1].next = (*_nodes)[n].first_arc;
   2.432 +            (*_nodes)[n].first_arc = (index << 1);
   2.433 +            _arcs[index << 1].target = m;
   2.434 +
   2.435 +            (*_nodes)[m].curr_arc = (index << 1);
   2.436 +            
   2.437 +            _arcs[(index << 1) | 1].prev = -1;
   2.438 +            if ((*_nodes)[m].first_arc != -1) {
   2.439 +              _arcs[(*_nodes)[m].first_arc].prev = ((index << 1) | 1);
   2.440 +            }
   2.441 +            _arcs[(index << 1) | 1].next = (*_nodes)[m].first_arc;
   2.442 +            (*_nodes)[m].first_arc = ((index << 1) | 1);
   2.443 +            _arcs[(index << 1) | 1].target = n;
   2.444 +            
   2.445 +            ++index;
   2.446 +          }
   2.447 +        }
   2.448 +      }
   2.449 +
   2.450 +      typename Graph::Node cut_node = INVALID;
   2.451 +      _min_cut = std::numeric_limits<Value>::max();
   2.452 +
   2.453 +      for (typename Graph::Node n = _first_node; 
   2.454 +           n != INVALID; n = (*_nodes)[n].next) {
   2.455 +        if ((*_nodes)[n].sum < _min_cut) {
   2.456 +          cut_node = n;
   2.457 +          _min_cut = (*_nodes)[n].sum;
   2.458 +        }
   2.459 +      }
   2.460 +      (*_cut_map)[cut_node] = true;
   2.461 +      if (_min_cut == 0) {
   2.462 +        _first_node = INVALID;
   2.463 +      }
   2.464 +    }
   2.465 +
   2.466 +  public:
   2.467 +
   2.468 +    /// \brief Processes the next phase
   2.469 +    ///
   2.470 +    /// Processes the next phase in the algorithm. It must be called
   2.471 +    /// at most one less the number of the nodes in the graph.
   2.472 +    ///
   2.473 +    ///\return %True when the algorithm finished.
   2.474 +    bool processNextPhase() {
   2.475 +      if (_first_node == INVALID) return true;
   2.476 +
   2.477 +      _heap->clear();
   2.478 +      for (typename Graph::Node n = _first_node; 
   2.479 +           n != INVALID; n = (*_nodes)[n].next) {
   2.480 +        (*_heap_cross_ref)[n] = Heap::PRE_HEAP;
   2.481 +      }
   2.482 +
   2.483 +      std::vector<typename Graph::Node> order;
   2.484 +      order.reserve(_node_num);
   2.485 +      int sep = 0;
   2.486 +
   2.487 +      Value alpha = 0;
   2.488 +      Value pmc = std::numeric_limits<Value>::max();
   2.489 +
   2.490 +      _heap->push(_first_node, static_cast<Value>(0));
   2.491 +      while (!_heap->empty()) {
   2.492 +        typename Graph::Node n = _heap->top();
   2.493 +        Value v = _heap->prio();
   2.494 +
   2.495 +        _heap->pop();
   2.496 +        for (int a = (*_nodes)[n].first_arc; a != -1; a = _arcs[a].next) {
   2.497 +          switch (_heap->state(_arcs[a].target)) {
   2.498 +          case Heap::PRE_HEAP: 
   2.499 +            {
   2.500 +              Value nv = _edges[a >> 1].capacity;
   2.501 +              _heap->push(_arcs[a].target, nv);
   2.502 +              _edges[a >> 1].cut = nv;
   2.503 +            } break;
   2.504 +          case Heap::IN_HEAP:
   2.505 +            {
   2.506 +              Value nv = _edges[a >> 1].capacity + (*_heap)[_arcs[a].target];
   2.507 +              _heap->decrease(_arcs[a].target, nv);
   2.508 +              _edges[a >> 1].cut = nv;
   2.509 +            } break;
   2.510 +          case Heap::POST_HEAP:
   2.511 +            break;
   2.512 +          }
   2.513 +        }
   2.514 +
   2.515 +        alpha += (*_nodes)[n].sum;
   2.516 +        alpha -= 2 * v;
   2.517 +
   2.518 +        order.push_back(n);
   2.519 +        if (!_heap->empty()) {
   2.520 +          if (alpha < pmc) {
   2.521 +            pmc = alpha;
   2.522 +            sep = order.size();
   2.523 +          }
   2.524 +        }
   2.525 +      }
   2.526 +
   2.527 +      if (static_cast<int>(order.size()) < _node_num) {
   2.528 +        _first_node = INVALID;
   2.529 +        for (typename Graph::NodeIt n(_graph); n != INVALID; ++n) {
   2.530 +          (*_cut_map)[n] = false;
   2.531 +        }
   2.532 +        for (int i = 0; i < static_cast<int>(order.size()); ++i) {
   2.533 +          typename Graph::Node n = order[i];
   2.534 +          while (n != INVALID) {
   2.535 +            (*_cut_map)[n] = true;
   2.536 +            n = (*_next_rep)[n];
   2.537 +          }
   2.538 +        }
   2.539 +        _min_cut = 0;
   2.540 +        return true;
   2.541 +      }
   2.542 +
   2.543 +      if (pmc < _min_cut) {
   2.544 +        for (typename Graph::NodeIt n(_graph); n != INVALID; ++n) {
   2.545 +          (*_cut_map)[n] = false;
   2.546 +        }
   2.547 +        for (int i = 0; i < sep; ++i) {
   2.548 +          typename Graph::Node n = order[i];
   2.549 +          while (n != INVALID) {
   2.550 +            (*_cut_map)[n] = true;
   2.551 +            n = (*_next_rep)[n];
   2.552 +          }
   2.553 +        }
   2.554 +        _min_cut = pmc;
   2.555 +      }
   2.556 +
   2.557 +      for (typename Graph::Node n = _first_node;
   2.558 +           n != INVALID; n = (*_nodes)[n].next) {
   2.559 +        bool merged = false;
   2.560 +        for (int a = (*_nodes)[n].first_arc; a != -1; a = _arcs[a].next) {
   2.561 +          if (!(_edges[a >> 1].cut < pmc)) {
   2.562 +            if (!merged) {
   2.563 +              for (int b = (*_nodes)[n].first_arc; b != -1; b = _arcs[b].next) {
   2.564 +                (*_nodes)[_arcs[b].target].curr_arc = b;          
   2.565 +              }
   2.566 +              merged = true;
   2.567 +            }
   2.568 +            typename Graph::Node m = _arcs[a].target;
   2.569 +            int nb = 0;
   2.570 +            for (int b = (*_nodes)[m].first_arc; b != -1; b = nb) {
   2.571 +              nb = _arcs[b].next;
   2.572 +              if ((b ^ a) == 1) continue;
   2.573 +              typename Graph::Node o = _arcs[b].target;
   2.574 +              int c = (*_nodes)[o].curr_arc; 
   2.575 +              if (c != -1 && _arcs[c ^ 1].target == n) {
   2.576 +                _edges[c >> 1].capacity += _edges[b >> 1].capacity;
   2.577 +                (*_nodes)[n].sum += _edges[b >> 1].capacity;
   2.578 +                if (_edges[b >> 1].cut < _edges[c >> 1].cut) {
   2.579 +                  _edges[b >> 1].cut = _edges[c >> 1].cut;
   2.580 +                }
   2.581 +                if (_arcs[b ^ 1].prev != -1) {
   2.582 +                  _arcs[_arcs[b ^ 1].prev].next = _arcs[b ^ 1].next;
   2.583 +                } else {
   2.584 +                  (*_nodes)[o].first_arc = _arcs[b ^ 1].next;
   2.585 +                }
   2.586 +                if (_arcs[b ^ 1].next != -1) {
   2.587 +                  _arcs[_arcs[b ^ 1].next].prev = _arcs[b ^ 1].prev;
   2.588 +                }
   2.589 +              } else {
   2.590 +                if (_arcs[a].next != -1) {
   2.591 +                  _arcs[_arcs[a].next].prev = b;
   2.592 +                }
   2.593 +                _arcs[b].next = _arcs[a].next;
   2.594 +                _arcs[b].prev = a;
   2.595 +                _arcs[a].next = b;
   2.596 +                _arcs[b ^ 1].target = n;
   2.597 +
   2.598 +                (*_nodes)[n].sum += _edges[b >> 1].capacity;
   2.599 +                (*_nodes)[o].curr_arc = b;
   2.600 +              }
   2.601 +            }
   2.602 +
   2.603 +            if (_arcs[a].prev != -1) {
   2.604 +              _arcs[_arcs[a].prev].next = _arcs[a].next;
   2.605 +            } else {
   2.606 +              (*_nodes)[n].first_arc = _arcs[a].next;
   2.607 +            }            
   2.608 +            if (_arcs[a].next != -1) {
   2.609 +              _arcs[_arcs[a].next].prev = _arcs[a].prev;
   2.610 +            }
   2.611 +
   2.612 +            (*_nodes)[n].sum -= _edges[a >> 1].capacity;
   2.613 +            (*_next_rep)[(*_nodes)[n].last_rep] = m;
   2.614 +            (*_nodes)[n].last_rep = (*_nodes)[m].last_rep;
   2.615 +            
   2.616 +            if ((*_nodes)[m].prev != INVALID) {
   2.617 +              (*_nodes)[(*_nodes)[m].prev].next = (*_nodes)[m].next;
   2.618 +            } else{
   2.619 +              _first_node = (*_nodes)[m].next;
   2.620 +            }
   2.621 +            if ((*_nodes)[m].next != INVALID) {
   2.622 +              (*_nodes)[(*_nodes)[m].next].prev = (*_nodes)[m].prev;
   2.623 +            }
   2.624 +            --_node_num;
   2.625 +          }
   2.626 +        }
   2.627 +      }
   2.628 +
   2.629 +      if (_node_num == 1) {
   2.630 +        _first_node = INVALID;
   2.631 +        return true;
   2.632 +      }
   2.633 +
   2.634 +      return false;
   2.635 +    }
   2.636 +
   2.637 +    /// \brief Executes the algorithm.
   2.638 +    ///
   2.639 +    /// Executes the algorithm.
   2.640 +    ///
   2.641 +    /// \pre init() must be called
   2.642 +    void start() {
   2.643 +      while (!processNextPhase()) {}
   2.644 +    }
   2.645 +
   2.646 +
   2.647 +    /// \brief Runs %NagamochiIbaraki algorithm.
   2.648 +    ///
   2.649 +    /// This method runs the %Min cut algorithm
   2.650 +    ///
   2.651 +    /// \note mc.run(s) is just a shortcut of the following code.
   2.652 +    ///\code
   2.653 +    ///  mc.init();
   2.654 +    ///  mc.start();
   2.655 +    ///\endcode
   2.656 +    void run() {
   2.657 +      init();
   2.658 +      start();
   2.659 +    }
   2.660 +
   2.661 +    ///@}
   2.662 +
   2.663 +    /// \name Query Functions
   2.664 +    ///
   2.665 +    /// The result of the %NagamochiIbaraki
   2.666 +    /// algorithm can be obtained using these functions.\n
   2.667 +    /// Before the use of these functions, either run() or start()
   2.668 +    /// must be called.
   2.669 +
   2.670 +    ///@{
   2.671 +
   2.672 +    /// \brief Returns the min cut value.
   2.673 +    ///
   2.674 +    /// Returns the min cut value if the algorithm finished.
   2.675 +    /// After the first processNextPhase() it is a value of a
   2.676 +    /// valid cut in the graph.
   2.677 +    Value minCutValue() const {
   2.678 +      return _min_cut;
   2.679 +    }
   2.680 +
   2.681 +    /// \brief Returns a min cut in a NodeMap.
   2.682 +    ///
   2.683 +    /// It sets the nodes of one of the two partitions to true and
   2.684 +    /// the other partition to false.
   2.685 +    /// \param cutMap A \ref concepts::WriteMap "writable" node map with
   2.686 +    /// \c bool (or convertible) value type.
   2.687 +    template <typename CutMap>
   2.688 +    Value minCutMap(CutMap& cutMap) const {
   2.689 +      for (typename Graph::NodeIt n(_graph); n != INVALID; ++n) {
   2.690 +        cutMap.set(n, (*_cut_map)[n]);
   2.691 +      }
   2.692 +      return minCutValue();
   2.693 +    }
   2.694 +
   2.695 +    ///@}
   2.696 +
   2.697 +  };
   2.698 +}
   2.699 +
   2.700 +#endif
     3.1 --- a/test/CMakeLists.txt	Wed Sep 22 09:38:23 2010 +0200
     3.2 +++ b/test/CMakeLists.txt	Sun Nov 14 09:25:03 2010 +0100
     3.3 @@ -35,6 +35,7 @@
     3.4    min_cost_arborescence_test
     3.5    min_cost_flow_test
     3.6    min_mean_cycle_test
     3.7 +  nagamochi_ibaraki_test
     3.8    path_test
     3.9    planarity_test
    3.10    preflow_test
     4.1 --- a/test/Makefile.am	Wed Sep 22 09:38:23 2010 +0200
     4.2 +++ b/test/Makefile.am	Sun Nov 14 09:25:03 2010 +0100
     4.3 @@ -37,6 +37,7 @@
     4.4  	test/min_cost_arborescence_test \
     4.5  	test/min_cost_flow_test \
     4.6  	test/min_mean_cycle_test \
     4.7 +	test/nagamochi_ibaraki_test \
     4.8  	test/path_test \
     4.9  	test/planarity_test \
    4.10  	test/preflow_test \
    4.11 @@ -89,6 +90,7 @@
    4.12  test_min_cost_arborescence_test_SOURCES = test/min_cost_arborescence_test.cc
    4.13  test_min_cost_flow_test_SOURCES = test/min_cost_flow_test.cc
    4.14  test_min_mean_cycle_test_SOURCES = test/min_mean_cycle_test.cc
    4.15 +test_nagamochi_ibaraki_test_SOURCES = test/nagamochi_ibaraki_test.cc
    4.16  test_path_test_SOURCES = test/path_test.cc
    4.17  test_planarity_test_SOURCES = test/planarity_test.cc
    4.18  test_preflow_test_SOURCES = test/preflow_test.cc
     5.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     5.2 +++ b/test/nagamochi_ibaraki_test.cc	Sun Nov 14 09:25:03 2010 +0100
     5.3 @@ -0,0 +1,141 @@
     5.4 +/* -*- mode: C++; indent-tabs-mode: nil; -*-
     5.5 + *
     5.6 + * This file is a part of LEMON, a generic C++ optimization library.
     5.7 + *
     5.8 + * Copyright (C) 2003-2010
     5.9 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
    5.10 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
    5.11 + *
    5.12 + * Permission to use, modify and distribute this software is granted
    5.13 + * provided that this copyright notice appears in all copies. For
    5.14 + * precise terms see the accompanying LICENSE file.
    5.15 + *
    5.16 + * This software is provided "AS IS" with no warranty of any kind,
    5.17 + * express or implied, and with no claim as to its suitability for any
    5.18 + * purpose.
    5.19 + *
    5.20 + */
    5.21 +
    5.22 +#include <sstream>
    5.23 +
    5.24 +#include <lemon/smart_graph.h>
    5.25 +#include <lemon/adaptors.h>
    5.26 +#include <lemon/concepts/graph.h>
    5.27 +#include <lemon/concepts/maps.h>
    5.28 +#include <lemon/lgf_reader.h>
    5.29 +#include <lemon/nagamochi_ibaraki.h>
    5.30 +
    5.31 +#include "test_tools.h"
    5.32 +
    5.33 +using namespace lemon;
    5.34 +using namespace std;
    5.35 +
    5.36 +const std::string lgf =
    5.37 +  "@nodes\n"
    5.38 +  "label\n"
    5.39 +  "0\n"
    5.40 +  "1\n"
    5.41 +  "2\n"
    5.42 +  "3\n"
    5.43 +  "4\n"
    5.44 +  "5\n"
    5.45 +  "@edges\n"
    5.46 +  "     cap1 cap2 cap3\n"
    5.47 +  "0 1  1    1    1   \n"
    5.48 +  "0 2  2    2    4   \n"
    5.49 +  "1 2  4    4    4   \n"
    5.50 +  "3 4  1    1    1   \n"
    5.51 +  "3 5  2    2    4   \n"
    5.52 +  "4 5  4    4    4   \n"
    5.53 +  "2 3  1    6    6   \n";
    5.54 +
    5.55 +void checkNagamochiIbarakiCompile()
    5.56 +{
    5.57 +  typedef int Value;
    5.58 +  typedef concepts::Graph Graph;
    5.59 +
    5.60 +  typedef Graph::Node Node;
    5.61 +  typedef Graph::Edge Edge;
    5.62 +  typedef concepts::ReadMap<Edge, Value> CapMap;
    5.63 +  typedef concepts::WriteMap<Node, bool> CutMap;
    5.64 +
    5.65 +  Graph g;
    5.66 +  Node n;
    5.67 +  CapMap cap;
    5.68 +  CutMap cut;
    5.69 +  Value v;
    5.70 +  bool b;
    5.71 +
    5.72 +  NagamochiIbaraki<Graph, CapMap> ni_test(g, cap);
    5.73 +  const NagamochiIbaraki<Graph, CapMap>& const_ni_test = ni_test;
    5.74 +
    5.75 +  ni_test.init();
    5.76 +  ni_test.start();
    5.77 +  b = ni_test.processNextPhase();
    5.78 +  ni_test.run();
    5.79 +
    5.80 +  v = const_ni_test.minCutValue();
    5.81 +  v = const_ni_test.minCutMap(cut);
    5.82 +}
    5.83 +
    5.84 +template <typename Graph, typename CapMap, typename CutMap>
    5.85 +typename CapMap::Value
    5.86 +  cutValue(const Graph& graph, const CapMap& cap, const CutMap& cut)
    5.87 +{
    5.88 +  typename CapMap::Value sum = 0;
    5.89 +  for (typename Graph::EdgeIt e(graph); e != INVALID; ++e) {
    5.90 +    if (cut[graph.u(e)] != cut[graph.v(e)]) {
    5.91 +      sum += cap[e];
    5.92 +    }
    5.93 +  }
    5.94 +  return sum;
    5.95 +}
    5.96 +
    5.97 +int main() {
    5.98 +  SmartGraph graph;
    5.99 +  SmartGraph::EdgeMap<int> cap1(graph), cap2(graph), cap3(graph);
   5.100 +  SmartGraph::NodeMap<bool> cut(graph);
   5.101 +
   5.102 +  istringstream input(lgf);
   5.103 +  graphReader(graph, input)
   5.104 +    .edgeMap("cap1", cap1)
   5.105 +    .edgeMap("cap2", cap2)
   5.106 +    .edgeMap("cap3", cap3)
   5.107 +    .run();
   5.108 +
   5.109 +  {
   5.110 +    NagamochiIbaraki<SmartGraph> ni(graph, cap1);
   5.111 +    ni.run();
   5.112 +    ni.minCutMap(cut);
   5.113 +
   5.114 +    check(ni.minCutValue() == 1, "Wrong cut value");
   5.115 +    check(ni.minCutValue() == cutValue(graph, cap1, cut), "Wrong cut value");
   5.116 +  }
   5.117 +  {
   5.118 +    NagamochiIbaraki<SmartGraph> ni(graph, cap2);
   5.119 +    ni.run();
   5.120 +    ni.minCutMap(cut);
   5.121 +
   5.122 +    check(ni.minCutValue() == 3, "Wrong cut value");
   5.123 +    check(ni.minCutValue() == cutValue(graph, cap2, cut), "Wrong cut value");
   5.124 +  }
   5.125 +  {
   5.126 +    NagamochiIbaraki<SmartGraph> ni(graph, cap3);
   5.127 +    ni.run();
   5.128 +    ni.minCutMap(cut);
   5.129 +
   5.130 +    check(ni.minCutValue() == 5, "Wrong cut value");
   5.131 +    check(ni.minCutValue() == cutValue(graph, cap3, cut), "Wrong cut value");
   5.132 +  }
   5.133 +  {
   5.134 +    NagamochiIbaraki<SmartGraph>::SetUnitCapacity::Create ni(graph);
   5.135 +    ni.run();
   5.136 +    ni.minCutMap(cut);
   5.137 +
   5.138 +    ConstMap<SmartGraph::Edge, int> cap4(1);
   5.139 +    check(ni.minCutValue() == 1, "Wrong cut value");
   5.140 +    check(ni.minCutValue() == cutValue(graph, cap4, cut), "Wrong cut value");
   5.141 +  }
   5.142 +
   5.143 +  return 0;
   5.144 +}