# HG changeset patch
# User deba
# Date 1168549500 0
# Node ID 9c3d44ac39fbc30ad9002c0a3a247823ce32c991
# Parent  215a6f3e33c906fca146ecf88c42d678f8c48102
Adding two heuristics
Based on:
http://www.avglab.com/andrew/pub/neci-tr-96-132.ps

diff -r 215a6f3e33c9 -r 9c3d44ac39fb lemon/nagamochi_ibaraki.h
--- a/lemon/nagamochi_ibaraki.h	Tue Jan 09 11:42:43 2007 +0000
+++ b/lemon/nagamochi_ibaraki.h	Thu Jan 11 21:05:00 2007 +0000
@@ -14,24 +14,31 @@
  *
  */
 
-#ifndef LEMON_MIN_CUT_H
-#define LEMON_MIN_CUT_H
+#ifndef LEMON_NAGAMOCHI_IBARAKI_H
+#define LEMON_NAGAMOCHI_IBARAKI_H
 
 
-/// \ingroup topology
-/// \file
-/// \brief Maximum cardinality search and min cut in undirected graphs.
+/// \ingroup topology 
+/// \file 
+/// \brief Maximum cardinality search and minimum cut in undirected
+/// graphs.
 
 #include <lemon/list_graph.h>
 #include <lemon/bin_heap.h>
 #include <lemon/bucket_heap.h>
 
+#include <lemon/unionfind.h>
+#include <lemon/topology.h>
+
 #include <lemon/bits/invalid.h>
 #include <lemon/error.h>
 #include <lemon/maps.h>
 
 #include <functional>
 
+#include <lemon/graph_writer.h>
+#include <lemon/time_measure.h>
+
 namespace lemon {
 
   namespace _min_cut_bits {
@@ -146,6 +153,7 @@
       return new CardinalityMap(graph);
     }
 
+
   };
   
   /// \ingroup topology
@@ -665,6 +673,12 @@
     /// of this funcion is undefined.
     Value cardinality(Node node) const { return (*_cardinality)[node]; }
 
+    /// \brief The current cardinality of a node.
+    ///
+    /// Returns the current cardinality of a node.
+    /// \pre the given node should be reached but not processed
+    Value currentCardinality(Node node) const { return (*_heap)[node]; }
+
     /// \brief Returns a reference to the NodeMap of cardinalities.
     ///
     /// Returns a reference to the NodeMap of cardinalities. \pre \ref run() 
@@ -729,9 +743,21 @@
       throw UninitializedParameter();
     }
 
+    /// \brief The CutValueMap type
+    ///
+    /// The type of the map that stores the cut value of a node.
+    typedef AuxGraph::NodeMap<Value> AuxCutValueMap;
+
+    /// \brief Instantiates a AuxCutValueMap.
+    ///
+    /// This function instantiates a \ref AuxCutValueMap. 
+    static AuxCutValueMap *createAuxCutValueMap(const AuxGraph& graph) {
+      return new AuxCutValueMap(graph);
+    }
+
     /// \brief The AuxCapacityMap type
     ///
-    /// The type of the map that stores the auxing edge capacities.
+    /// The type of the map that stores the auxiliary edge capacities.
     typedef AuxGraph::UEdgeMap<Value> AuxCapacityMap;
 
     /// \brief Instantiates a AuxCapacityMap.
@@ -806,28 +832,6 @@
 
   };
 
-  namespace _min_cut_bits {
-    template <typename _Key>
-    class LastTwoMap {
-    public:
-      typedef _Key Key;
-      typedef bool Value;
-
-      LastTwoMap(int _num) : num(_num) {}
-      void set(const Key& key, bool val) {
-        if (!val) return;
-        --num;
-        if (num > 1) return;
-        keys[num] = key;
-      }
-      
-      Key operator[](int index) const { return keys[index]; }
-    private:
-      Key keys[2];
-      int num;
-    };
-  }
-
   /// \ingroup topology
   ///
   /// \brief Calculates the minimum cut in an undirected graph.
@@ -841,9 +845,12 @@
   /// distinict subnetwork.
   ///
   /// The complexity of the algorithm is \f$ O(ne\log(n)) \f$ but with
-  /// Fibonacci heap it can be decreased to \f$ O(ne+n^2\log(n))
-  /// \f$. When capacity map is neutral then it uses BucketHeap which
+  /// Fibonacci heap it can be decreased to \f$ O(ne+n^2\log(n)) \f$. 
+  /// When capacity map is neutral then it uses BucketHeap which
   /// results \f$ O(ne) \f$ time complexity.
+  ///
+  /// \warning The value type of the capacity map should be able to hold
+  /// any cut value of the graph, otherwise the result can overflow.
 #ifdef DOXYGEN
   template <typename _Graph, typename _CapacityMap, typename _Traits>
 #else
@@ -880,6 +887,8 @@
     typedef typename Traits::AuxGraph AuxGraph;
     /// The type of the aux capacity map
     typedef typename Traits::AuxCapacityMap AuxCapacityMap;
+    /// The type of the aux cut value map
+    typedef typename Traits::AuxCutValueMap AuxCutValueMap;
     /// The cross reference type used for the current heap.
     typedef typename Traits::HeapCrossRef HeapCrossRef;
     /// The heap type used by the max cardinality algorithm.
@@ -988,6 +997,11 @@
     /// \brief Indicates if \ref _aux_capacity is locally allocated 
     /// (\c true) or not.
     bool local_aux_capacity;
+    /// Pointer to the aux cut value map
+    AuxCutValueMap *_aux_cut_value;
+    /// \brief Indicates if \ref _aux_cut_value is locally allocated 
+    /// (\c true) or not.
+    bool local_aux_cut_value;
     /// Pointer to the heap cross references.
     HeapCrossRef *_heap_cross_ref;
     /// \brief Indicates if \ref _heap_cross_ref is locally allocated 
@@ -1002,8 +1016,8 @@
     Value _min_cut;
     /// The number of the nodes of the aux graph.
     int _node_num;
-    /// The first and last node of the min cut in the next list;
-    typename Graph::Node _first_node, _last_node;
+    /// The first and last node of the min cut in the next list.
+    std::vector<typename Graph::Node> _cut;
 
     /// \brief The first and last element in the list associated
     /// to the aux graph node.
@@ -1011,7 +1025,7 @@
     /// \brief The next node in the node lists.
     ListRefMap *_next;
     
-    void create_structures() {
+    void createStructures() {
       if (!_capacity) {
         local_capacity = true;
         _capacity = Traits::createCapacityMap(*_graph);
@@ -1024,31 +1038,16 @@
 	local_aux_capacity = true;
 	_aux_capacity = Traits::createAuxCapacityMap(*_aux_graph);
       }
+      if(!_aux_cut_value) {
+	local_aux_cut_value = true;
+	_aux_cut_value = Traits::createAuxCutValueMap(*_aux_graph);
+      }
 
       _first = Traits::createNodeRefMap(*_aux_graph);
       _last = Traits::createNodeRefMap(*_aux_graph);
 
       _next = Traits::createListRefMap(*_graph);
 
-      typename Graph::template NodeMap<typename AuxGraph::Node> ref(*_graph);
-
-      for (typename Graph::NodeIt it(*_graph); it != INVALID; ++it) {
-        _next->set(it, INVALID);
-        typename AuxGraph::Node node = _aux_graph->addNode();
-        _first->set(node, it);
-        _last->set(node, it);
-        ref.set(it, node);
-      }
-
-      for (typename Graph::UEdgeIt it(*_graph); it != INVALID; ++it) {
-        if (_graph->source(it) == _graph->target(it)) continue;
-        typename AuxGraph::UEdge uedge = 
-          _aux_graph->addEdge(ref[_graph->source(it)], 
-                               ref[_graph->target(it)]);
-        _aux_capacity->set(uedge, (*_capacity)[it]);
-        
-      }
-
       if (!_heap_cross_ref) {
 	local_heap_cross_ref = true;
 	_heap_cross_ref = Traits::createHeapCrossRef(*_aux_graph);
@@ -1059,6 +1058,57 @@
       }
     }
 
+    void createAuxGraph() {
+      typename Graph::template NodeMap<typename AuxGraph::Node> ref(*_graph);
+
+      for (typename Graph::NodeIt n(*_graph); n != INVALID; ++n) {
+        _next->set(n, INVALID);
+        typename AuxGraph::Node node = _aux_graph->addNode();
+        _first->set(node, n);
+        _last->set(node, n);
+        ref.set(n, node);
+      }
+
+      typename AuxGraph::template NodeMap<typename AuxGraph::UEdge> 
+      edges(*_aux_graph, INVALID);
+
+      for (typename Graph::NodeIt n(*_graph); n != INVALID; ++n) {
+        for (typename Graph::IncEdgeIt e(*_graph, n); e != INVALID; ++e) {
+          typename Graph::Node tn = _graph->runningNode(e);
+          if (n < tn || n == tn) continue;
+          if (edges[ref[tn]] != INVALID) {
+            Value value = 
+              (*_aux_capacity)[edges[ref[tn]]] + (*_capacity)[e];
+            _aux_capacity->set(edges[ref[tn]], value);
+          } else {
+            edges.set(ref[tn], _aux_graph->addEdge(ref[n], ref[tn]));
+            Value value = (*_capacity)[e];
+            _aux_capacity->set(edges[ref[tn]], value);            
+          }
+        }
+        for (typename Graph::IncEdgeIt e(*_graph, n); e != INVALID; ++e) {
+          typename Graph::Node tn = _graph->runningNode(e);
+          edges.set(ref[tn], INVALID);
+        }
+      }
+
+      _cut.resize(1, INVALID);
+      _min_cut = std::numeric_limits<Value>::max();
+      for (typename AuxGraph::NodeIt n(*_aux_graph); n != INVALID; ++n) {
+        Value value = 0;
+        for (typename AuxGraph::IncEdgeIt e(*_aux_graph, n); 
+             e != INVALID; ++e) {
+          value += (*_aux_capacity)[e];
+        }
+        if (_min_cut > value) {
+          _min_cut = value;
+          _cut[0] = (*_first)[n];
+        } 
+        (*_aux_cut_value)[n] = value;
+      }
+    }
+    
+
   public :
 
     typedef NagamochiIbaraki Create;
@@ -1073,6 +1123,7 @@
         _capacity(&capacity), local_capacity(false),
         _aux_graph(0), local_aux_graph(false),
         _aux_capacity(0), local_aux_capacity(false),
+        _aux_cut_value(0), local_aux_cut_value(false),
         _heap_cross_ref(0), local_heap_cross_ref(false),
         _heap(0), local_heap(false),
         _first(0), _last(0), _next(0) {}
@@ -1090,6 +1141,7 @@
         _capacity(0), local_capacity(false),
         _aux_graph(0), local_aux_graph(false),
         _aux_capacity(0), local_aux_capacity(false),
+        _aux_cut_value(0), local_aux_cut_value(false),
         _heap_cross_ref(0), local_heap_cross_ref(false),
         _heap(0), local_heap(false),
         _first(0), _last(0), _next(0) {}
@@ -1104,6 +1156,7 @@
       if (_last) delete _last;
       if (_next) delete _next;
       if (local_aux_capacity) delete _aux_capacity;
+      if (local_aux_cut_value) delete _aux_cut_value;
       if (local_aux_graph) delete _aux_graph;
       if (local_capacity) delete _capacity;
     }
@@ -1178,91 +1231,221 @@
     ///
     /// Initializes the internal data structures.
     void init() {
-      create_structures();
-      _first_node = _last_node = INVALID;
       _node_num = countNodes(*_graph);
+      createStructures();
+      createAuxGraph();
     }
 
+  private:
+
+    struct EdgeInfo {
+      typename AuxGraph::Node source, target;
+      Value capacity;
+    };
+    
+    struct EdgeInfoLess {
+      bool operator()(const EdgeInfo& left, const EdgeInfo& right) const {
+        return (left.source < right.source) || 
+          (left.source == right.source && left.target < right.target);
+      }
+    };
+
+  public:
+
+
     /// \brief Processes the next phase
     ///
-    /// Processes the next phase in the algorithm. The function
-    /// should be called countNodes(graph) - 1 times to get
-    /// surely the min cut in the graph. The  
+    /// Processes the next phase in the algorithm. The function should
+    /// be called at most countNodes(graph) - 1 times to get surely
+    /// the min cut in the graph.
     ///
     ///\return %True when the algorithm finished.
     bool processNextPhase() {
       if (_node_num <= 1) return true;
-      using namespace _min_cut_bits;
 
       typedef typename AuxGraph::Node Node;
       typedef typename AuxGraph::NodeIt NodeIt;
       typedef typename AuxGraph::UEdge UEdge;
+      typedef typename AuxGraph::UEdgeIt UEdgeIt;
       typedef typename AuxGraph::IncEdgeIt IncEdgeIt;
       
-      typedef typename MaxCardinalitySearch<AuxGraph, AuxCapacityMap>::
-      template DefHeap<Heap, HeapCrossRef>::
-      template DefCardinalityMap<NullMap<Node, Value> >::
-      template DefProcessedMap<LastTwoMap<Node> >::
-      Create MaxCardinalitySearch;
-      
-      MaxCardinalitySearch mcs(*_aux_graph, *_aux_capacity);
-      for (NodeIt it(*_aux_graph); it != INVALID; ++it) {
-        _heap_cross_ref->set(it, Heap::PRE_HEAP);
+      typename AuxGraph::template UEdgeMap<Value> _edge_cut(*_aux_graph);
+
+
+      for (NodeIt n(*_aux_graph); n != INVALID; ++n) {
+        _heap_cross_ref->set(n, Heap::PRE_HEAP);
       }
-      mcs.heap(*_heap, *_heap_cross_ref);
 
-      LastTwoMap<Node> last_two_nodes(_node_num);
-      mcs.processedMap(last_two_nodes);
+      std::vector<Node> nodes;
+      nodes.reserve(_node_num);
+      int sep = 0;
 
-      NullMap<Node, Value> cardinality;
-      mcs.cardinalityMap(cardinality);
+      Value alpha = 0;
+      Value emc = std::numeric_limits<Value>::max();
 
-      mcs.run();
+      _heap->push(typename AuxGraph::NodeIt(*_aux_graph), 0);
+      while (!_heap->empty()) {
+        Node node = _heap->top();
+        Value value = _heap->prio();
+        
+        _heap->pop();
+        for (typename AuxGraph::IncEdgeIt e(*_aux_graph, node); 
+             e != INVALID; ++e) {
+          Node tn = _aux_graph->runningNode(e);
+          switch (_heap->state(tn)) {
+          case Heap::PRE_HEAP:
+            _heap->push(tn, (*_aux_capacity)[e]);
+            _edge_cut[e] = (*_heap)[tn];
+            break;
+          case Heap::IN_HEAP:
+            _heap->decrease(tn, (*_aux_capacity)[e] + (*_heap)[tn]);
+            _edge_cut[e] = (*_heap)[tn];
+            break;
+          case Heap::POST_HEAP:
+            break;
+          }
+        }
 
-      Node new_node = _aux_graph->addNode();
+        alpha += (*_aux_cut_value)[node];
+        alpha -= 2 * value;
 
-      typename AuxGraph::template NodeMap<UEdge> edges(*_aux_graph, INVALID);
-      
-      Node first_node = last_two_nodes[0];
-      Node second_node = last_two_nodes[1];
-      
-      _next->set((*_last)[first_node], (*_first)[second_node]);
-      _first->set(new_node, (*_first)[first_node]);
-      _last->set(new_node, (*_last)[second_node]);
-
-      Value current_cut = 0;      
-      for (IncEdgeIt it(*_aux_graph, first_node); it != INVALID; ++it) {
-        Node node = _aux_graph->runningNode(it);
-        current_cut += (*_aux_capacity)[it];
-        if (node == second_node) continue;
-        if (edges[node] == INVALID) {
-          edges[node] = _aux_graph->addEdge(new_node, node);
-          (*_aux_capacity)[edges[node]] = (*_aux_capacity)[it];
-        } else {
-          (*_aux_capacity)[edges[node]] += (*_aux_capacity)[it];          
+        nodes.push_back(node);
+        if (!_heap->empty()) {
+          if (alpha < emc) {
+            emc = alpha;
+            sep = nodes.size();
+          }
         }
       }
 
-      if (_first_node == INVALID || current_cut < _min_cut) {
-        _first_node = (*_first)[first_node];
-        _last_node = (*_last)[first_node];
-        _min_cut = current_cut;
+      if ((int)nodes.size() < _node_num) {
+        _aux_graph->clear();
+        _node_num = 1;
+        _cut.clear();
+        for (int i = 0; i < (int)nodes.size(); ++i) {
+          typename Graph::Node n = (*_first)[nodes[i]];
+          while (n != INVALID) {
+            _cut.push_back(n);
+            n = (*_next)[n];
+          }
+        }
+        _min_cut = 0;
+        return true;
       }
 
-      _aux_graph->erase(first_node);
+      if (emc < _min_cut) {
+        _cut.clear();
+        for (int i = 0; i < sep; ++i) {
+          typename Graph::Node n = (*_first)[nodes[i]];
+          while (n != INVALID) {
+            _cut.push_back(n);
+            n = (*_next)[n];
+          }
+        }
+        _min_cut = emc;
+      }
 
-      for (IncEdgeIt it(*_aux_graph, second_node); it != INVALID; ++it) {
-        Node node = _aux_graph->runningNode(it);
-        if (edges[node] == INVALID) {
-          edges[node] = _aux_graph->addEdge(new_node, node);
-          (*_aux_capacity)[edges[node]] = (*_aux_capacity)[it];
-        } else {
-          (*_aux_capacity)[edges[node]] += (*_aux_capacity)[it];          
+      typedef typename AuxGraph::template NodeMap<int> UfeCr;
+      UfeCr ufecr(*_aux_graph);
+      typedef UnionFindEnum<UfeCr> Ufe; 
+      Ufe ufe(ufecr);
+
+      for (typename AuxGraph::NodeIt n(*_aux_graph); n != INVALID; ++n) {
+        ufe.insert(n);
+      }
+
+      for (typename AuxGraph::UEdgeIt e(*_aux_graph); e != INVALID; ++e) {
+        if (_edge_cut[e] >= emc) {
+          ufe.join(_aux_graph->source(e), _aux_graph->target(e));
         }
       }
-      _aux_graph->erase(second_node);
 
-      --_node_num;
+      if (ufe.size((typename Ufe::ClassIt)(ufe)) == _node_num) {
+        _aux_graph->clear();
+        _node_num = 1;
+        return true;
+      }
+      
+      std::vector<typename AuxGraph::Node> remnodes;
+
+      typename AuxGraph::template NodeMap<UEdge> edges(*_aux_graph, INVALID);
+      for (typename Ufe::ClassIt c(ufe); c != INVALID; ++c) {
+        if (ufe.size(c) == 1) continue;
+        for (typename Ufe::ItemIt r(ufe, c); r != INVALID; ++r) {
+          if ((Node)r == (Node)c) continue;
+          _next->set((*_last)[c], (*_first)[r]);
+          _last->set(c, (*_last)[r]);
+          remnodes.push_back(r);
+          --_node_num;
+        }
+      }
+
+      std::vector<EdgeInfo> addedges;
+      std::vector<UEdge> remedges;
+
+      for (typename AuxGraph::UEdgeIt e(*_aux_graph);
+           e != INVALID; ++e) {
+        Node sn = ufe.find(_aux_graph->source(e));
+        Node tn = ufe.find(_aux_graph->target(e));
+        if ((ufe.size(sn) == 1 && ufe.size(tn) == 1)) {
+          continue;
+        }
+        if (sn == tn) {
+          remedges.push_back(e);
+          continue;
+        }
+        EdgeInfo info;
+        if (sn < tn) {
+          info.source = sn;
+          info.target = tn;
+        } else {
+          info.source = tn;
+          info.target = sn;
+        }
+        info.capacity = (*_aux_capacity)[e];
+        addedges.push_back(info);
+        remedges.push_back(e);
+      }
+
+      for (int i = 0; i < (int)remedges.size(); ++i) {
+        _aux_graph->erase(remedges[i]);
+      }
+
+      sort(addedges.begin(), addedges.end(), EdgeInfoLess());
+
+      {
+        int i = 0;
+        while (i < (int)addedges.size()) {
+          Node sn = addedges[i].source;
+          Node tn = addedges[i].target;
+          Value ec = addedges[i].capacity;
+          ++i;
+          while (i < (int)addedges.size() && 
+                 sn == addedges[i].source && tn == addedges[i].target) {
+            ec += addedges[i].capacity;
+            ++i;
+          }
+          typename AuxGraph::UEdge ne = _aux_graph->addEdge(sn, tn);
+          (*_aux_capacity)[ne] = ec;
+        }
+      }
+
+      for (typename Ufe::ClassIt c(ufe); c != INVALID; ++c) {
+        if (ufe.size(c) == 1) continue;
+        Value cutvalue = 0;
+        for (typename AuxGraph::IncEdgeIt e(*_aux_graph, c);
+             e != INVALID; ++e) {
+          cutvalue += (*_aux_capacity)[e];
+        }
+        
+        (*_aux_cut_value)[c] = cutvalue;
+        
+      }
+
+      for (int i = 0; i < (int)remnodes.size(); ++i) {
+        _aux_graph->erase(remnodes[i]);
+      }
+
       return _node_num == 1;
     }
 
@@ -1317,11 +1500,9 @@
     /// map have been set false previously. 
     template <typename NodeMap>
     Value quickMinCut(NodeMap& nodeMap) const { 
-      for (typename Graph::Node it = _first_node; 
-           it != _last_node; it = (*_next)[it]) {
-             nodeMap.set(it, true);
-           }
-      nodeMap.set(_last_node, true);
+      for (int i = 0; i < (int)_cut.size(); ++i) {
+        nodeMap.set(_cut[i], true);
+      }
       return minCut();
     }
 
@@ -1355,6 +1536,9 @@
 
     ///@}
 
+  private:
+
+
   };
 }